Statistics
| Revision:

root / src / egrad_siesta.f90 @ 7

History | View | Annotate | Download (5.4 kB)

1
subroutine egrad_siesta(e,geomcart,gradcart)
2

    
3
  ! This routines calculates the energy and the gradient of 
4
  ! a molecule, using Siesta 
5

    
6
  use Path_module, only : Nat, renum,Order,OrderInv, Coord, ProgExe
7
  use Io_module
8

    
9
  IMPLICIT NONE
10

    
11
  INTERFACE
12
     function valid(string) result (isValid)
13
       CHARACTER(*), intent(in) :: string
14
       logical                  :: isValid
15
     END function VALID
16

    
17

    
18
    SUBROUTINE die(routine, msg, file, line, unit)
19

    
20
      Use VarTypes
21
      Use io_module
22

    
23
      implicit none
24
!--------------------------------------------------------------- Input Variables
25
      character(len=*), intent(in)           :: routine, msg
26
      character(len=*), intent(in), optional :: file
27
      integer(KINT), intent(in), optional      :: line, unit
28

    
29
    END SUBROUTINE die
30

    
31
    SUBROUTINE Warning(routine, msg, file, line, unit)
32

    
33
      Use VarTypes
34
      Use io_module
35

    
36
      implicit none
37

    
38
      character(len=*), intent(in)           :: routine, msg
39
      character(len=*), intent(in), optional :: file
40
      integer(KINT), intent(in), optional      :: line, unit
41

    
42
    END SUBROUTINE Warning
43

    
44

    
45
 SUBROUTINE WriteList(Input,Unit)
46

    
47
! This routine reads an input template for Siesta
48

    
49
  use VarTypes
50
  use Io_module
51

    
52
  IMPLICIT NONE
53

    
54
! Input
55
      TYPE (Input_line), POINTER, INTENT(IN) :: Input
56
      INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit
57

    
58
    END SUBROUTINE WriteList
59
  END INTERFACE
60

    
61

    
62
  ! Energy 
63
  REAL(KREAL), INTENT (OUT) :: e
64
  ! Geometry in cartesian coordinates
65
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
66
  ! Gradient calculated at Geom geometry
67
  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)
68

    
69
  ! ======================================================================
70

    
71
  character(LCHARS) :: LINE
72

    
73
  logical           :: debug, TChk
74
  LOGICAL, SAVE :: first=.true.
75

    
76
  REAL(KREAL) :: Pi
77

    
78
  INTEGER(KINT) :: iat, i, n3at,Idx
79

    
80
  CHARACTER(LCHARS) :: FileIn, FileOut, FileForces,FileE
81

    
82
  CHARACTER(3) :: SepIn=' < ', SepOut=' > '
83
  CHARACTER(VLCHARS), SAVE :: RunCommand
84

    
85
  ! ======================================================================
86

    
87
  Pi=dacos(-1.0d0)
88
  n3at=3*nat
89

    
90
  debug=valid('EGRAD').OR.valid('egrad_siesta')
91

    
92
  if (debug) Call Header("Entering Egrad_siesta")
93

    
94
  if (first) THEN
95
     First=.FALSE.
96
     call system("uname > Path.tmp")
97
     OPEN(IOTMP,FILE='Path.tmp')
98
     READ(IOTMP,'(A)') Line
99
     CLOSE(IOTMP)
100
     IF (index(Line,'CYGWIN').NE.0) THEN
101
        SepIn='   '
102
        SepOut='   '
103
     END IF
104
  END IF
105

    
106
  gradcart=0.
107
  E=0.
108
  FileIn=Trim(CalcName) // "." // Trim(ISuffix)
109
  FileOut=Trim(CalcName) // "." // Trim(OSuffix)
110
  FileForces=Trim(Siesta_Label) // ".FA"
111
  FileE='BASIS_ENTHALPY'
112

    
113
  RunCommand=Trim(ProgExe) // SepIn // Trim(FileIn) // SepOut // Trim(FileOut)
114

    
115
  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)
116

    
117
  ! we create the input file
118

    
119
  OPEN(IOTMP,File=FileIn)
120

    
121
  Call WriteList(Siesta_input,Unit=IOTMP)
122

    
123
  WRITE(IOTMP,*) 
124

    
125
  WRITE(IOTMP,'(1X,A)')  '%block AtomicCoordinatesAndAtomicSpecies' 
126

    
127
  DO I=1,Nat
128
     If (renum) THEN
129
        Iat=Order(I)
130
        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(Iat,:)/Siesta_Unit_Read, IdxSpecies(Iat),TRIM(Siesta_Paste(I))
131
     ELSE
132
        Iat=OrderInv(I)
133
        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(I,:)/Siesta_Unit_Read, IdxSpecies(Iat), TRIM(Siesta_Paste(Iat))
134
     END IF
135
  END DO
136

    
137
  WRITE(IOTMP,'(1X,A)')  '%endblock AtomicCoordinatesAndAtomicSpecies'
138
  WRITE(IOTMP,*) 
139

    
140
  CLOSE(IOTMP)
141

    
142
  IF (debug) THEN
143
     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)
144
     call system('cat ' // Trim(FileIn))
145
  END IF
146

    
147
  call system(RunCommand)
148

    
149
  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'
150

    
151
  ! Whatever the coordinate system, we read the cartesian gradient
152
  ! and we convert it to Zmat or other if needed
153

    
154
! in Siesta, the forces are in the file 'Label'.FA so this is easy to read :)
155
  OPEN(IOTMP,FILE=FileForces, STATUS='old')
156
! firt line is the atom number
157
  READ(IOTMP,'(A)') LINE
158

    
159
  if (debug) WRITE(*,*) "GradCart reading - renum=",Renum
160
  DO I=1,Nat
161
     Iat=I
162
     IF (renum) Iat=Order(I)
163
     READ(IOTMP,'(6X,3(F12.6))')  GradCart(3*Iat-2:3*Iat)
164
     if (debug) WRITE(*,*) "I,Iat,GradCart:",I,Iat,GradCart(3*Iat-2:3*Iat)
165
  END DO
166

    
167
! Siesta gives the Forces in eV/Ang, so we convert it into the 
168
! gradient in au/Angstrom 
169
  gradCart=-1.0d0*ev2au*gradCart
170

    
171
  CLOSE(IOTMP)
172

    
173
! We now have to read the energy
174

    
175
  INQUIRE(File=FileE,EXIST=TChk)
176
  If (Tchk) THEN
177
! We have a BASIS_ENTHALPY file. Great. 
178
! Energy is the last number on the second line
179
     OPEN(IOTMP,File=FileE,Status='Old')
180
     READ(IOTMP,*) Line
181
     READ(IOTMP,'(A)') Line
182
     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)
183
     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)
184
     Line=Line(Idx+1:)
185
     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)
186
     Read(Line,*) E
187
  ELSE
188
! We do not have this file... we have to look for it in the .out file
189
     FileE=FileOut
190
     OPEN(IOTMP,File=FileE,Status='Old')
191
! We search for  "siesta:         Total ="
192
     Line=""
193
     DO WHILE (INDEX(Line,"siesta:         Total =")/=1)
194
        READ(IOTMP,*) Line
195
     END DO
196

    
197
     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)
198
     Line=Line(Idx+1:)
199
     Read(Line,*) E
200
  END IF
201
  CLOSE (IOTMP)
202

    
203
  WRITE(*,*) 'E=',E
204
  WRITE(*,*) 'Gradcart='
205
  DO I=1,Nat
206
     WRITE(*,'(1X,3(1X,F15.8))') GradCart(3*I-2:3*I)
207
  END DO
208

    
209
!  Call Die('Egrad_Siesta','Ah ah.',UNIT=IOOUT)
210
  if (debug) Call header("Egrad_siesta Over")
211

    
212
  RETURN
213

    
214
  ! ======================================================================
215
end subroutine egrad_siesta