root / src / EgradPath.f90 @ 8
History  View  Annotate  Download (27.7 kB)
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=NGeomF1 
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 NONFIRST OPTIMIZATION. 
121 
! THIS WILL HAVE EFFECT ON XyzGeomF UPDATE IN BAKER CASE. 
122 
WRITE(Line,'(I5)') IGeom1 
123 
Idx=2Len_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*I2,IGeom) = XyzGeomF(IGeom,1,I) 
147 
PathCart(3*I1,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*I2,IGeom) = x(I) 
161 
PathCart(3*I1,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=NGeomF1 
223 
Grad=0.d0 
224  
225 
Line2="00" 
226 
DO IGeom=IGeom0,IGeomF 
227 
WRITE(Line,'(I5)') IGeom1 
228 
Idx=2Len_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,'TOTALFORCE')==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*Iat2: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*Iat2: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*I2: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*i2) 
303 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
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*I2: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*i2) 
351 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
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=NGeomF1 
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*I2,IGeom) = XyzGeomF(IGeom,1,I) 
399 
PathCart(3*I1,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*I2,IGeom) = x(I) 
414 
PathCart(3*I1,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*Iat2: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 Bmat 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*i2) 
571 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
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*I2: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*i2) 
619 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
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=NGeomF1 
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 