Statistics
| Revision:

root / src / egrad_siesta.f90 @ 7

History | View | Annotate | Download (5.4 kB)

1 5 pfleura2
subroutine egrad_siesta(e,geomcart,gradcart)
2 5 pfleura2
3 5 pfleura2
  ! This routines calculates the energy and the gradient of
4 5 pfleura2
  ! a molecule, using Siesta
5 5 pfleura2
6 5 pfleura2
  use Path_module, only : Nat, renum,Order,OrderInv, Coord, ProgExe
7 5 pfleura2
  use Io_module
8 5 pfleura2
9 5 pfleura2
  IMPLICIT NONE
10 5 pfleura2
11 5 pfleura2
  INTERFACE
12 5 pfleura2
     function valid(string) result (isValid)
13 5 pfleura2
       CHARACTER(*), intent(in) :: string
14 5 pfleura2
       logical                  :: isValid
15 5 pfleura2
     END function VALID
16 5 pfleura2
17 5 pfleura2
18 5 pfleura2
    SUBROUTINE die(routine, msg, file, line, unit)
19 5 pfleura2
20 5 pfleura2
      Use VarTypes
21 5 pfleura2
      Use io_module
22 5 pfleura2
23 5 pfleura2
      implicit none
24 5 pfleura2
!--------------------------------------------------------------- Input Variables
25 5 pfleura2
      character(len=*), intent(in)           :: routine, msg
26 5 pfleura2
      character(len=*), intent(in), optional :: file
27 5 pfleura2
      integer(KINT), intent(in), optional      :: line, unit
28 5 pfleura2
29 5 pfleura2
    END SUBROUTINE die
30 5 pfleura2
31 5 pfleura2
    SUBROUTINE Warning(routine, msg, file, line, unit)
32 5 pfleura2
33 5 pfleura2
      Use VarTypes
34 5 pfleura2
      Use io_module
35 5 pfleura2
36 5 pfleura2
      implicit none
37 5 pfleura2
38 5 pfleura2
      character(len=*), intent(in)           :: routine, msg
39 5 pfleura2
      character(len=*), intent(in), optional :: file
40 5 pfleura2
      integer(KINT), intent(in), optional      :: line, unit
41 5 pfleura2
42 5 pfleura2
    END SUBROUTINE Warning
43 5 pfleura2
44 5 pfleura2
45 5 pfleura2
 SUBROUTINE WriteList(Input,Unit)
46 5 pfleura2
47 5 pfleura2
! This routine reads an input template for Siesta
48 5 pfleura2
49 5 pfleura2
  use VarTypes
50 5 pfleura2
  use Io_module
51 5 pfleura2
52 5 pfleura2
  IMPLICIT NONE
53 5 pfleura2
54 5 pfleura2
! Input
55 5 pfleura2
      TYPE (Input_line), POINTER, INTENT(IN) :: Input
56 5 pfleura2
      INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit
57 5 pfleura2
58 5 pfleura2
    END SUBROUTINE WriteList
59 5 pfleura2
  END INTERFACE
60 5 pfleura2
61 5 pfleura2
62 5 pfleura2
  ! Energy
63 5 pfleura2
  REAL(KREAL), INTENT (OUT) :: e
64 5 pfleura2
  ! Geometry in cartesian coordinates
65 5 pfleura2
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
66 5 pfleura2
  ! Gradient calculated at Geom geometry
67 5 pfleura2
  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)
68 5 pfleura2
69 5 pfleura2
  ! ======================================================================
70 5 pfleura2
71 5 pfleura2
  character(LCHARS) :: LINE
72 5 pfleura2
73 5 pfleura2
  logical           :: debug, TChk
74 5 pfleura2
  LOGICAL, SAVE :: first=.true.
75 5 pfleura2
76 5 pfleura2
  REAL(KREAL) :: Pi
77 5 pfleura2
78 5 pfleura2
  INTEGER(KINT) :: iat, i, n3at,Idx
79 5 pfleura2
80 5 pfleura2
  CHARACTER(LCHARS) :: FileIn, FileOut, FileForces,FileE
81 5 pfleura2
82 5 pfleura2
  CHARACTER(3) :: SepIn=' < ', SepOut=' > '
83 5 pfleura2
  CHARACTER(VLCHARS), SAVE :: RunCommand
84 5 pfleura2
85 5 pfleura2
  ! ======================================================================
86 5 pfleura2
87 5 pfleura2
  Pi=dacos(-1.0d0)
88 5 pfleura2
  n3at=3*nat
89 5 pfleura2
90 5 pfleura2
  debug=valid('EGRAD').OR.valid('egrad_siesta')
91 5 pfleura2
92 5 pfleura2
  if (debug) Call Header("Entering Egrad_siesta")
93 5 pfleura2
94 5 pfleura2
  if (first) THEN
95 5 pfleura2
     First=.FALSE.
96 5 pfleura2
     call system("uname > Path.tmp")
97 5 pfleura2
     OPEN(IOTMP,FILE='Path.tmp')
98 5 pfleura2
     READ(IOTMP,'(A)') Line
99 5 pfleura2
     CLOSE(IOTMP)
100 5 pfleura2
     IF (index(Line,'CYGWIN').NE.0) THEN
101 5 pfleura2
        SepIn='   '
102 5 pfleura2
        SepOut='   '
103 5 pfleura2
     END IF
104 5 pfleura2
  END IF
105 5 pfleura2
106 5 pfleura2
  gradcart=0.
107 5 pfleura2
  E=0.
108 5 pfleura2
  FileIn=Trim(CalcName) // "." // Trim(ISuffix)
109 5 pfleura2
  FileOut=Trim(CalcName) // "." // Trim(OSuffix)
110 5 pfleura2
  FileForces=Trim(Siesta_Label) // ".FA"
111 5 pfleura2
  FileE='BASIS_ENTHALPY'
112 5 pfleura2
113 5 pfleura2
  RunCommand=Trim(ProgExe) // SepIn // Trim(FileIn) // SepOut // Trim(FileOut)
114 5 pfleura2
115 5 pfleura2
  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)
116 5 pfleura2
117 5 pfleura2
  ! we create the input file
118 5 pfleura2
119 5 pfleura2
  OPEN(IOTMP,File=FileIn)
120 5 pfleura2
121 5 pfleura2
  Call WriteList(Siesta_input,Unit=IOTMP)
122 5 pfleura2
123 5 pfleura2
  WRITE(IOTMP,*)
124 5 pfleura2
125 5 pfleura2
  WRITE(IOTMP,'(1X,A)')  '%block AtomicCoordinatesAndAtomicSpecies'
126 5 pfleura2
127 5 pfleura2
  DO I=1,Nat
128 5 pfleura2
     If (renum) THEN
129 5 pfleura2
        Iat=Order(I)
130 5 pfleura2
        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(Iat,:)/Siesta_Unit_Read, IdxSpecies(Iat),TRIM(Siesta_Paste(I))
131 5 pfleura2
     ELSE
132 5 pfleura2
        Iat=OrderInv(I)
133 5 pfleura2
        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(I,:)/Siesta_Unit_Read, IdxSpecies(Iat), TRIM(Siesta_Paste(Iat))
134 5 pfleura2
     END IF
135 5 pfleura2
  END DO
136 5 pfleura2
137 5 pfleura2
  WRITE(IOTMP,'(1X,A)')  '%endblock AtomicCoordinatesAndAtomicSpecies'
138 5 pfleura2
  WRITE(IOTMP,*)
139 5 pfleura2
140 5 pfleura2
  CLOSE(IOTMP)
141 5 pfleura2
142 5 pfleura2
  IF (debug) THEN
143 5 pfleura2
     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)
144 5 pfleura2
     call system('cat ' // Trim(FileIn))
145 5 pfleura2
  END IF
146 5 pfleura2
147 5 pfleura2
  call system(RunCommand)
148 5 pfleura2
149 5 pfleura2
  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'
150 5 pfleura2
151 5 pfleura2
  ! Whatever the coordinate system, we read the cartesian gradient
152 5 pfleura2
  ! and we convert it to Zmat or other if needed
153 5 pfleura2
154 5 pfleura2
! in Siesta, the forces are in the file 'Label'.FA so this is easy to read :)
155 5 pfleura2
  OPEN(IOTMP,FILE=FileForces, STATUS='old')
156 5 pfleura2
! firt line is the atom number
157 5 pfleura2
  READ(IOTMP,'(A)') LINE
158 5 pfleura2
159 5 pfleura2
  if (debug) WRITE(*,*) "GradCart reading - renum=",Renum
160 5 pfleura2
  DO I=1,Nat
161 5 pfleura2
     Iat=I
162 5 pfleura2
     IF (renum) Iat=Order(I)
163 5 pfleura2
     READ(IOTMP,'(6X,3(F12.6))')  GradCart(3*Iat-2:3*Iat)
164 5 pfleura2
     if (debug) WRITE(*,*) "I,Iat,GradCart:",I,Iat,GradCart(3*Iat-2:3*Iat)
165 5 pfleura2
  END DO
166 5 pfleura2
167 5 pfleura2
! Siesta gives the Forces in eV/Ang, so we convert it into the
168 5 pfleura2
! gradient in au/Angstrom
169 5 pfleura2
  gradCart=-1.0d0*ev2au*gradCart
170 5 pfleura2
171 5 pfleura2
  CLOSE(IOTMP)
172 5 pfleura2
173 5 pfleura2
! We now have to read the energy
174 5 pfleura2
175 5 pfleura2
  INQUIRE(File=FileE,EXIST=TChk)
176 5 pfleura2
  If (Tchk) THEN
177 5 pfleura2
! We have a BASIS_ENTHALPY file. Great.
178 5 pfleura2
! Energy is the last number on the second line
179 5 pfleura2
     OPEN(IOTMP,File=FileE,Status='Old')
180 5 pfleura2
     READ(IOTMP,*) Line
181 5 pfleura2
     READ(IOTMP,'(A)') Line
182 5 pfleura2
     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)
183 5 pfleura2
     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)
184 5 pfleura2
     Line=Line(Idx+1:)
185 5 pfleura2
     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)
186 5 pfleura2
     Read(Line,*) E
187 5 pfleura2
  ELSE
188 5 pfleura2
! We do not have this file... we have to look for it in the .out file
189 5 pfleura2
     FileE=FileOut
190 5 pfleura2
     OPEN(IOTMP,File=FileE,Status='Old')
191 5 pfleura2
! We search for  "siesta:         Total ="
192 5 pfleura2
     Line=""
193 5 pfleura2
     DO WHILE (INDEX(Line,"siesta:         Total =")/=1)
194 5 pfleura2
        READ(IOTMP,*) Line
195 5 pfleura2
     END DO
196 5 pfleura2
197 5 pfleura2
     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)
198 5 pfleura2
     Line=Line(Idx+1:)
199 5 pfleura2
     Read(Line,*) E
200 5 pfleura2
  END IF
201 5 pfleura2
  CLOSE (IOTMP)
202 5 pfleura2
203 5 pfleura2
  WRITE(*,*) 'E=',E
204 5 pfleura2
  WRITE(*,*) 'Gradcart='
205 5 pfleura2
  DO I=1,Nat
206 5 pfleura2
     WRITE(*,'(1X,3(1X,F15.8))') GradCart(3*I-2:3*I)
207 5 pfleura2
  END DO
208 5 pfleura2
209 5 pfleura2
!  Call Die('Egrad_Siesta','Ah ah.',UNIT=IOOUT)
210 5 pfleura2
  if (debug) Call header("Egrad_siesta Over")
211 5 pfleura2
212 5 pfleura2
  RETURN
213 5 pfleura2
214 5 pfleura2
  ! ======================================================================
215 5 pfleura2
end subroutine egrad_siesta