Statistics
| Revision:

root / src / EgradPath.f90 @ 8

History | View | Annotate | Download (27.7 kB)

1 2 pfleura2
SUBROUTINE EGradPath(Iopt)
2 1 pfleura2
3 1 pfleura2
  ! This subroutine calculates the energy+gradient for all points of
4 1 pfleura2
  ! the path
5 1 pfleura2
6 1 pfleura2
  use Path_module, only : ProgExe, NCoord, PROG, Nat, NGeomF, Coord, IntCoordF, &
7 1 pfleura2
       IndZmat, XyzGeomF, Atome, Order, OrderInv, Latr, FrozAtoms, &
8 1 pfleura2
       Grad, Energies, renum, dzdc, massat, Pi, NCART, OptReac, &
9 8 pfleura2
       OptProd,  AtName, Unit, CalcEReac, CalcEProd
10 1 pfleura2
  ! GeomOld_all(NGeomF,NCoord), BTransInv (3*Nat,NCoord): used for Baker case
11 1 pfleura2
  use Io_module
12 1 pfleura2
13 1 pfleura2
  IMPLICIT NONE
14 1 pfleura2
15 1 pfleura2
  INTEGER(KINT), INTENT(IN) :: Iopt
16 1 pfleura2
17 1 pfleura2
  INTEGER(KINT) :: IGeom,IGeom0,IGeomF
18 1 pfleura2
  CHARACTER(LCHARS) :: Line,FileInp,RunCommand
19 1 pfleura2
  CHARACTER(SCHARS) :: Line2
20 1 pfleura2
21 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:) ! NCoord
22 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: GradCart(:) ! 3*Nat
23 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: PathCart(:,:) ! (3*Nat,NGeomF)
24 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
25 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: x_k(:), y_k(:), z_k(:)
26 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3)
27 1 pfleura2
  REAL(KREAL) :: E
28 1 pfleura2
  LOGICAL :: Debug
29 1 pfleura2
  LOGICAL, SAVE :: First=.TRUE.
30 1 pfleura2
31 1 pfleura2
  INTEGER(KINT) :: I, Iat, IBeg, Idx
32 1 pfleura2
  REAL(KREAL) :: RTmp1, RTmp2, RTmp3
33 1 pfleura2
  INTEGER(KINT) :: Itmp1, ITmp2
34 1 pfleura2
35 1 pfleura2
  INTERFACE
36 1 pfleura2
     function valid(string) result (isValid)
37 1 pfleura2
       CHARACTER(*), intent(in) :: string
38 1 pfleura2
       logical                  :: isValid
39 1 pfleura2
     END function VALID
40 1 pfleura2
41 1 pfleura2
     subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old)
42 1 pfleura2
43 1 pfleura2
  ! This routines calculates the energy E and the gradient Grad of
44 1 pfleura2
  ! a molecule with Geometry Geom (may be in internal coordinates),
45 1 pfleura2
  ! using for now, either Gaussian or Ext, more general later.
46 1 pfleura2
47 1 pfleura2
  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
48 1 pfleura2
       prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &
49 1 pfleura2
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
50 1 pfleura2
       , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
51 1 pfleura2
      , Order,OrderInv, XPrimitiveF
52 1 pfleura2
53 1 pfleura2
  ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
54 1 pfleura2
  ! allocated in Path.f90
55 1 pfleura2
56 1 pfleura2
  use Io_module
57 1 pfleura2
58 1 pfleura2
  ! Energy (calculated if F300K=.F., else estimated)
59 1 pfleura2
  REAL(KREAL), INTENT (OUT) :: E
60 1 pfleura2
  ! NCoord: Number of the degrees of freedom
61 1 pfleura2
  ! IGeom: index of the geometry.
62 1 pfleura2
  INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
63 1 pfleura2
  ! Geometry at which gradient is calculated (cf Factual also):
64 1 pfleura2
  REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
65 1 pfleura2
  ! Gradient calculated at Geom geometry:
66 1 pfleura2
  REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
67 1 pfleura2
  ! Cartesian geometry corresponding to (Internal Geometry) Geom:
68 1 pfleura2
  REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)
69 1 pfleura2
!!! Optional, just for geometry optimization with Baker coordinates
70 1 pfleura2
  REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
71 1 pfleura2
  REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
72 1 pfleura2
! FOptGeom is a flag indicating if we are doing a geom optimization
73 1 pfleura2
! it can be omitted so that we use a local flag: Flag_Opt_Geom
74 1 pfleura2
  LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
75 1 pfleura2
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
76 1 pfleura2
  LOGICAL  :: Flag_Opt_Geom
77 1 pfleura2
78 1 pfleura2
 END subroutine Egrad
79 1 pfleura2
80 1 pfleura2
  END INTERFACE
81 1 pfleura2
82 1 pfleura2
  debug=valid('EGradPath')
83 1 pfleura2
84 1 pfleura2
  if (debug) Call header("Entering EgradPath")
85 1 pfleura2
86 1 pfleura2
  ALLOCATE(GeomTmp(NCoord))
87 1 pfleura2
  ALLOCATE(x(Nat),y(Nat),z(Nat))
88 1 pfleura2
  ALLOCATE(x_k(Nat),y_k(Nat),z_k(Nat))
89 1 pfleura2
90 1 pfleura2
91 1 pfleura2
  IF (RunMode=="PARA") THEN ! matches at L315
92 1 pfleura2
93 1 pfleura2
     ALLOCATE(GradCart(3*Nat),PathCart(3*Nat,NGeomF))
94 1 pfleura2
95 1 pfleura2
96 1 pfleura2
     SELECT CASE(Prog)
97 1 pfleura2
     CASE ("VASP")
98 1 pfleura2
        ! We will use the NEB routine of VASP to get all forces at once on Para8 queue
99 1 pfleura2
100 1 pfleura2
        ! First, we create all the POSCAR into the 0x directories
101 1 pfleura2
102 1 pfleura2
        ! With RunMode="Para" one cannot optimize reactants or products (due to Vasp constraints)
103 1 pfleura2
        ! so this has been done beforehand :-)
104 1 pfleura2
        IGeom0=2
105 1 pfleura2
        IGeomF=NGeomF-1
106 1 pfleura2
        PathCart=0.d0
107 1 pfleura2
108 1 pfleura2
        IF (First) THEN
109 1 pfleura2
           ! For the first iteration, we do write a POSCAR into 00 and 0Final
110 1 pfleura2
           ! Vasp needs it !
111 1 pfleura2
           First=.FALSE.
112 1 pfleura2
           IGeom0=1
113 1 pfleura2
           IGeomF=NGeomF
114 1 pfleura2
        END IF
115 1 pfleura2
116 1 pfleura2
        Line2="00"
117 1 pfleura2
        DO IGeom=IGeom0,IGeomF
118 1 pfleura2
119 1 pfleura2
           ! NOTE THAT HERE IGeom0 (IGeomF) IS NOT 1 (NGeomF) FOR NON-FIRST OPTIMIZATION.
120 1 pfleura2
           ! THIS WILL HAVE EFFECT ON XyzGeomF UPDATE IN BAKER CASE.
121 1 pfleura2
           WRITE(Line,'(I5)') IGeom-1
122 1 pfleura2
           Idx=2-Len_TRIM(ADJUSTL(Line))
123 1 pfleura2
           FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/POSCAR"
124 1 pfleura2
           if (debug) WRITE(*,*) "Creating ",TRIM(FileInp)
125 1 pfleura2
           OPEN(IOTMP,File=TRIM(FileInp))
126 1 pfleura2
127 1 pfleura2
           WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt
128 1 pfleura2
           if (debug) WRITE(*,*) Line
129 1 pfleura2
           SELECT CASE(Coord)
130 1 pfleura2
           CASE ('ZMAT')
131 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
132 1 pfleura2
              Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
133 1 pfleura2
           CASE ('BAKER')
134 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
135 1 pfleura2
              !
136 1 pfleura2
              IF (IOpt .EQ. 0) THEN
137 1 pfleura2
                 ! EgradPath(...) is called only after the call of PathCreate(...)
138 1 pfleura2
                 ! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always
139 1 pfleura2
                 ! true if IOpt==0), THEN ExtraPol_baker(...) is called.
140 1 pfleura2
                 ! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting
141 1 pfleura2
                 ! the interpolated internal coordinates into cartesian coordinates.
142 1 pfleura2
                 ! Thus here we don't need to reconvert the internal coordinates again
143 1 pfleura2
                 ! to the cartesian coordinates.
144 1 pfleura2
                 DO I=1,Nat
145 1 pfleura2
                    PathCart(3*I-2,IGeom) = XyzGeomF(IGeom,1,I)
146 1 pfleura2
                    PathCart(3*I-1,IGeom) = XyzGeomF(IGeom,2,I)
147 1 pfleura2
                    PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I)
148 1 pfleura2
                 END DO
149 1 pfleura2
              ELSE
150 1 pfleura2
                 ! XyzGeomF(NGeomF,3,Nat).
151 1 pfleura2
                 x_k(:) = XyzGeomF(IGeom,1,:)
152 1 pfleura2
                 y_k(:) = XyzGeomF(IGeom,2,:)
153 1 pfleura2
                 z_k(:) = XyzGeomF(IGeom,3,:)
154 1 pfleura2
155 1 pfleura2
                 !Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), &
156 1 pfleura2
                 !                               x_k,y_k,z_k,x,y,z)
157 1 pfleura2
                 ! PathCart(3*Nat,NGeomF)
158 1 pfleura2
                 DO I=1,Nat
159 1 pfleura2
                    PathCart(3*I-2,IGeom) = x(I)
160 1 pfleura2
                    PathCart(3*I-1,IGeom) = y(I)
161 1 pfleura2
                    PathCart(3*I,IGeom) = z(I)
162 1 pfleura2
                 END DO
163 1 pfleura2
164 1 pfleura2
                 XyzGeomF(IGeom,1,:)=x(:)
165 1 pfleura2
                 XyzGeomF(IGeom,2,:)=y(:)
166 1 pfleura2
                 XyzGeomF(IGeom,3,:)=z(:)
167 1 pfleura2
              END IF
168 1 pfleura2
           CASE ('MIXED')
169 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
170 1 pfleura2
              Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
171 1 pfleura2
           CASE ('CART','HYBRID')
172 1 pfleura2
!!! CAUTION : PFL 29.AUG.2008 ->
173 1 pfleura2
              ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
174 1 pfleura2
              ! This should be made  consistent !!!!!!!!!!!!!!!
175 1 pfleura2
              GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
176 1 pfleura2
!!! <- CAUTION : PFL 29.AUG.2008
177 1 pfleura2
            PathCart(:,IGeom)=GeomTmp
178 1 pfleura2
          CASE DEFAULT
179 1 pfleura2
            WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L134. STOP"
180 1 pfleura2
            STOP
181 1 pfleura2
        END SELECT
182 1 pfleura2
183 1 pfleura2
           IF (debug) THEN
184 1 pfleura2
              WRITE(*,*) "Calling PrintGeomVasp for geom",Igeom
185 1 pfleura2
              Call PrintGeom(Line,Nat,NCoord,PathCart(1,Igeom),Coord,IoOut,Atome,Order,OrderInv,IndZmat)
186 1 pfleura2
           END IF
187 1 pfleura2
           Call PrintGeomVasp(Line,PathCart(1,IGeom),Latr,FrozAtoms,IoTmp)
188 1 pfleura2
189 1 pfleura2
           !           IF (debug) THEN
190 1 pfleura2
           !              WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD), ' Igeom=',IGeom
191 1 pfleura2
           !              WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
192 1 pfleura2
           !              WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
193 1 pfleura2
           !              DO Iat=1,Nat
194 1 pfleura2
           !                 WRITE(*,*) XyzGeomF(IGeom,1:3,Iat)
195 1 pfleura2
           !              END DO
196 1 pfleura2
           !           END IF
197 1 pfleura2
198 1 pfleura2
           IF (debug) THEN
199 1 pfleura2
              WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD),' Igeom=',IGeom,' Iopt=',Iopt
200 1 pfleura2
              WRITE(*,*) GeomTmp(1:min(NCoord,9))
201 1 pfleura2
           END IF
202 1 pfleura2
203 1 pfleura2
           CLOSE(IOTMP)
204 1 pfleura2
        END DO ! matches DO IGeom=IGeom0,IGeomF
205 1 pfleura2
206 1 pfleura2
!!!!!!!!!!!!!!!!!!!
207 1 pfleura2
!
208 1 pfleura2
        ! We calculate the energies and gradients
209 1 pfleura2
!
210 1 pfleura2
211 1 pfleura2
        RunCommand=Trim(ProgExe)
212 1 pfleura2
        if (debug) WRITE(*,*) "Launching:", TRIM(RunCommand)
213 1 pfleura2
        call system(RunCommand)
214 1 pfleura2
        if (debug) WRITE(*,*) "Back from:", TRIM(RunCommand)
215 1 pfleura2
216 1 pfleura2
!!!!!!!!!!!!!
217 1 pfleura2
!
218 1 pfleura2
        ! We read the gradients and energies
219 1 pfleura2
!
220 1 pfleura2
        IGeom0=2
221 1 pfleura2
        IGeomF=NGeomF-1
222 1 pfleura2
        Grad=0.d0
223 1 pfleura2
224 1 pfleura2
        Line2="00"
225 1 pfleura2
        DO IGeom=IGeom0,IGeomF
226 1 pfleura2
           WRITE(Line,'(I5)') IGeom-1
227 1 pfleura2
           Idx=2-Len_TRIM(ADJUSTL(Line))
228 1 pfleura2
           FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/OUTCAR"
229 1 pfleura2
           if (debug) WRITE(*,*) "Reading E and forces from ", TRIM(FileInp)
230 1 pfleura2
           OPEN(IOTMP,File=TRIM(FileInp))
231 1 pfleura2
232 1 pfleura2
           ! We first search for the forces
233 1 pfleura2
           READ(IOTMP,'(A)',END=999,ERR=999) LINE
234 1 pfleura2
           ! We search for the  last part of the OUTCAR file, after wavefunction is converged
235 1 pfleura2
           DO WHILE (INDEX(LINE,'EDIFF is reached')==0)
236 1 pfleura2
              READ(IOTMP,'(A)',END=999,ERR=999) LINE
237 1 pfleura2
           END DO
238 1 pfleura2
           READ(IOTMP,'(A)',END=999,ERR=999) LINE
239 1 pfleura2
           DO WHILE (INDEX(LINE,'TOTEN')==0)
240 1 pfleura2
              READ(IOTMP,'(A)',END=999,ERR=999) LINE
241 1 pfleura2
           END DO
242 1 pfleura2
           Line=Line(26:)
243 1 pfleura2
           READ(LINE,*) E
244 1 pfleura2
           if (debug) WRITE(*,'(1X,A,I5,A,F20.6)') "DBG Egrad_VASP E(",Igeom,")=",E
245 1 pfleura2
           Energies(Igeom)=E
246 1 pfleura2
247 1 pfleura2
           ! We search for the forces
248 1 pfleura2
           DO WHILE (INDEX(LINE,'TOTAL-FORCE')==0)
249 1 pfleura2
              READ(IOTMP,'(A)',END=999,ERR=999) LINE
250 1 pfleura2
           END DO
251 1 pfleura2
           READ(IOTMP,'(A)',END=999,ERR=999) LINE
252 1 pfleura2
           DO I=1,Nat
253 1 pfleura2
              Iat=I
254 1 pfleura2
              IF (renum) Iat=Order(I)
255 1 pfleura2
              READ(IOTMP,*,END=999,ERR=999)  Rtmp1,RTmp2,RTmp3,GradCart(3*Iat-2:3*Iat)
256 1 pfleura2
           END DO
257 1 pfleura2
258 1 pfleura2
           Close(IoTmp)
259 1 pfleura2
260 1 pfleura2
           IF (debug) THEN
261 1 pfleura2
              WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - original order"
262 1 pfleura2
              DO I=1,Nat
263 1 pfleura2
                 Iat=Order(I)
264 1 pfleura2
                 WRITE(*,'(1X,I5,3(1X,F20.6))') I,GradCart(3*Iat-2:3*Iat)
265 1 pfleura2
              END DO
266 1 pfleura2
              WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - internal order"
267 1 pfleura2
              DO I=1,Nat
268 1 pfleura2
                 Iat=Order(I)
269 1 pfleura2
                 WRITE(*,'(3(1X,I5),3(1X,F20.6))') I,Iat,OrderInv(I),GradCart(3*I-2:3*I)
270 1 pfleura2
              END DO
271 1 pfleura2
           END IF
272 1 pfleura2
273 1 pfleura2
           ! VASP gives the Forces in eV/Angstrom, so we convert it into the
274 1 pfleura2
           ! gradient in ua/Angstrom
275 1 pfleura2
           GradCart=-1.0d0*ev2Au*GradCart
276 1 pfleura2
277 1 pfleura2
           ! We convert it into the right coordinate system
278 1 pfleura2
           ! This should be put into a subroutine (because it is also in egradient)
279 1 pfleura2
           SELECT CASE (COORD)
280 1 pfleura2
           CASE ("ZMAT")
281 1 pfleura2
              ALLOCATE(GradTmp(3*Nat))
282 1 pfleura2
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
283 1 pfleura2
              CALL CalcBmat_int(nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome)
284 1 pfleura2
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
285 1 pfleura2
286 1 pfleura2
              Grad(IGeom,1)=GradTmp(4)
287 1 pfleura2
              Grad(IGeom,2)=GradTmp(7)
288 1 pfleura2
              ! We might have troubles whith angles that are not in the [0:pi] range because,
289 1 pfleura2
              ! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian...
290 1 pfleura2
              ! so that the derivative is not good, and a multiplicative factor should be added...
291 1 pfleura2
              !
292 1 pfleura2
              ! This in fact should be taken care of in calc Bmat ...
293 1 pfleura2
              !
294 1 pfleura2
              IF ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN
295 1 pfleura2
                 Grad(IGeom,3)=-1.0d0*GradTmp(8)
296 1 pfleura2
              ELSE
297 1 pfleura2
                 Grad(IGeom,3)=GradTmp(8)
298 1 pfleura2
              END IF
299 1 pfleura2
              Idx=4
300 1 pfleura2
              DO I=4,Nat
301 1 pfleura2
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
302 1 pfleura2
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
303 1 pfleura2
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
304 1 pfleura2
                 Idx=Idx+3
305 1 pfleura2
              END DO
306 1 pfleura2
              DEALLOCATE(GradTmp)
307 1 pfleura2
              ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
308 1 pfleura2
              ! but we have to convert it into internal coordinates if COORD=Baker
309 1 pfleura2
           CASE ('BAKER')
310 1 pfleura2
              ALLOCATE(GradTmp(3*Nat))
311 1 pfleura2
              GradTmp=0.d0
312 1 pfleura2
              ! Below is to be corrected.
313 1 pfleura2
              ! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN'
314 1 pfleura2
              !DO I=1, 3*Nat
315 1 pfleura2
              ! here GradTmp and Grad are gradients in Baker coordinates
316 1 pfleura2
              ! Not implemented  IF ((PROG=="VASP").AND.(RunMode=="PARA"))
317 1 pfleura2
              ! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I)
318 1 pfleura2
              !END DO
319 1 pfleura2
              Grad(IGeom,:) = GradTmp
320 1 pfleura2
              DEALLOCATE(GradTmp)
321 1 pfleura2
           CASE ('MIXED')
322 1 pfleura2
              ALLOCATE(GradTmp(3*Nat))
323 1 pfleura2
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
324 1 pfleura2
              CALL CalcBMat_mixed (nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome)
325 1 pfleura2
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
326 1 pfleura2
              DO I=1,Nat
327 1 pfleura2
                 WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I-2:3*I)
328 1 pfleura2
              END DO
329 1 pfleura2
              SELECT CASE (NCART)
330 1 pfleura2
              CASE (1)
331 1 pfleura2
                 Grad(IGeom,1:3)=GradTmp(1:3)
332 1 pfleura2
                 Grad(IGeom,4)=GradTmp(4)
333 1 pfleura2
                 Grad(IGeom,5)=GradTmp(7)
334 1 pfleura2
                 Grad(IGeom,6)=GradTmp(8)
335 1 pfleura2
                 Idx=7
336 1 pfleura2
                 IBeg=4
337 1 pfleura2
              CASE(2)
338 1 pfleura2
                 Grad(IGeom,1:3)=GradTmp(1:3)
339 1 pfleura2
                 Grad(IGeom,4:6)=GradTmp(4:6)
340 1 pfleura2
                 Grad(IGeom,7)=GradTmp(7)
341 1 pfleura2
                 Grad(IGeom,8)=GradTmp(8)
342 1 pfleura2
                 Idx=9
343 1 pfleura2
                 IBeg=4
344 1 pfleura2
              CASE DEFAULT
345 1 pfleura2
                 Idx=1
346 1 pfleura2
                 IBeg=1
347 1 pfleura2
              END SELECT
348 1 pfleura2
              DO I=IBeg,Nat
349 1 pfleura2
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
350 1 pfleura2
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
351 1 pfleura2
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
352 1 pfleura2
                 Idx=Idx+3
353 1 pfleura2
              END DO
354 1 pfleura2
              DEALLOCATE(GradTmp)
355 1 pfleura2
           CASE ("CART","HYBRID")
356 1 pfleura2
              Grad(IGeom,:)=GradCart
357 1 pfleura2
           CASE DEFAULT
358 1 pfleura2
              WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP"
359 1 pfleura2
              STOP
360 1 pfleura2
           END SELECT
361 1 pfleura2
362 1 pfleura2
           IF (debug) THEN
363 1 pfleura2
              Line='DBG Path, GradTmp'
364 1 pfleura2
              Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat)
365 1 pfleura2
           END IF
366 1 pfleura2
        END DO ! matches  DO IGeom=IGeom0,IGeomF
367 1 pfleura2
368 1 pfleura2
     CASE ("GAUSSIAN")
369 1 pfleura2
        ! We create the output files.
370 1 pfleura2
        IGeom0=2
371 1 pfleura2
        IGeomF=NGeomF-1
372 1 pfleura2
        PathCart=0.d0
373 1 pfleura2
374 1 pfleura2
        If (OptReac) IGeom0=1
375 1 pfleura2
        If (OptProd) IGeomF=NGeomF
376 1 pfleura2
377 1 pfleura2
        OPEN(IOTMP2,File="ListJob",STATUS="REPLACE")
378 1 pfleura2
379 1 pfleura2
        DO Igeom=IGeom0,IGeomF
380 1 pfleura2
           ! We first convert the geometry into a Cartesian one
381 1 pfleura2
           SELECT CASE(Coord)
382 1 pfleura2
           CASE ('ZMAT')
383 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
384 1 pfleura2
              Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
385 1 pfleura2
           CASE ('BAKER')
386 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
387 1 pfleura2
              !
388 1 pfleura2
              IF (IOpt .EQ. 0) THEN
389 1 pfleura2
                 ! EgradPath(...) is called only after the call of PathCreate(...)
390 1 pfleura2
                 ! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always
391 1 pfleura2
                 ! true if IOpt==0), THEN ExtraPol_baker(...) is called.
392 1 pfleura2
                 ! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting
393 1 pfleura2
                 ! the interpolated internal coordinates into cartesian coordinates.
394 1 pfleura2
                 ! Thus here we don't need to reconvert the internal coordinates again
395 1 pfleura2
                 ! to the cartesian coordinates.
396 1 pfleura2
                 DO I=1,Nat
397 1 pfleura2
                    PathCart(3*I-2,IGeom) = XyzGeomF(IGeom,1,I)
398 1 pfleura2
                    PathCart(3*I-1,IGeom) = XyzGeomF(IGeom,2,I)
399 1 pfleura2
                    PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I)
400 1 pfleura2
                 END DO
401 1 pfleura2
              ELSE
402 1 pfleura2
403 1 pfleura2
                 ! XyzGeomF(NGeomF,3,Nat).
404 1 pfleura2
                 x_k(:) = XyzGeomF(IGeom,1,:)
405 1 pfleura2
                 y_k(:) = XyzGeomF(IGeom,2,:)
406 1 pfleura2
                 z_k(:) = XyzGeomF(IGeom,3,:)
407 1 pfleura2
408 1 pfleura2
                 !Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), &
409 1 pfleura2
                 !                               x_k,y_k,z_k,x,y,z)
410 1 pfleura2
                 ! PathCart(3*Nat,NGeomF)
411 1 pfleura2
                 DO I=1,Nat
412 1 pfleura2
                    PathCart(3*I-2,IGeom) = x(I)
413 1 pfleura2
                    PathCart(3*I-1,IGeom) = y(I)
414 1 pfleura2
                    PathCart(3*I,IGeom) = z(I)
415 1 pfleura2
                 END DO
416 1 pfleura2
417 1 pfleura2
                 XyzGeomF(IGeom,1,:)=x(:)
418 1 pfleura2
                 XyzGeomF(IGeom,2,:)=y(:)
419 1 pfleura2
                 XyzGeomF(IGeom,3,:)=z(:)
420 1 pfleura2
              END IF
421 1 pfleura2
           CASE ('MIXED')
422 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
423 1 pfleura2
              Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
424 1 pfleura2
           CASE ('CART','HYBRID')
425 1 pfleura2
!!! CAUTION : PFL 29.AUG.2008 ->
426 1 pfleura2
              ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
427 1 pfleura2
              ! This should be made  consistent !!!!!!!!!!!!!!!
428 1 pfleura2
              GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
429 1 pfleura2
!!! <- CAUTION : PFL 29.AUG.2008
430 1 pfleura2
              PathCart(:,IGeom)=GeomTmp
431 1 pfleura2
           CASE DEFAULT
432 1 pfleura2
              WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L75.STOP"
433 1 pfleura2
              STOP
434 1 pfleura2
           END SELECT
435 1 pfleura2
436 1 pfleura2
           WRITE(Line,'(I5)') IGeom
437 1 pfleura2
           FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".com"
438 1 pfleura2
           WRITE(IOTMP2,'(1X,A)') TRIM(FileInp)
439 1 pfleura2
440 1 pfleura2
           OPEN(IOTMP,File=TRIM(FileInp))
441 1 pfleura2
           Current => Gauss_root
442 1 pfleura2
           DO WHILE (ASSOCIATED(Current%next))
443 1 pfleura2
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
444 1 pfleura2
              WRITE(*,'(1X,A)') Trim(current%line)
445 1 pfleura2
              Current => current%next
446 1 pfleura2
           END DO
447 1 pfleura2
448 1 pfleura2
           WRITE(IOTMP,*)
449 1 pfleura2
450 1 pfleura2
           Current => Gauss_comment
451 1 pfleura2
           DO WHILE (ASSOCIATED(Current%next))
452 1 pfleura2
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
453 1 pfleura2
              Current => current%next
454 1 pfleura2
           END DO
455 1 pfleura2
           !  WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)"),Igeom,NgeomF,Iopt
456 1 pfleura2
457 1 pfleura2
           WRITE(IOTMP,*)
458 1 pfleura2
           WRITE(IOTMP,*) Trim(Gauss_charge)
459 1 pfleura2
460 1 pfleura2
           DO I=1,Nat
461 1 pfleura2
              If (renum) THEN
462 1 pfleura2
                 Iat=Order(I)
463 1 pfleura2
                 WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)) &
464 1 pfleura2
                      ,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom)  &
465 1 pfleura2
                      ,TRIM(Gauss_Paste(I))
466 1 pfleura2
              ELSE
467 1 pfleura2
                 Iat=OrderInv(I)
468 1 pfleura2
                 WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)) &
469 1 pfleura2
                 ,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom)  &
470 1 pfleura2
                      , TRIM(Gauss_Paste(Iat))
471 1 pfleura2
              END IF
472 1 pfleura2
           END DO
473 1 pfleura2
474 1 pfleura2
           WRITE(IOTMP,*)
475 1 pfleura2
476 1 pfleura2
           Current => Gauss_End
477 1 pfleura2
           DO WHILE (ASSOCIATED(Current%next))
478 1 pfleura2
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
479 1 pfleura2
              Current => current%next
480 1 pfleura2
           END DO
481 1 pfleura2
482 1 pfleura2
           WRITE(IOTMP,*)
483 1 pfleura2
           CLOSE(IOTMP)
484 1 pfleura2
485 1 pfleura2
        END DO
486 1 pfleura2
487 1 pfleura2
        Close(IOTMP2)
488 1 pfleura2
489 1 pfleura2
! We launch the parallel calculations
490 1 pfleura2
491 1 pfleura2
        Line=Trim(ProgExe)
492 1 pfleura2
        IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(Line)
493 1 pfleura2
        call system(Line)
494 1 pfleura2
        IF (DEBUG) WRITE(*,*)'Back from RunCommand:'
495 1 pfleura2
496 1 pfleura2
! We read energies and gradients from the log files
497 1 pfleura2
498 1 pfleura2
        DO IGeom=IGeom0,IGeomF
499 1 pfleura2
           WRITE(Line,'(I5)') IGeom
500 1 pfleura2
           FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".log"
501 1 pfleura2
           OPEN(IOTMP,File=TRIM(FileInp))
502 1 pfleura2
503 1 pfleura2
  ! We first search for the forces
504 1 pfleura2
  READ(IOTMP,'(A)') LINE
505 1 pfleura2
  DO WHILE (INDEX(LINE,'Forces (Hartrees/Bohr)')==0)
506 1 pfleura2
     ! here, we have a problem, because there are so many ways to write E for Gaussian :-)
507 1 pfleura2
     ! but we do not really need E except if we want to weight the different points
508 1 pfleura2
     ! on the path... that will be done later :-)
509 1 pfleura2
     ! And I should use ConvGaussXmol ...
510 1 pfleura2
! PFL 3rd Jan 2008 ->
511 1 pfleura2
! I have finally upgraded the energy reading phase so that it should be able to read
512 1 pfleura2
! many formats corresponding to HF, DFT, MP2, AM1, PM3, ONIOM with or w/o PCM
513 1 pfleura2
!
514 1 pfleura2
515 1 pfleura2
! For MP2, energy is after the last =
516 1 pfleura2
     IF ((Line(2:9)=="E2 =    ")) THEN
517 1 pfleura2
      Itmp1=Index(LINE,"=",BACK=.TRUE.)+1
518 1 pfleura2
    READ(LINE(Itmp1:),*) e
519 1 pfleura2
      END IF
520 1 pfleura2
! For other methods, it is after the first =
521 1 pfleura2
     IF ((LINE(2:9)=='Energy= ').OR.(Line(2:9)=="SCF Done").OR.(Line(2:9)=="ONIOM: e") &
522 1 pfleura2
   .OR.(Line(2:9)==" with al") &
523 1 pfleura2
         .OR.(Line(2:31)=="Counterpoise: corrected energy")) THEN
524 1 pfleura2
      Itmp1=Index(LINE,"=")+1
525 1 pfleura2
    READ(LINE(Itmp1:),*) e
526 1 pfleura2
  END IF
527 1 pfleura2
! <- PFL 3 Jan 2008
528 1 pfleura2
     READ(IOTMP,'(A)') LINE
529 1 pfleura2
  END DO
530 1 pfleura2
  READ(IOTMP,'(A)') LINE
531 1 pfleura2
  READ(IOTMP,'(A)') LINE
532 1 pfleura2
  DO I=1,Nat
533 1 pfleura2
     Iat=I
534 1 pfleura2
     IF (renum) Iat=Order(I)
535 1 pfleura2
     READ(IOTMP,*)  Itmp1,ITmp2, GradCart(3*Iat-2:3*Iat)
536 1 pfleura2
  END DO
537 1 pfleura2
538 1 pfleura2
! Gaussian gives the Forces in ua/a0, so we convert it into the
539 1 pfleura2
! gradient in ua/Angstrom
540 1 pfleura2
  gradCart=-1.0d0*Unit*gradCart
541 1 pfleura2
542 1 pfleura2
  CLOSE(IOTMP)
543 1 pfleura2
  Energies(IGeom)=E
544 1 pfleura2
545 1 pfleura2
           ! We convert it into the right coordinate system
546 1 pfleura2
           ! This should be put into a subroutine (because it is also in egradient)
547 1 pfleura2
           SELECT CASE (COORD)
548 1 pfleura2
           CASE ("ZMAT")
549 1 pfleura2
              ALLOCATE(GradTmp(3*Nat))
550 1 pfleura2
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
551 1 pfleura2
              CALL CalcBMat_int (nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome)
552 1 pfleura2
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
553 1 pfleura2
554 1 pfleura2
              Grad(IGeom,1)=GradTmp(4)
555 1 pfleura2
              Grad(IGeom,2)=GradTmp(7)
556 1 pfleura2
              ! We might have troubles whith angles that are not in the [0:pi] range because,
557 1 pfleura2
              ! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian...
558 1 pfleura2
              ! so that the derivative is not good, and a multiplicative factor should be added...
559 1 pfleura2
              !
560 1 pfleura2
              ! This in fact should be taken care of in B-mat calculation...
561 1 pfleura2
              !
562 1 pfleura2
              IF ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN
563 1 pfleura2
                 Grad(IGeom,3)=-1.0d0*GradTmp(8)
564 1 pfleura2
              ELSE
565 1 pfleura2
                 Grad(IGeom,3)=GradTmp(8)
566 1 pfleura2
              END IF
567 1 pfleura2
              Idx=4
568 1 pfleura2
              DO I=4,Nat
569 1 pfleura2
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
570 1 pfleura2
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
571 1 pfleura2
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
572 1 pfleura2
                 Idx=Idx+3
573 1 pfleura2
              END DO
574 1 pfleura2
              DEALLOCATE(GradTmp)
575 1 pfleura2
              ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
576 1 pfleura2
              ! but we have to convert it into internal coordinates if COORD=Baker
577 1 pfleura2
           CASE ('BAKER')
578 1 pfleura2
              ALLOCATE(GradTmp(3*Nat))
579 1 pfleura2
              GradTmp=0.d0
580 1 pfleura2
              ! Below is to be corrected.
581 1 pfleura2
              ! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN'
582 1 pfleura2
              !DO I=1, 3*Nat
583 1 pfleura2
              ! here GradTmp and Grad are gradients in Baker coordinates
584 1 pfleura2
              ! Not implemented  IF ((PROG=="VASP").AND.(RunMode=="PARA"))
585 1 pfleura2
              ! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I)
586 1 pfleura2
              !END DO
587 1 pfleura2
              Grad(IGeom,:) = GradTmp
588 1 pfleura2
              DEALLOCATE(GradTmp)
589 1 pfleura2
           CASE ('MIXED')
590 1 pfleura2
              ALLOCATE(GradTmp(3*Nat))
591 1 pfleura2
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
592 1 pfleura2
              CALL CalcBMat_mixed(nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome)
593 1 pfleura2
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
594 1 pfleura2
              DO I=1,Nat
595 1 pfleura2
                 WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I-2:3*I)
596 1 pfleura2
              END DO
597 1 pfleura2
              SELECT CASE (NCART)
598 1 pfleura2
              CASE (1)
599 1 pfleura2
                 Grad(IGeom,1:3)=GradTmp(1:3)
600 1 pfleura2
                 Grad(IGeom,4)=GradTmp(4)
601 1 pfleura2
                 Grad(IGeom,5)=GradTmp(7)
602 1 pfleura2
                 Grad(IGeom,6)=GradTmp(8)
603 1 pfleura2
                 Idx=7
604 1 pfleura2
                 IBeg=4
605 1 pfleura2
              CASE(2)
606 1 pfleura2
                 Grad(IGeom,1:3)=GradTmp(1:3)
607 1 pfleura2
                 Grad(IGeom,4:6)=GradTmp(4:6)
608 1 pfleura2
                 Grad(IGeom,7)=GradTmp(7)
609 1 pfleura2
                 Grad(IGeom,8)=GradTmp(8)
610 1 pfleura2
                 Idx=9
611 1 pfleura2
                 IBeg=4
612 1 pfleura2
              CASE DEFAULT
613 1 pfleura2
                 Idx=1
614 1 pfleura2
                 IBeg=1
615 1 pfleura2
              END SELECT
616 1 pfleura2
              DO I=IBeg,Nat
617 1 pfleura2
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
618 1 pfleura2
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
619 1 pfleura2
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
620 1 pfleura2
                 Idx=Idx+3
621 1 pfleura2
              END DO
622 1 pfleura2
              DEALLOCATE(GradTmp)
623 1 pfleura2
           CASE ("CART","HYBRID")
624 1 pfleura2
              Grad(IGeom,:)=GradCart
625 1 pfleura2
           CASE DEFAULT
626 1 pfleura2
              WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP"
627 1 pfleura2
              STOP
628 1 pfleura2
           END SELECT
629 1 pfleura2
630 1 pfleura2
           IF (debug) THEN
631 1 pfleura2
              Line='DBG Path, GradTmp'
632 1 pfleura2
              Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat)
633 1 pfleura2
           END IF
634 1 pfleura2
635 1 pfleura2
        END DO
636 1 pfleura2
637 1 pfleura2
     CASE DEFAULT
638 1 pfleura2
        WRITE(*,*) "Prog=",TRIM(Prog)," not yet supported for Para mode."
639 1 pfleura2
        STOP
640 1 pfleura2
     END SELECT
641 1 pfleura2
642 1 pfleura2
     DEALLOCATE(GradCart,PathCart)
643 1 pfleura2
644 4 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
645 4 pfleura2
!
646 4 pfleura2
! Serial executions
647 4 pfleura2
!
648 4 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
649 1 pfleura2
  ELSE ! matches IF (RunMode=="PARA") THEN
650 1 pfleura2
 ! We will launch all calculations sequentially
651 1 pfleura2
     ALLOCATE(GradTmp(NCoord))
652 1 pfleura2
     ! We have the new path, we calculate its energy and gradient
653 1 pfleura2
     IGeom0=2
654 1 pfleura2
     IGeomF=NGeomF-1
655 1 pfleura2
     IF (OptReac.OR.CalcEReac) IGeom0=1
656 1 pfleura2
     IF (OptProd.OR.CalcEProd) IGeomF=NGeomF
657 1 pfleura2
     CalcEReac=.FALSE.
658 1 pfleura2
     CalcEProd=.FALSE.
659 1 pfleura2
     ALLOCATE(GeomCart(Nat,3))
660 1 pfleura2
661 1 pfleura2
     DO IGeom=IGeom0,IGeomF
662 1 pfleura2
        WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt
663 1 pfleura2
        if (debug) WRITE(*,*) Line
664 1 pfleura2
        SELECT CASE(Coord)
665 1 pfleura2
        CASE ('ZMAT','MIXED')
666 1 pfleura2
           GeomTmp=IntCoordF(IGeom,:)
667 1 pfleura2
        CASE ('BAKER')
668 1 pfleura2
           GeomTmp=IntCoordF(IGeom,:)
669 1 pfleura2
        CASE ('CART','HYBRID')
670 1 pfleura2
!!! CAUTION : PFL 29.AUG.2008 ->
671 1 pfleura2
           ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
672 1 pfleura2
           ! This should be made  consistent !!!!!!!!!!!!!!!
673 1 pfleura2
           GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
674 1 pfleura2
           !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
675 1 pfleura2
!!! <- CAUTION : PFL 29.AUG.2008
676 1 pfleura2
        CASE DEFAULT
677 1 pfleura2
           WRITE(*,*) "Coord=",TRIM(Coord)," not recognized in EgradPath. L299"
678 1 pfleura2
           STOP
679 1 pfleura2
        END SELECT
680 1 pfleura2
681 1 pfleura2
        !if (debug) WRITE(*,*) "Calling PrintGeom"
682 1 pfleura2
        !Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
683 1 pfleura2
684 1 pfleura2
        !           IF (debug) THEN
685 1 pfleura2
        !              WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD), ' Igeom=',IGeom
686 1 pfleura2
        !              WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
687 1 pfleura2
        !              WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
688 1 pfleura2
        !              DO Iat=1,Nat
689 1 pfleura2
        !                 WRITE(*,*) XyzGeomF(IGeom,1:3,Iat)
690 1 pfleura2
        !              END DO
691 1 pfleura2
        !           END IF
692 1 pfleura2
693 1 pfleura2
        IF (debug) THEN
694 1 pfleura2
           WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD),' IGeom=',IGeom,' Iopt=',Iopt
695 1 pfleura2
           WRITE(*,*) GeomTmp(1:min(NCoord,12))
696 1 pfleura2
        END IF
697 1 pfleura2
698 1 pfleura2
        GeomCart=0.d0
699 1 pfleura2
        IF(IOpt .EQ. 1) Then
700 1 pfleura2
           Print *, 'GeomTmp='
701 1 pfleura2
           WRITE(*,'(12(1X,F6.3))') GeomTmp
702 1 pfleura2
           Print *, 'GeomCart='
703 1 pfleura2
           WRITE(*,'(12(1X,F6.3))') GeomCart
704 1 pfleura2
      END IF
705 1 pfleura2
706 1 pfleura2
       ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
707 1 pfleura2
       ! GeomCart has INTENT(OUT)
708 1 pfleura2
       Call EGrad(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart) !GeomTmp=IntCoordF(IGeom,:)
709 1 pfleura2
710 1 pfleura2
    ! Egrad calls ConvertBakerInternal_cart to convert GeomTmp=IntCoordF(IGeom,:) into
711 1 pfleura2
    ! GeomCart (Cartesian Coordinates) so that energy and gradients can be calculated
712 1 pfleura2
    ! for each geometry.
713 1 pfleura2
        XyzGeomF(IGeom,1,:)=GeomCart(:,1)
714 1 pfleura2
        XyzGeomF(IGeom,2,:)=GeomCart(:,2)
715 1 pfleura2
        XyzGeomF(IGeom,3,:)=GeomCart(:,3)
716 1 pfleura2
717 1 pfleura2
        Energies(IGeom)=E
718 1 pfleura2
719 1 pfleura2
        if (debug) THEN
720 1 pfleura2
           Line='DBG Path, GradTmp'
721 1 pfleura2
           Call PrintGeom(Line,Nat,NCoord,GRadTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
722 1 pfleura2
        END IF
723 1 pfleura2
        Grad(IGeom,:)=GradTmp
724 1 pfleura2
     END DO ! matches DO IGeom=IGeom0,IGeomF
725 1 pfleura2
     DEALLOCATE(GradTmp,GeomCart)
726 1 pfleura2
  END IF ! matches IF (RunMode=="PARA") THEN AFTER ELSE
727 1 pfleura2
728 1 pfleura2
  DEALLOCATE(GeomTmp)
729 1 pfleura2
  DEALLOCATE(x,y,z)
730 1 pfleura2
  DEALLOCATE(x_k,y_k,z_k)
731 1 pfleura2
  if (debug) Call header('Exiting EgradPath')
732 1 pfleura2
  RETURN
733 1 pfleura2
999 WRITE(*,*) "EgradPath : We should not be here !"
734 1 pfleura2
  STOP
735 1 pfleura2
736 1 pfleura2
END SUBROUTINE EGRADPATH