Statistiques
| Révision :

root / src / EgradPath.f90 @ 8

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

1
SUBROUTINE EGradPath(Iopt,Flag_Opt_Geom)
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, IntCoordI, GeomOld_all, 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
  LOGICAL, INTENT(IN) :: Flag_Opt_Geom
17

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

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

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

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

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

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

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

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

    
57
  use Io_module
58

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

    
79
 END subroutine Egrad
80

    
81
  END INTERFACE
82

    
83
  debug=valid('EGradPath')
84

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

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

    
91

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

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

    
96

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    
259
           Close(IoTmp)
260

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    
449
           WRITE(IOTMP,*) 
450

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

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

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

    
475
           WRITE(IOTMP,*) 
476

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

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

    
486
        END DO
487

    
488
        Close(IOTMP2)
489

    
490
! We launch the parallel calculations
491

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

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

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

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

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

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

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

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

    
636
        END DO
637

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

    
643
     DEALLOCATE(GradCart,PathCart)
644

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

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

    
677
        !if (debug) WRITE(*,*) "Calling PrintGeom"
678
        !Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
679

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

    
689
        IF (debug) THEN 
690
           WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD),' IGeom=',IGeom,' Iopt=',Iopt
691
           WRITE(*,*) GeomTmp(1:min(NCoord,12))
692
        END IF
693

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

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

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

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

    
732
END SUBROUTINE EGRADPATH