Dear siesta users;
I used constr.f routine to implement symmetry in my calculation.
I generate this constr.f by TETR utility. My calculations finishes
successfully,
however I got the same results without any constraints in my calculations.
It always does in siesta-trunk-301 or siesta-trunk-320 version
( the results are the same with and without any constraints).
If I do my calculation in siesta2.0.1 everything works fine
(my calculation with and without constraints are different, I mean forces
and stress.)
Do you have any ideas why new version of siesta can't sees this routine.
I'm not sure about the ncon variable in constr.f
( in siesta2.0.1 there was no need to write ncon variable).
Please help me if you are experienced at writing consrt.f
My constr.f routine is as following:
!
! This file is part of the SIESTA package.
!
! Copyright (c) Fundacion General Universidad Autonoma de Madrid:
! E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal
! and J.M.Soler, 1996-2006.
!
! Use of this software constitutes agreement with the full conditions
! given in the SIESTA license, as signed by all legitimate users.
!
c $Id: constr.f,v 1.6 2003/06/23 09:46:16 ordejon Exp $
subroutine constr( cell, na, isa, amass, xa, stress, fa, ntcon )
c *****************************************************************
c User-written routine to implement specific geometric constraints,
c by orthogonalizing the forces and stress to undesired changes.
c Arguments:
c real*8 cell(3,3) : input lattice vectors (Bohr)
c integer na : input number of atoms
c integer isa(na) : input species indexes
c real*8 amass(na) : input atomic masses
c real*8 xa(3,na) : input atomic cartesian coordinates (Bohr)
c real*8 stress( 3,3) : input/output stress tensor (Ry/Bohr**3)
c real*8 fa(3,na) : input/output atomic forces (Ry/Bohr)
c integer ntcon : total number of positions constr. imposed
c *****************************************************************
implicit none
integer na, isa(na), ntcon
double precision amass(na), cell(3,3), fa(3,na),
. stress(3,3), xa(3,na)
c Write here your problem-specific code.
c..operations: 1 5 11 31 13 28 16 15 32 14 47 17
44 7 33 30 45 46 12 42 18 29 9 43
c UNIT C2Z C3D S4Z3 C3A S4Z -C3C -C3D S4X3 C3B P6 -C3A
P3 C2X S4Y3 S4Y P4 P5 C3C P1 -C3B S4X C2Y P2
c=======> orbit= 1, 1st atom = 1 <=======
c...local: C2Z C3D S4Z3 C3A S4Z -C3C -C3D S4X3 C3B P6
c...local: -C3A P3 C2X S4Y3 S4Y P4 P5 C3C P1 -C3B
c...local: S4X C2Y P2
c type=I, eigenvector = 0.00000 0.00000 1.00000
fa(1,1)=0.0
fa(2,1)=0.0
fa(3,1)=0.0
c=======> orbit= 2, 1st atom = 2 <=======
c...local: C2Z S4Z3 S4Z C2X P1 C2Y P2
c type=I, eigenvector = 0.00000 0.00000 1.00000
fa(1,2)=0.0
fa(2,2)=0.0
fa(3,2)=0.0
c=======> orbit= 3, 1st atom = 5 <=======
c...local: -C3C P3 P5 C3C P1
c type=C, eigenvector = -0.57735 -0.57735 -0.57735
fa(2,5)= 1.00000*fa(1,5)
fa(3,5)= 1.00000*fa(1,5)
c....... force on Ga atom 3 via atom 2 and operation C3D
c fa(1,3)=C(1,1,11)*fa(1,2)+C(1,2,11)*fa(2,2)+C(1,3,11)*fa(3,2)
c fa(2,3)=C(2,1,11)*fa(1,2)+C(2,2,11)*fa(2,2)+C(2,3,11)*fa(3,2)
c fa(3,3)=C(3,1,11)*fa(1,2)+C(3,2,11)*fa(2,2)+C(3,3,11)*fa(3,2)
fa(1,3)=-1.00000*fa(2,2)
fa(2,3)=+1.00000*fa(3,2)
fa(3,3)=-1.00000*fa(1,2)
c
c....... force on Ga atom 4 via atom 2 and operation -C3D
c fa(1,4)=C(1,1,15)*fa(1,2)+C(1,2,15)*fa(2,2)+C(1,3,15)*fa(3,2)
c fa(2,4)=C(2,1,15)*fa(1,2)+C(2,2,15)*fa(2,2)+C(2,3,15)*fa(3,2)
c fa(3,4)=C(3,1,15)*fa(1,2)+C(3,2,15)*fa(2,2)+C(3,3,15)*fa(3,2)
fa(1,4)=-1.00000*fa(3,2)
fa(2,4)=-1.00000*fa(1,2)
fa(3,4)=+1.00000*fa(2,2)
c
c....... force on As atom 6 via atom 5 and operation C2Z
c fa(1,6)=C(1,1,5)*fa(1,5)+C(1,2,5)*fa(2,5)+C(1,3,5)*fa(3,5)
c fa(2,6)=C(2,1,5)*fa(1,5)+C(2,2,5)*fa(2,5)+C(2,3,5)*fa(3,5)
c fa(3,6)=C(3,1,5)*fa(1,5)+C(3,2,5)*fa(2,5)+C(3,3,5)*fa(3,5)
fa(1,6)=-1.00000*fa(1,5)
fa(2,6)=-1.00000*fa(2,5)
fa(3,6)=+1.00000*fa(3,5)
c
c....... force on As atom 7 via atom 5 and operation C3D
c fa(1,7)=C(1,1,11)*fa(1,5)+C(1,2,11)*fa(2,5)+C(1,3,11)*fa(3,5)
c fa(2,7)=C(2,1,11)*fa(1,5)+C(2,2,11)*fa(2,5)+C(2,3,11)*fa(3,5)
c fa(3,7)=C(3,1,11)*fa(1,5)+C(3,2,11)*fa(2,5)+C(3,3,11)*fa(3,5)
fa(1,7)=-1.00000*fa(2,5)
fa(2,7)=+1.00000*fa(3,5)
fa(3,7)=-1.00000*fa(1,5)
c
c....... force on As atom 8 via atom 5 and operation C3A
c fa(1,8)=C(1,1,13)*fa(1,5)+C(1,2,13)*fa(2,5)+C(1,3,13)*fa(3,5)
c fa(2,8)=C(2,1,13)*fa(1,5)+C(2,2,13)*fa(2,5)+C(2,3,13)*fa(3,5)
c fa(3,8)=C(3,1,13)*fa(1,5)+C(3,2,13)*fa(2,5)+C(3,3,13)*fa(3,5)
fa(1,8)=+1.00000*fa(2,5)
fa(2,8)=-1.00000*fa(3,5)
fa(3,8)=-1.00000*fa(1,5)
c
c... STRESS tensor: crystal class Td , dimension = 3D
stress(1,3)=0.0
stress(2,3)=0.0
stress(1,2)=0.0
stress(2,2)=stress(1,1)
stress(3,3)=stress(1,1)
stress(2,1)=stress(1,2)
stress(3,1)=stress(1,3)
stress(3,2)=stress(2,3)
ntcon=9
end
Many thanks in advance for any suggestions.
Best regards,
Magda Birowska
