Statistiques
| Révision :

root / src / egrad_gaussian.f90 @ 8

Historique | Voir | Annoter | Télécharger (5,85 ko)

1
subroutine egrad_gaussian(e,geomcart,gradcart)
2

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

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

    
9
  IMPLICIT NONE
10

    
11
  ! Energy (calculated if F300K=.F., else estimated)
12
  REAL(KREAL), INTENT (OUT) :: e
13
  ! Nb degree of freedom
14
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
15
  ! Gradient calculated at Geom geometry
16
  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)
17

    
18
  ! ======================================================================
19

    
20
  character(LCHARS) :: LINE
21

    
22
  logical           :: debug, G03Ok
23
  LOGICAL, SAVE :: first=.true.
24

    
25
  REAL(KREAL) :: Pi
26

    
27
  INTEGER(KINT) :: ItG03
28
  INTEGER(KINT) :: iat, i, n3at
29
  INTEGER(KINT) :: ITmp1, ITmp2
30

    
31
  CHARACTER(132) :: FileIn, FileOut, FileChk
32

    
33
  CHARACTER(3) :: SepIn=' < ', SepOut=' > '
34
  CHARACTER(VLCHARS), SAVE :: RunCommand
35
!  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.
36

    
37
  ! ======================================================================
38

    
39
  LOGICAL, EXTERNAL :: valid
40

    
41
  ! ======================================================================
42

    
43
  Pi=dacos(-1.0d0)
44
  n3at=3*nat
45

    
46
  debug=valid('EGRAD')
47
  if (debug) Call Header("Entering Egrad_gaussian")
48

    
49
  if (first) THEN
50
     First=.FALSE.
51
     call system("uname > Path.tmp")
52
     OPEN(IOTMP,FILE='Path.tmp')
53
     READ(IOTMP,'(A)') Line
54
     CLOSE(IOTMP)
55
     IF (index(Line,'CYGWIN').NE.0) THEN
56
        SepIn='   '
57
        SepOut='   '
58
     END IF
59
  END IF
60

    
61
  gradcart=0.
62
  E=0.
63
  FileIn=Trim(CalcName) // "." // Trim(ISuffix)
64
  FileOut=Trim(CalcName) // "." // Trim(OSuffix)
65
  FileChk=Trim(CalcName) // ".chk"
66
  RunCommand=Trim(ProgExe) // SepIn // Trim(FileIn) // SepOut // Trim(FileOut)
67

    
68
  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)
69

    
70
  ! we create the input file
71

    
72
  OPEN(IOTMP,File=FileIn)
73

    
74
  Current => Gauss_root
75
  DO WHILE (ASSOCIATED(Current%next))
76
     WRITE(IOTMP,'(1X,A)') Trim(current%line)
77
     WRITE(*,'(1X,A)') Trim(current%line)
78
     Current => current%next
79
  END DO
80

    
81
  WRITE(IOTMP,*) 
82

    
83
  Current => Gauss_comment
84
  DO WHILE (ASSOCIATED(Current%next))
85
     WRITE(IOTMP,'(1X,A)') Trim(current%line)
86
     Current => current%next
87
  END DO
88
  !  WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)"),Igeom,NgeomF,Iopt
89

    
90
  WRITE(IOTMP,*) 
91
  WRITE(IOTMP,*) Trim(Gauss_charge)
92

    
93
  DO I=1,Nat
94
     If (renum) THEN
95
        Iat=Order(I)
96
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),GeomCart(Iat,:), TRIM(Gauss_Paste(I))
97
     ELSE
98
        Iat=OrderInv(I)
99
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,:), TRIM(Gauss_Paste(Iat))
100
     END IF
101
  END DO
102

    
103
  WRITE(IOTMP,*) 
104

    
105
  Current => Gauss_End
106
  DO WHILE (ASSOCIATED(Current%next))
107
     WRITE(IOTMP,'(1X,A)') Trim(current%line)
108
     Current => current%next
109
  END DO
110

    
111
  WRITE(IOTMP,*) 
112
  CLOSE(IOTMP)
113

    
114
  IF (debug) THEN
115
     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)
116
     call system('cat ' // Trim(FileIn))
117
  END IF
118

    
119
!!! PFL 15th March 2008 ->
120
! Sometimes for unknown reasons G03 fails.
121
! We thus add one test to check that Gaussian has finisehd properly.
122
! We also check that we do not run an infinite loop...
123

    
124
 G03Ok=.FALSE.
125
 ItG03=0
126

    
127
  DO WHILE ((.NOT.G03Ok).AND.(ItG03<3))
128
   ItG03=ItG03+1
129
  call system(RunCommand)
130
 OPEN(IOTMP,FILE=FileOut, STATUS='old', POSITION='APPEND')
131
 BACKSPACE(IOTMP)
132
 READ(IOTMP,'(A)') Line
133
 IF (LINE(2:19).EQ.'Normal termination')  G03Ok=.TRUE.
134
  END DO
135
  CLOSE(IOTMP)
136
  IF (.NOT.G03Ok) THEN
137
     WRITE(*,*) "Egrad_gaussian Failed..."
138
     WRITE(*,*) "Wrong termination of Gaussian after 3 tries. Check Input and geometry"
139
     STOP
140
  END IF
141

    
142
  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'
143

    
144
  ! Whatever the coordinate system, we read the cartesian gradient
145
  ! and we convert it to Zmat or other if needed
146

    
147
  OPEN(IOTMP,FILE=FileOut, STATUS='old')
148
  ! We first search for the forces
149
  READ(IOTMP,'(A)') LINE
150
  DO WHILE (INDEX(LINE,'Forces (Hartrees/Bohr)')==0)
151
     ! here, we have a problem, because there are so many ways to write E for Gaussian :-)
152
     ! but we do not really need E except if we want to weight the different points
153
     ! on the path... that will be done later :-)
154
     ! And I should use ConvGaussXmol ...
155
! PFL 3rd Jan 2008 ->
156
! I have finally upgraded the energy reading phase so that it should be able to read 
157
! many formats corresponding to HF, DFT, MP2, AM1, PM3, ONIOM with or w/o PCM
158
     
159
! For MP2, energy is after the last =
160
     IF ((Line(2:9)=="E2 =    ")) THEN
161
      Itmp1=Index(LINE,"=",BACK=.TRUE.)+1
162
    READ(LINE(Itmp1:),*) e
163
      END IF
164
! For other methods, it is after the first =
165
     IF ((LINE(2:9)=='Energy= ').OR.(Line(2:9)=="SCF Done").OR.(Line(2:9)=="ONIOM: e") &
166
   .OR.(Line(2:9)==" with al") &
167
         .OR.(Line(2:31)=="Counterpoise: corrected energy")) THEN
168
      Itmp1=Index(LINE,"=")+1
169
    READ(LINE(Itmp1:),*) e
170
  END IF
171
! <- PFL 3 Jan 2008
172
     READ(IOTMP,'(A)') LINE
173
  END DO
174
  if (debug) WRITE(*,*) "GradCart reading - renum=",Renum
175
  READ(IOTMP,'(A)') LINE
176
  READ(IOTMP,'(A)') LINE
177
  DO I=1,Nat
178
     Iat=I
179
     IF (renum) Iat=Order(I)
180
     READ(IOTMP,*)  Itmp1,ITmp2, GradCart(3*Iat-2:3*Iat)
181
     if (debug) WRITE(*,*) "I,Iat,GradCart:",I,Iat,GradCart(3*Iat-2:3*Iat)
182
  END DO
183

    
184
! Gaussian gives the Forces in ua/a0, so we convert it into the 
185
! gradient in ua/Angstrom 
186
  gradCart=-1.0d0*gradCart*Unit
187

    
188
!  if (debug.AND.(COORD.EQ.'ZMAT')) THEN
189
!     DO WHILE (INDEX(LINE,'Internal Coordinate Forces')==0)
190
!        READ(IOTMP,'(A)') LINE
191
!!        WRITE(*,*) 'Pas bon:' // TRIM(LINE)
192
!     END DO
193
!     WRITE(*,*) TRIM(LINE)
194
!     DO I=1,Nat+2
195
!        READ(IOTMP,'(A)') LINE
196
!        WRITE(*,*) TRIM(LINE)
197
!     END DO
198
!  END IF
199

    
200
  CLOSE(IOTMP)
201

    
202
  if (debug) Call header("Egrad_gaussian Over")
203

    
204
  RETURN
205

    
206
!999 CONTINUE
207
!  if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'
208
!  STOP
209
  ! ======================================================================
210
end subroutine egrad_gaussian