Statistiques
| Révision :

root / src / EgradPath.f90 @ 9

Historique | Voir | Annoter | Télécharger (27,75 ko)

1
SUBROUTINE EGradPath(Iopt)
2

    
3
  ! This subroutine calculates the energy+gradient for all points of
4
  ! the path
5

    
6
  use Path_module, only : ProgExe, NCoord, PROG, Nat, NGeomF, Coord, IntCoordF, &
7
       IndZmat, XyzGeomF, Atome, Order, OrderInv, Latr, FrozAtoms, &
8
       Grad, Energies, renum, dzdc, massat, Pi, NCART, OptReac, &
9
       OptProd,  AtName, Unit, CalcEReac, CalcEProd
10
  ! GeomOld_all(NGeomF,NCoord), BTransInv (3*Nat,NCoord): used for Baker case
11
  use Io_module
12

    
13
  IMPLICIT NONE
14

    
15
  INTEGER(KINT), INTENT(IN) :: Iopt
16

    
17
  INTEGER(KINT) :: IGeom,IGeom0,IGeomF
18
  CHARACTER(LCHARS) :: Line,FileInp,RunCommand
19
  CHARACTER(SCHARS) :: Line2
20

    
21
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:) ! NCoord
22
  REAL(KREAL), ALLOCATABLE :: GradCart(:) ! 3*Nat
23
  REAL(KREAL), ALLOCATABLE :: PathCart(:,:) ! (3*Nat,NGeomF)
24
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
25
  REAL(KREAL), ALLOCATABLE :: x_k(:), y_k(:), z_k(:)
26
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3)
27
  REAL(KREAL) :: E
28
  LOGICAL :: Debug
29
  LOGICAL, SAVE :: First=.TRUE.
30

    
31
  INTEGER(KINT) :: I, Iat, IBeg, Idx
32
  REAL(KREAL) :: RTmp1, RTmp2, RTmp3
33
  INTEGER(KINT) :: Itmp1, ITmp2
34

    
35
  INTERFACE
36
     function valid(string) result (isValid)
37
       CHARACTER(*), intent(in) :: string
38
       logical                  :: isValid
39
     END function VALID
40

    
41
     subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old)
42

    
43
  ! This routines calculates the energy E and the gradient Grad of 
44
  ! a molecule with Geometry Geom (may be in internal coordinates),
45
  ! using for now, either Gaussian or Ext, more general later.
46

    
47
  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
48
       prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &
49
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
50
       , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
51
      , Order,OrderInv, XPrimitiveF
52

    
53
  ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
54
  ! allocated in Path.f90
55

    
56
  use Io_module
57

    
58
  ! Energy (calculated if F300K=.F., else estimated)
59
  REAL(KREAL), INTENT (OUT) :: E
60
  ! NCoord: Number of the degrees of freedom
61
  ! IGeom: index of the geometry.
62
  INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
63
  ! Geometry at which gradient is calculated (cf Factual also):
64
  REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
65
  ! Gradient calculated at Geom geometry:
66
  REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
67
  ! Cartesian geometry corresponding to (Internal Geometry) Geom:
68
  REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)
69
!!! Optional, just for geometry optimization with Baker coordinates
70
  REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
71
  REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
72
! FOptGeom is a flag indicating if we are doing a geom optimization
73
! it can be omitted so that we use a local flag: Flag_Opt_Geom
74
  LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
75
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
76
  LOGICAL  :: Flag_Opt_Geom
77

    
78
 END subroutine Egrad
79

    
80
  END INTERFACE
81

    
82
  debug=valid('EGradPath')
83

    
84
  if (debug) Call header("Entering EgradPath")
85

    
86
  ALLOCATE(GeomTmp(NCoord))
87
  ALLOCATE(x(Nat),y(Nat),z(Nat))
88
  ALLOCATE(x_k(Nat),y_k(Nat),z_k(Nat))
89

    
90

    
91
  IF (RunMode=="PARA") THEN ! matches at L315
92

    
93
     ALLOCATE(GradCart(3*Nat),PathCart(3*Nat,NGeomF))
94

    
95

    
96
     SELECT CASE(Prog)
97
     CASE ("VASP")
98
        ! We will use the NEB routine of VASP to get all forces at once on Para8 queue
99

    
100
        ! First, we create all the POSCAR into the 0x directories
101

    
102
        ! With RunMode="Para" one cannot optimize reactants or products (due to Vasp constraints)
103
        ! so this has been done beforehand :-)
104
        IGeom0=2
105
        IGeomF=NGeomF-1
106
        PathCart=0.d0
107

    
108
        IF (First) THEN
109
           ! For the first iteration, we do write a POSCAR into 00 and 0Final
110
           ! Vasp needs it !
111
           First=.FALSE.
112
           IGeom0=1
113
           IGeomF=NGeomF
114
        END IF
115

    
116
        Line2="00"
117
        DO IGeom=IGeom0,IGeomF
118

    
119
           ! NOTE THAT HERE IGeom0 (IGeomF) IS NOT 1 (NGeomF) FOR NON-FIRST OPTIMIZATION.
120
           ! THIS WILL HAVE EFFECT ON XyzGeomF UPDATE IN BAKER CASE.
121
           WRITE(Line,'(I5)') IGeom-1
122
           Idx=2-Len_TRIM(ADJUSTL(Line))
123
           FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/POSCAR"
124
           if (debug) WRITE(*,*) "Creating ",TRIM(FileInp)
125
           OPEN(IOTMP,File=TRIM(FileInp))
126

    
127
           WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt
128
           if (debug) WRITE(*,*) Line
129
           SELECT CASE(Coord)
130
           CASE ('ZMAT')
131
              GeomTmp=IntCoordF(IGeom,:)
132
              Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
133
           CASE ('BAKER')
134
              GeomTmp=IntCoordF(IGeom,:)
135
              !
136
              IF (IOpt .EQ. 0) THEN
137
                 ! EgradPath(...) is called only after the call of PathCreate(...)
138
                 ! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always 
139
                 ! true if IOpt==0), THEN ExtraPol_baker(...) is called.
140
                 ! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting
141
                 ! the interpolated internal coordinates into cartesian coordinates.
142
                 ! Thus here we don't need to reconvert the internal coordinates again
143
                 ! to the cartesian coordinates. 
144
                 DO I=1,Nat
145
                    PathCart(3*I-2,IGeom) = XyzGeomF(IGeom,1,I)
146
                    PathCart(3*I-1,IGeom) = XyzGeomF(IGeom,2,I)
147
                    PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I)
148
                 END DO
149
              ELSE
150
                 ! XyzGeomF(NGeomF,3,Nat).
151
                 x_k(:) = XyzGeomF(IGeom,1,:)
152
                 y_k(:) = XyzGeomF(IGeom,2,:)
153
                 z_k(:) = XyzGeomF(IGeom,3,:)
154

    
155
                 !Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), &
156
                 !                               x_k,y_k,z_k,x,y,z)
157
                 ! PathCart(3*Nat,NGeomF)
158
                 DO I=1,Nat
159
                    PathCart(3*I-2,IGeom) = x(I)
160
                    PathCart(3*I-1,IGeom) = y(I)
161
                    PathCart(3*I,IGeom) = z(I)
162
                 END DO
163

    
164
                 XyzGeomF(IGeom,1,:)=x(:)
165
                 XyzGeomF(IGeom,2,:)=y(:)
166
                 XyzGeomF(IGeom,3,:)=z(:)
167
              END IF
168
           CASE ('MIXED')
169
              GeomTmp=IntCoordF(IGeom,:)
170
              Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
171
           CASE ('CART','HYBRID')
172
!!! CAUTION : PFL 29.AUG.2008 ->
173
              ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
174
              ! This should be made  consistent !!!!!!!!!!!!!!!
175
              GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
176
!!! <- CAUTION : PFL 29.AUG.2008 
177
            PathCart(:,IGeom)=GeomTmp
178
          CASE DEFAULT
179
            WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L134. STOP"
180
            STOP
181
        END SELECT
182

    
183
           IF (debug) THEN
184
              WRITE(*,*) "Calling PrintGeomVasp for geom",Igeom
185
              Call PrintGeom(Line,Nat,NCoord,PathCart(1,Igeom),Coord,IoOut,Atome,Order,OrderInv,IndZmat)
186
           END IF
187
           Call PrintGeomVasp(Line,PathCart(1,IGeom),Latr,FrozAtoms,IoTmp)
188

    
189
           !           IF (debug) THEN 
190
           !              WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD), ' Igeom=',IGeom
191
           !              WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
192
           !              WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
193
           !              DO Iat=1,Nat
194
           !                 WRITE(*,*) XyzGeomF(IGeom,1:3,Iat)
195
           !              END DO
196
           !           END IF
197

    
198
           IF (debug) THEN 
199
              WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD),' Igeom=',IGeom,' Iopt=',Iopt
200
              WRITE(*,*) GeomTmp(1:min(NCoord,9))
201
           END IF
202

    
203
           CLOSE(IOTMP)
204
        END DO ! matches DO IGeom=IGeom0,IGeomF
205

    
206
!!!!!!!!!!!!!!!!!!!
207
!
208
        ! We calculate the energies and gradients
209
!
210

    
211
        RunCommand=Trim(ProgExe)
212
        if (debug) WRITE(*,*) "Launching:", TRIM(RunCommand)
213
        call system(RunCommand)
214
        if (debug) WRITE(*,*) "Back from:", TRIM(RunCommand)
215

    
216
!!!!!!!!!!!!!
217
!
218
        ! We read the gradients and energies
219
!
220
        IGeom0=2
221
        IGeomF=NGeomF-1
222
        Grad=0.d0
223

    
224
        Line2="00"
225
        DO IGeom=IGeom0,IGeomF
226
           WRITE(Line,'(I5)') IGeom-1
227
           Idx=2-Len_TRIM(ADJUSTL(Line))
228
           FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/OUTCAR"
229
           if (debug) WRITE(*,*) "Reading E and forces from ", TRIM(FileInp)
230
           OPEN(IOTMP,File=TRIM(FileInp))
231

    
232
           ! We first search for the forces
233
           READ(IOTMP,'(A)',END=999,ERR=999) LINE
234
           ! We search for the  last part of the OUTCAR file, after wavefunction is converged
235
           DO WHILE (INDEX(LINE,'EDIFF is reached')==0)
236
              READ(IOTMP,'(A)',END=999,ERR=999) LINE
237
           END DO
238
           READ(IOTMP,'(A)',END=999,ERR=999) LINE
239
           DO WHILE (INDEX(LINE,'TOTEN')==0)
240
              READ(IOTMP,'(A)',END=999,ERR=999) LINE
241
           END DO
242
           Line=Line(26:)
243
           READ(LINE,*) E
244
           if (debug) WRITE(*,'(1X,A,I5,A,F20.6)') "DBG Egrad_VASP E(",Igeom,")=",E
245
           Energies(Igeom)=E
246

    
247
           ! We search for the forces
248
           DO WHILE (INDEX(LINE,'TOTAL-FORCE')==0)
249
              READ(IOTMP,'(A)',END=999,ERR=999) LINE
250
           END DO
251
           READ(IOTMP,'(A)',END=999,ERR=999) LINE
252
           DO I=1,Nat
253
              Iat=I
254
              IF (renum) Iat=Order(I)
255
              READ(IOTMP,*,END=999,ERR=999)  Rtmp1,RTmp2,RTmp3,GradCart(3*Iat-2:3*Iat)
256
           END DO
257

    
258
           Close(IoTmp)
259

    
260
           IF (debug) THEN
261
              WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - original order"
262
              DO I=1,Nat
263
                 Iat=Order(I)
264
                 WRITE(*,'(1X,I5,3(1X,F20.6))') I,GradCart(3*Iat-2:3*Iat)
265
              END DO
266
              WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - internal order"
267
              DO I=1,Nat
268
                 Iat=Order(I)
269
                 WRITE(*,'(3(1X,I5),3(1X,F20.6))') I,Iat,OrderInv(I),GradCart(3*I-2:3*I)
270
              END DO
271
           END IF
272

    
273
           ! VASP gives the Forces in eV/Angstrom, so we convert it into the
274
           ! gradient in ua/Angstrom
275
           GradCart=-1.0d0*ev2Au*GradCart
276

    
277
           ! We convert it into the right coordinate system
278
           ! This should be put into a subroutine (because it is also in egradient)
279
           SELECT CASE (COORD)
280
           CASE ("ZMAT")
281
              ALLOCATE(GradTmp(3*Nat))
282
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
283
              CALL CalcBmat_int(nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome)
284
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
285

    
286
              Grad(IGeom,1)=GradTmp(4)
287
              Grad(IGeom,2)=GradTmp(7)
288
              ! We might have troubles whith angles that are not in the [0:pi] range because,
289
              ! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian...
290
              ! so that the derivative is not good, and a multiplicative factor should be added...
291
              !
292
              ! This in fact should be taken care of in calc Bmat ...
293
              !
294
              IF ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN
295
                 Grad(IGeom,3)=-1.0d0*GradTmp(8)
296
              ELSE
297
                 Grad(IGeom,3)=GradTmp(8)
298
              END IF
299
              Idx=4
300
              DO I=4,Nat
301
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
302
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
303
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
304
                 Idx=Idx+3
305
              END DO
306
              DEALLOCATE(GradTmp)
307
              ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
308
              ! but we have to convert it into internal coordinates if COORD=Baker      
309
           CASE ('BAKER')
310
              ALLOCATE(GradTmp(3*Nat))
311
              GradTmp=0.d0
312
              ! Below is to be corrected.
313
              ! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN'
314
              !DO I=1, 3*Nat
315
              ! here GradTmp and Grad are gradients in Baker coordinates
316
              ! Not implemented  IF ((PROG=="VASP").AND.(RunMode=="PARA"))
317
              ! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I)
318
              !END DO
319
              Grad(IGeom,:) = GradTmp
320
              DEALLOCATE(GradTmp)
321
           CASE ('MIXED')
322
              ALLOCATE(GradTmp(3*Nat))
323
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
324
              CALL CalcBMat_mixed (nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome)
325
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
326
              DO I=1,Nat
327
                 WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I-2:3*I)
328
              END DO
329
              SELECT CASE (NCART)
330
              CASE (1)
331
                 Grad(IGeom,1:3)=GradTmp(1:3)
332
                 Grad(IGeom,4)=GradTmp(4)
333
                 Grad(IGeom,5)=GradTmp(7)
334
                 Grad(IGeom,6)=GradTmp(8)
335
                 Idx=7
336
                 IBeg=4
337
              CASE(2)
338
                 Grad(IGeom,1:3)=GradTmp(1:3)
339
                 Grad(IGeom,4:6)=GradTmp(4:6)
340
                 Grad(IGeom,7)=GradTmp(7)
341
                 Grad(IGeom,8)=GradTmp(8)
342
                 Idx=9
343
                 IBeg=4
344
              CASE DEFAULT
345
                 Idx=1
346
                 IBeg=1
347
              END SELECT
348
              DO I=IBeg,Nat
349
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
350
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
351
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
352
                 Idx=Idx+3
353
              END DO
354
              DEALLOCATE(GradTmp)
355
           CASE ("CART","HYBRID")
356
              Grad(IGeom,:)=GradCart
357
           CASE DEFAULT
358
              WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP"
359
              STOP
360
           END SELECT
361

    
362
           IF (debug) THEN
363
              Line='DBG Path, GradTmp'
364
              Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat)
365
           END IF
366
        END DO ! matches  DO IGeom=IGeom0,IGeomF
367

    
368
     CASE ("GAUSSIAN")
369
        ! We create the output files.
370
        IGeom0=2
371
        IGeomF=NGeomF-1
372
        PathCart=0.d0
373

    
374
        If (OptReac) IGeom0=1
375
        If (OptProd) IGeomF=NGeomF
376

    
377
        OPEN(IOTMP2,File="ListJob",STATUS="REPLACE")
378

    
379
        DO Igeom=IGeom0,IGeomF
380
           ! We first convert the geometry into a Cartesian one
381
           SELECT CASE(Coord)
382
           CASE ('ZMAT')
383
              GeomTmp=IntCoordF(IGeom,:)
384
              Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
385
           CASE ('BAKER')
386
              GeomTmp=IntCoordF(IGeom,:)
387
              !
388
              IF (IOpt .EQ. 0) THEN
389
                 ! EgradPath(...) is called only after the call of PathCreate(...)
390
                 ! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always 
391
                 ! true if IOpt==0), THEN ExtraPol_baker(...) is called.
392
                 ! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting
393
                 ! the interpolated internal coordinates into cartesian coordinates.
394
                 ! Thus here we don't need to reconvert the internal coordinates again
395
                 ! to the cartesian coordinates. 
396
                 DO I=1,Nat
397
                    PathCart(3*I-2,IGeom) = XyzGeomF(IGeom,1,I)
398
                    PathCart(3*I-1,IGeom) = XyzGeomF(IGeom,2,I)
399
                    PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I)
400
                 END DO
401
              ELSE
402

    
403
                 ! XyzGeomF(NGeomF,3,Nat).
404
                 x_k(:) = XyzGeomF(IGeom,1,:)
405
                 y_k(:) = XyzGeomF(IGeom,2,:)
406
                 z_k(:) = XyzGeomF(IGeom,3,:)
407

    
408
                 !Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), &
409
                 !                               x_k,y_k,z_k,x,y,z)
410
                 ! PathCart(3*Nat,NGeomF)
411
                 DO I=1,Nat
412
                    PathCart(3*I-2,IGeom) = x(I)
413
                    PathCart(3*I-1,IGeom) = y(I)
414
                    PathCart(3*I,IGeom) = z(I)
415
                 END DO
416

    
417
                 XyzGeomF(IGeom,1,:)=x(:)
418
                 XyzGeomF(IGeom,2,:)=y(:)
419
                 XyzGeomF(IGeom,3,:)=z(:)
420
              END IF
421
           CASE ('MIXED')
422
              GeomTmp=IntCoordF(IGeom,:)
423
              Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
424
           CASE ('CART','HYBRID')
425
!!! CAUTION : PFL 29.AUG.2008 ->
426
              ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
427
              ! This should be made  consistent !!!!!!!!!!!!!!!
428
              GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
429
!!! <- CAUTION : PFL 29.AUG.2008 
430
              PathCart(:,IGeom)=GeomTmp
431
           CASE DEFAULT
432
              WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L75.STOP"
433
              STOP
434
           END SELECT
435

    
436
           WRITE(Line,'(I5)') IGeom
437
           FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".com"
438
           WRITE(IOTMP2,'(1X,A)') TRIM(FileInp)
439

    
440
           OPEN(IOTMP,File=TRIM(FileInp))
441
           Current => Gauss_root
442
           DO WHILE (ASSOCIATED(Current%next))
443
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
444
              WRITE(*,'(1X,A)') Trim(current%line)
445
              Current => current%next
446
           END DO
447

    
448
           WRITE(IOTMP,*) 
449

    
450
           Current => Gauss_comment
451
           DO WHILE (ASSOCIATED(Current%next))
452
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
453
              Current => current%next
454
           END DO
455
           !  WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)"),Igeom,NgeomF,Iopt
456

    
457
           WRITE(IOTMP,*) 
458
           WRITE(IOTMP,*) Trim(Gauss_charge)
459

    
460
           DO I=1,Nat
461
              If (renum) THEN
462
                 Iat=Order(I)
463
                 WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)) &
464
                      ,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom)  &
465
                      ,TRIM(Gauss_Paste(I))
466
              ELSE
467
                 Iat=OrderInv(I)
468
                 WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)) &
469
                 ,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom)  &
470
                      , TRIM(Gauss_Paste(Iat))
471
              END IF
472
           END DO
473

    
474
           WRITE(IOTMP,*) 
475

    
476
           Current => Gauss_End
477
           DO WHILE (ASSOCIATED(Current%next))
478
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
479
              Current => current%next
480
           END DO
481

    
482
           WRITE(IOTMP,*) 
483
           CLOSE(IOTMP)
484

    
485
        END DO
486

    
487
        Close(IOTMP2)
488

    
489
! We launch the parallel calculations
490

    
491
        Line=Trim(ProgExe) 
492
        IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(Line)
493
        call system(Line)
494
        IF (DEBUG) WRITE(*,*)'Back from RunCommand:'
495

    
496
! We read energies and gradients from the log files
497

    
498
        DO IGeom=IGeom0,IGeomF
499
           WRITE(Line,'(I5)') IGeom
500
           FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".log"
501
           OPEN(IOTMP,File=TRIM(FileInp))
502

    
503
  ! We first search for the forces
504
  READ(IOTMP,'(A)') LINE
505
  DO WHILE (INDEX(LINE,'Forces (Hartrees/Bohr)')==0)
506
     ! here, we have a problem, because there are so many ways to write E for Gaussian :-)
507
     ! but we do not really need E except if we want to weight the different points
508
     ! on the path... that will be done later :-)
509
     ! And I should use ConvGaussXmol ...
510
! PFL 3rd Jan 2008 ->
511
! I have finally upgraded the energy reading phase so that it should be able to read 
512
! many formats corresponding to HF, DFT, MP2, AM1, PM3, ONIOM with or w/o PCM
513
!
514
     
515
! For MP2, energy is after the last =
516
     IF ((Line(2:9)=="E2 =    ")) THEN
517
      Itmp1=Index(LINE,"=",BACK=.TRUE.)+1
518
    READ(LINE(Itmp1:),*) e
519
      END IF
520
! For other methods, it is after the first =
521
     IF ((LINE(2:9)=='Energy= ').OR.(Line(2:9)=="SCF Done").OR.(Line(2:9)=="ONIOM: e") &
522
   .OR.(Line(2:9)==" with al") &
523
         .OR.(Line(2:31)=="Counterpoise: corrected energy")) THEN
524
      Itmp1=Index(LINE,"=")+1
525
    READ(LINE(Itmp1:),*) e
526
  END IF
527
! <- PFL 3 Jan 2008
528
     READ(IOTMP,'(A)') LINE
529
  END DO
530
  READ(IOTMP,'(A)') LINE
531
  READ(IOTMP,'(A)') LINE
532
  DO I=1,Nat
533
     Iat=I
534
     IF (renum) Iat=Order(I)
535
     READ(IOTMP,*)  Itmp1,ITmp2, GradCart(3*Iat-2:3*Iat)
536
  END DO
537

    
538
! Gaussian gives the Forces in ua/a0, so we convert it into the 
539
! gradient in ua/Angstrom 
540
  gradCart=-1.0d0*Unit*gradCart
541

    
542
  CLOSE(IOTMP)
543
  Energies(IGeom)=E
544
  
545
           ! We convert it into the right coordinate system
546
           ! This should be put into a subroutine (because it is also in egradient)
547
           SELECT CASE (COORD)
548
           CASE ("ZMAT")
549
              ALLOCATE(GradTmp(3*Nat))
550
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
551
              CALL CalcBMat_int (nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome)
552
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
553

    
554
              Grad(IGeom,1)=GradTmp(4)
555
              Grad(IGeom,2)=GradTmp(7)
556
              ! We might have troubles whith angles that are not in the [0:pi] range because,
557
              ! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian...
558
              ! so that the derivative is not good, and a multiplicative factor should be added...
559
              !
560
              ! This in fact should be taken care of in B-mat calculation...
561
              !
562
              IF ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN
563
                 Grad(IGeom,3)=-1.0d0*GradTmp(8)
564
              ELSE
565
                 Grad(IGeom,3)=GradTmp(8)
566
              END IF
567
              Idx=4
568
              DO I=4,Nat
569
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
570
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
571
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
572
                 Idx=Idx+3
573
              END DO
574
              DEALLOCATE(GradTmp)
575
              ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
576
              ! but we have to convert it into internal coordinates if COORD=Baker      
577
           CASE ('BAKER')
578
              ALLOCATE(GradTmp(3*Nat))
579
              GradTmp=0.d0
580
              ! Below is to be corrected.
581
              ! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN'
582
              !DO I=1, 3*Nat
583
              ! here GradTmp and Grad are gradients in Baker coordinates
584
              ! Not implemented  IF ((PROG=="VASP").AND.(RunMode=="PARA"))
585
              ! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I)
586
              !END DO
587
              Grad(IGeom,:) = GradTmp
588
              DEALLOCATE(GradTmp)
589
           CASE ('MIXED')
590
              ALLOCATE(GradTmp(3*Nat))
591
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
592
              CALL CalcBMat_mixed(nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome)
593
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
594
              DO I=1,Nat
595
                 WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I-2:3*I)
596
              END DO
597
              SELECT CASE (NCART)
598
              CASE (1)
599
                 Grad(IGeom,1:3)=GradTmp(1:3)
600
                 Grad(IGeom,4)=GradTmp(4)
601
                 Grad(IGeom,5)=GradTmp(7)
602
                 Grad(IGeom,6)=GradTmp(8)
603
                 Idx=7
604
                 IBeg=4
605
              CASE(2)
606
                 Grad(IGeom,1:3)=GradTmp(1:3)
607
                 Grad(IGeom,4:6)=GradTmp(4:6)
608
                 Grad(IGeom,7)=GradTmp(7)
609
                 Grad(IGeom,8)=GradTmp(8)
610
                 Idx=9
611
                 IBeg=4
612
              CASE DEFAULT
613
                 Idx=1
614
                 IBeg=1
615
              END SELECT
616
              DO I=IBeg,Nat
617
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
618
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
619
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
620
                 Idx=Idx+3
621
              END DO
622
              DEALLOCATE(GradTmp)
623
           CASE ("CART","HYBRID")
624
              Grad(IGeom,:)=GradCart
625
           CASE DEFAULT
626
              WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP"
627
              STOP
628
           END SELECT
629

    
630
           IF (debug) THEN
631
              Line='DBG Path, GradTmp'
632
              Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat)
633
           END IF
634

    
635
        END DO
636

    
637
     CASE DEFAULT
638
        WRITE(*,*) "Prog=",TRIM(Prog)," not yet supported for Para mode."
639
        STOP
640
     END SELECT
641

    
642
     DEALLOCATE(GradCart,PathCart)
643

    
644
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
645
!
646
! Serial executions
647
!
648
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
649
  ELSE ! matches IF (RunMode=="PARA") THEN
650
 ! We will launch all calculations sequentially
651
     ALLOCATE(GradTmp(NCoord))
652
     ! We have the new path, we calculate its energy and gradient
653
     IGeom0=2
654
     IGeomF=NGeomF-1
655
     IF (OptReac.OR.CalcEReac) IGeom0=1
656
     IF (OptProd.OR.CalcEProd) IGeomF=NGeomF
657
     CalcEReac=.FALSE.
658
     CalcEProd=.FALSE.
659
     ALLOCATE(GeomCart(Nat,3))
660

    
661
     DO IGeom=IGeom0,IGeomF
662
        WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt
663
        if (debug) WRITE(*,*) Line
664
        SELECT CASE(Coord)
665
        CASE ('ZMAT','MIXED')
666
           GeomTmp=IntCoordF(IGeom,:)
667
        CASE ('BAKER')
668
           GeomTmp=IntCoordF(IGeom,:)
669
        CASE ('CART','HYBRID')
670
!!! CAUTION : PFL 29.AUG.2008 ->
671
           ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
672
           ! This should be made  consistent !!!!!!!!!!!!!!!
673
           GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
674
           !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
675
!!! <- CAUTION : PFL 29.AUG.2008
676
        CASE DEFAULT
677
           WRITE(*,*) "Coord=",TRIM(Coord)," not recognized in EgradPath. L299"
678
           STOP
679
        END SELECT
680

    
681
        !if (debug) WRITE(*,*) "Calling PrintGeom"
682
        !Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
683

    
684
        !           IF (debug) THEN 
685
        !              WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD), ' Igeom=',IGeom
686
        !              WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
687
        !              WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
688
        !              DO Iat=1,Nat
689
        !                 WRITE(*,*) XyzGeomF(IGeom,1:3,Iat)
690
        !              END DO
691
        !           END IF
692

    
693
        IF (debug) THEN 
694
           WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD),' IGeom=',IGeom,' Iopt=',Iopt
695
           WRITE(*,*) GeomTmp(1:min(NCoord,12))
696
        END IF
697

    
698
        GeomCart=0.d0
699
!        IF(IOpt .EQ. 1) Then
700
!           Print *, 'GeomTmp='
701
!           WRITE(*,'(12(1X,F6.3))') GeomTmp
702
!           Print *, 'GeomCart='
703
!           WRITE(*,'(12(1X,F6.3))') GeomCart
704
!      END IF
705
    
706
       ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
707
       ! GeomCart has INTENT(OUT)
708
       Call EGrad(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart) !GeomTmp=IntCoordF(IGeom,:)
709

    
710
    ! Egrad calls ConvertBakerInternal_cart to convert GeomTmp=IntCoordF(IGeom,:) into 
711
    ! GeomCart (Cartesian Coordinates) so that energy and gradients can be calculated 
712
    ! for each geometry.
713
        XyzGeomF(IGeom,1,:)=GeomCart(:,1)
714
        XyzGeomF(IGeom,2,:)=GeomCart(:,2)
715
        XyzGeomF(IGeom,3,:)=GeomCart(:,3)
716
      
717
        Energies(IGeom)=E
718

    
719
        if (debug) THEN
720
           Line='DBG Path, GradTmp'
721
           Call PrintGeom(Line,Nat,NCoord,GRadTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
722
        END IF
723
        Grad(IGeom,:)=GradTmp
724
     END DO ! matches DO IGeom=IGeom0,IGeomF
725
     DEALLOCATE(GradTmp,GeomCart)
726
  END IF ! matches IF (RunMode=="PARA") THEN AFTER ELSE
727

    
728
  DEALLOCATE(GeomTmp)
729
  DEALLOCATE(x,y,z)
730
  DEALLOCATE(x_k,y_k,z_k)
731
  if (debug) Call header('Exiting EgradPath')
732
  RETURN
733
999 WRITE(*,*) "EgradPath : We should not be here !"
734
  STOP
735

    
736
END SUBROUTINE EGRADPATH