root / src / EgradPath.f90 @ 8
History  View  Annotate  Download (27.7 kB)
1 
SUBROUTINE EGradPath(Iopt) 

2  
3 
! This subroutine calculates the energy+gradient for all points of 
4 
! the path 
5  
6 
use Path_module, only : ProgExe, NCoord, PROG, Nat, NGeomF, Coord, IntCoordF, & 
7 
IndZmat, XyzGeomF, Atome, Order, OrderInv, Latr, FrozAtoms, & 
8 
Grad, Energies, renum, dzdc, massat, Pi, NCART, OptReac, & 
9 
OptProd, AtName, Unit, CalcEReac, CalcEProd 
10 
! GeomOld_all(NGeomF,NCoord), BTransInv (3*Nat,NCoord): used for Baker case 
11 
use Io_module 
12  
13 
IMPLICIT NONE 
14  
15 
INTEGER(KINT), INTENT(IN) :: Iopt 
16  
17 
INTEGER(KINT) :: IGeom,IGeom0,IGeomF 
18 
CHARACTER(LCHARS) :: Line,FileInp,RunCommand 
19 
CHARACTER(SCHARS) :: Line2 
20  
21 
REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:) ! NCoord 
22 
REAL(KREAL), ALLOCATABLE :: GradCart(:) ! 3*Nat 
23 
REAL(KREAL), ALLOCATABLE :: PathCart(:,:) ! (3*Nat,NGeomF) 
24 
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) 
25 
REAL(KREAL), ALLOCATABLE :: x_k(:), y_k(:), z_k(:) 
26 
REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3) 
27 
REAL(KREAL) :: E 
28 
LOGICAL :: Debug 
29 
LOGICAL, SAVE :: First=.TRUE. 
30  
31 
INTEGER(KINT) :: I, Iat, IBeg, Idx 
32 
REAL(KREAL) :: RTmp1, RTmp2, RTmp3 
33 
INTEGER(KINT) :: Itmp1, ITmp2 
34  
35 
INTERFACE 
36 
function valid(string) result (isValid) 
37 
CHARACTER(*), intent(in) :: string 
38 
logical :: isValid 
39 
END function VALID 
40  
41 
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old) 
42  
43 
! This routines calculates the energy E and the gradient Grad of 
44 
! a molecule with Geometry Geom (may be in internal coordinates), 
45 
! using for now, either Gaussian or Ext, more general later. 
46  
47 
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & 
48 
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, & 
49 
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & 
50 
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & 
51 
, Order,OrderInv, XPrimitiveF 
52  
53 
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) 
54 
! allocated in Path.f90 
55  
56 
use Io_module 
57  
58 
! Energy (calculated if F300K=.F., else estimated) 
59 
REAL(KREAL), INTENT (OUT) :: E 
60 
! NCoord: Number of the degrees of freedom 
61 
! IGeom: index of the geometry. 
62 
INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt 
63 
! Geometry at which gradient is calculated (cf Factual also): 
64 
REAL(KREAL), INTENT (INOUT) :: Geom(NCoord) 
65 
! Gradient calculated at Geom geometry: 
66 
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) 
67 
! Cartesian geometry corresponding to (Internal Geometry) Geom: 
68 
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) 
69 
!!! Optional, just for geometry optimization with Baker coordinates 
70 
REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3) 
71 
REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord) 
72 
! FOptGeom is a flag indicating if we are doing a geom optimization 
73 
! it can be omitted so that we use a local flag: Flag_Opt_Geom 
74 
LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom 
75 
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F 
76 
LOGICAL :: Flag_Opt_Geom 
77  
78 
END subroutine Egrad 
79  
80 
END INTERFACE 
81  
82 
debug=valid('EGradPath') 
83  
84 
if (debug) Call header("Entering EgradPath") 
85  
86 
ALLOCATE(GeomTmp(NCoord)) 
87 
ALLOCATE(x(Nat),y(Nat),z(Nat)) 
88 
ALLOCATE(x_k(Nat),y_k(Nat),z_k(Nat)) 
89  
90  
91 
IF (RunMode=="PARA") THEN ! matches at L315 
92  
93 
ALLOCATE(GradCart(3*Nat),PathCart(3*Nat,NGeomF)) 
94  
95  
96 
SELECT CASE(Prog) 
97 
CASE ("VASP") 
98 
! We will use the NEB routine of VASP to get all forces at once on Para8 queue 
99  
100 
! First, we create all the POSCAR into the 0x directories 
101  
102 
! With RunMode="Para" one cannot optimize reactants or products (due to Vasp constraints) 
103 
! so this has been done beforehand :) 
104 
IGeom0=2 
105 
IGeomF=NGeomF1 
106 
PathCart=0.d0 
107  
108 
IF (First) THEN 
109 
! For the first iteration, we do write a POSCAR into 00 and 0Final 
110 
! Vasp needs it ! 
111 
First=.FALSE. 
112 
IGeom0=1 
113 
IGeomF=NGeomF 
114 
END IF 
115  
116 
Line2="00" 
117 
DO IGeom=IGeom0,IGeomF 
118  
119 
! NOTE THAT HERE IGeom0 (IGeomF) IS NOT 1 (NGeomF) FOR NONFIRST OPTIMIZATION. 
120 
! THIS WILL HAVE EFFECT ON XyzGeomF UPDATE IN BAKER CASE. 
121 
WRITE(Line,'(I5)') IGeom1 
122 
Idx=2Len_TRIM(ADJUSTL(Line)) 
123 
FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/POSCAR" 
124 
if (debug) WRITE(*,*) "Creating ",TRIM(FileInp) 
125 
OPEN(IOTMP,File=TRIM(FileInp)) 
126  
127 
WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt 
128 
if (debug) WRITE(*,*) Line 
129 
SELECT CASE(Coord) 
130 
CASE ('ZMAT') 
131 
GeomTmp=IntCoordF(IGeom,:) 
132 
Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom)) 
133 
CASE ('BAKER') 
134 
GeomTmp=IntCoordF(IGeom,:) 
135 
! 
136 
IF (IOpt .EQ. 0) THEN 
137 
! EgradPath(...) is called only after the call of PathCreate(...) 
138 
! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always 
139 
! true if IOpt==0), THEN ExtraPol_baker(...) is called. 
140 
! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting 
141 
! the interpolated internal coordinates into cartesian coordinates. 
142 
! Thus here we don't need to reconvert the internal coordinates again 
143 
! to the cartesian coordinates. 
144 
DO I=1,Nat 
145 
PathCart(3*I2,IGeom) = XyzGeomF(IGeom,1,I) 
146 
PathCart(3*I1,IGeom) = XyzGeomF(IGeom,2,I) 
147 
PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I) 
148 
END DO 
149 
ELSE 
150 
! XyzGeomF(NGeomF,3,Nat). 
151 
x_k(:) = XyzGeomF(IGeom,1,:) 
152 
y_k(:) = XyzGeomF(IGeom,2,:) 
153 
z_k(:) = XyzGeomF(IGeom,3,:) 
154  
155 
!Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), & 
156 
! x_k,y_k,z_k,x,y,z) 
157 
! PathCart(3*Nat,NGeomF) 
158 
DO I=1,Nat 
159 
PathCart(3*I2,IGeom) = x(I) 
160 
PathCart(3*I1,IGeom) = y(I) 
161 
PathCart(3*I,IGeom) = z(I) 
162 
END DO 
163  
164 
XyzGeomF(IGeom,1,:)=x(:) 
165 
XyzGeomF(IGeom,2,:)=y(:) 
166 
XyzGeomF(IGeom,3,:)=z(:) 
167 
END IF 
168 
CASE ('MIXED') 
169 
GeomTmp=IntCoordF(IGeom,:) 
170 
Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom)) 
171 
CASE ('CART','HYBRID') 
172 
!!! CAUTION : PFL 29.AUG.2008 > 
173 
! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3) 
174 
! This should be made consistent !!!!!!!!!!!!!!! 
175 
GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 
176 
!!! < CAUTION : PFL 29.AUG.2008 
177 
PathCart(:,IGeom)=GeomTmp 
178 
CASE DEFAULT 
179 
WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L134. STOP" 
180 
STOP 
181 
END SELECT 
182  
183 
IF (debug) THEN 
184 
WRITE(*,*) "Calling PrintGeomVasp for geom",Igeom 
185 
Call PrintGeom(Line,Nat,NCoord,PathCart(1,Igeom),Coord,IoOut,Atome,Order,OrderInv,IndZmat) 
186 
END IF 
187 
Call PrintGeomVasp(Line,PathCart(1,IGeom),Latr,FrozAtoms,IoTmp) 
188  
189 
! IF (debug) THEN 
190 
! WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD), ' Igeom=',IGeom 
191 
! WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 
192 
! WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
193 
! DO Iat=1,Nat 
194 
! WRITE(*,*) XyzGeomF(IGeom,1:3,Iat) 
195 
! END DO 
196 
! END IF 
197  
198 
IF (debug) THEN 
199 
WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD),' Igeom=',IGeom,' Iopt=',Iopt 
200 
WRITE(*,*) GeomTmp(1:min(NCoord,9)) 
201 
END IF 
202  
203 
CLOSE(IOTMP) 
204 
END DO ! matches DO IGeom=IGeom0,IGeomF 
205  
206 
!!!!!!!!!!!!!!!!!!! 
207 
! 
208 
! We calculate the energies and gradients 
209 
! 
210  
211 
RunCommand=Trim(ProgExe) 
212 
if (debug) WRITE(*,*) "Launching:", TRIM(RunCommand) 
213 
call system(RunCommand) 
214 
if (debug) WRITE(*,*) "Back from:", TRIM(RunCommand) 
215  
216 
!!!!!!!!!!!!! 
217 
! 
218 
! We read the gradients and energies 
219 
! 
220 
IGeom0=2 
221 
IGeomF=NGeomF1 
222 
Grad=0.d0 
223  
224 
Line2="00" 
225 
DO IGeom=IGeom0,IGeomF 
226 
WRITE(Line,'(I5)') IGeom1 
227 
Idx=2Len_TRIM(ADJUSTL(Line)) 
228 
FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/OUTCAR" 
229 
if (debug) WRITE(*,*) "Reading E and forces from ", TRIM(FileInp) 
230 
OPEN(IOTMP,File=TRIM(FileInp)) 
231  
232 
! We first search for the forces 
233 
READ(IOTMP,'(A)',END=999,ERR=999) LINE 
234 
! We search for the last part of the OUTCAR file, after wavefunction is converged 
235 
DO WHILE (INDEX(LINE,'EDIFF is reached')==0) 
236 
READ(IOTMP,'(A)',END=999,ERR=999) LINE 
237 
END DO 
238 
READ(IOTMP,'(A)',END=999,ERR=999) LINE 
239 
DO WHILE (INDEX(LINE,'TOTEN')==0) 
240 
READ(IOTMP,'(A)',END=999,ERR=999) LINE 
241 
END DO 
242 
Line=Line(26:) 
243 
READ(LINE,*) E 
244 
if (debug) WRITE(*,'(1X,A,I5,A,F20.6)') "DBG Egrad_VASP E(",Igeom,")=",E 
245 
Energies(Igeom)=E 
246  
247 
! We search for the forces 
248 
DO WHILE (INDEX(LINE,'TOTALFORCE')==0) 
249 
READ(IOTMP,'(A)',END=999,ERR=999) LINE 
250 
END DO 
251 
READ(IOTMP,'(A)',END=999,ERR=999) LINE 
252 
DO I=1,Nat 
253 
Iat=I 
254 
IF (renum) Iat=Order(I) 
255 
READ(IOTMP,*,END=999,ERR=999) Rtmp1,RTmp2,RTmp3,GradCart(3*Iat2:3*Iat) 
256 
END DO 
257  
258 
Close(IoTmp) 
259  
260 
IF (debug) THEN 
261 
WRITE(*,*) "DBG EGRAD_VASP  GradCart=FORCES  original order" 
262 
DO I=1,Nat 
263 
Iat=Order(I) 
264 
WRITE(*,'(1X,I5,3(1X,F20.6))') I,GradCart(3*Iat2:3*Iat) 
265 
END DO 
266 
WRITE(*,*) "DBG EGRAD_VASP  GradCart=FORCES  internal order" 
267 
DO I=1,Nat 
268 
Iat=Order(I) 
269 
WRITE(*,'(3(1X,I5),3(1X,F20.6))') I,Iat,OrderInv(I),GradCart(3*I2:3*I) 
270 
END DO 
271 
END IF 
272  
273 
! VASP gives the Forces in eV/Angstrom, so we convert it into the 
274 
! gradient in ua/Angstrom 
275 
GradCart=1.0d0*ev2Au*GradCart 
276  
277 
! We convert it into the right coordinate system 
278 
! This should be put into a subroutine (because it is also in egradient) 
279 
SELECT CASE (COORD) 
280 
CASE ("ZMAT") 
281 
ALLOCATE(GradTmp(3*Nat)) 
282 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int" 
283 
CALL CalcBmat_int(nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome) 
284 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
285  
286 
Grad(IGeom,1)=GradTmp(4) 
287 
Grad(IGeom,2)=GradTmp(7) 
288 
! We might have troubles whith angles that are not in the [0:pi] range because, 
289 
! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian... 
290 
! so that the derivative is not good, and a multiplicative factor should be added... 
291 
! 
292 
! This in fact should be taken care of in calc Bmat ... 
293 
! 
294 
IF ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN 
295 
Grad(IGeom,3)=1.0d0*GradTmp(8) 
296 
ELSE 
297 
Grad(IGeom,3)=GradTmp(8) 
298 
END IF 
299 
Idx=4 
300 
DO I=4,Nat 
301 
Grad(IGeom,Idx)=GradTmp(3*i2) 
302 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
303 
Grad(IGeom,Idx+2)=GradTmp(3*i) 
304 
Idx=Idx+3 
305 
END DO 
306 
DEALLOCATE(GradTmp) 
307 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
308 
! but we have to convert it into internal coordinates if COORD=Baker 
309 
CASE ('BAKER') 
310 
ALLOCATE(GradTmp(3*Nat)) 
311 
GradTmp=0.d0 
312 
! Below is to be corrected. 
313 
! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN' 
314 
!DO I=1, 3*Nat 
315 
! here GradTmp and Grad are gradients in Baker coordinates 
316 
! Not implemented IF ((PROG=="VASP").AND.(RunMode=="PARA")) 
317 
! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I) 
318 
!END DO 
319 
Grad(IGeom,:) = GradTmp 
320 
DEALLOCATE(GradTmp) 
321 
CASE ('MIXED') 
322 
ALLOCATE(GradTmp(3*Nat)) 
323 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed" 
324 
CALL CalcBMat_mixed (nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome) 
325 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
326 
DO I=1,Nat 
327 
WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I2:3*I) 
328 
END DO 
329 
SELECT CASE (NCART) 
330 
CASE (1) 
331 
Grad(IGeom,1:3)=GradTmp(1:3) 
332 
Grad(IGeom,4)=GradTmp(4) 
333 
Grad(IGeom,5)=GradTmp(7) 
334 
Grad(IGeom,6)=GradTmp(8) 
335 
Idx=7 
336 
IBeg=4 
337 
CASE(2) 
338 
Grad(IGeom,1:3)=GradTmp(1:3) 
339 
Grad(IGeom,4:6)=GradTmp(4:6) 
340 
Grad(IGeom,7)=GradTmp(7) 
341 
Grad(IGeom,8)=GradTmp(8) 
342 
Idx=9 
343 
IBeg=4 
344 
CASE DEFAULT 
345 
Idx=1 
346 
IBeg=1 
347 
END SELECT 
348 
DO I=IBeg,Nat 
349 
Grad(IGeom,Idx)=GradTmp(3*i2) 
350 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
351 
Grad(IGeom,Idx+2)=GradTmp(3*i) 
352 
Idx=Idx+3 
353 
END DO 
354 
DEALLOCATE(GradTmp) 
355 
CASE ("CART","HYBRID") 
356 
Grad(IGeom,:)=GradCart 
357 
CASE DEFAULT 
358 
WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP" 
359 
STOP 
360 
END SELECT 
361  
362 
IF (debug) THEN 
363 
Line='DBG Path, GradTmp' 
364 
Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat) 
365 
END IF 
366 
END DO ! matches DO IGeom=IGeom0,IGeomF 
367  
368 
CASE ("GAUSSIAN") 
369 
! We create the output files. 
370 
IGeom0=2 
371 
IGeomF=NGeomF1 
372 
PathCart=0.d0 
373  
374 
If (OptReac) IGeom0=1 
375 
If (OptProd) IGeomF=NGeomF 
376  
377 
OPEN(IOTMP2,File="ListJob",STATUS="REPLACE") 
378  
379 
DO Igeom=IGeom0,IGeomF 
380 
! We first convert the geometry into a Cartesian one 
381 
SELECT CASE(Coord) 
382 
CASE ('ZMAT') 
383 
GeomTmp=IntCoordF(IGeom,:) 
384 
Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom)) 
385 
CASE ('BAKER') 
386 
GeomTmp=IntCoordF(IGeom,:) 
387 
! 
388 
IF (IOpt .EQ. 0) THEN 
389 
! EgradPath(...) is called only after the call of PathCreate(...) 
390 
! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always 
391 
! true if IOpt==0), THEN ExtraPol_baker(...) is called. 
392 
! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting 
393 
! the interpolated internal coordinates into cartesian coordinates. 
394 
! Thus here we don't need to reconvert the internal coordinates again 
395 
! to the cartesian coordinates. 
396 
DO I=1,Nat 
397 
PathCart(3*I2,IGeom) = XyzGeomF(IGeom,1,I) 
398 
PathCart(3*I1,IGeom) = XyzGeomF(IGeom,2,I) 
399 
PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I) 
400 
END DO 
401 
ELSE 
402  
403 
! XyzGeomF(NGeomF,3,Nat). 
404 
x_k(:) = XyzGeomF(IGeom,1,:) 
405 
y_k(:) = XyzGeomF(IGeom,2,:) 
406 
z_k(:) = XyzGeomF(IGeom,3,:) 
407  
408 
!Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), & 
409 
! x_k,y_k,z_k,x,y,z) 
410 
! PathCart(3*Nat,NGeomF) 
411 
DO I=1,Nat 
412 
PathCart(3*I2,IGeom) = x(I) 
413 
PathCart(3*I1,IGeom) = y(I) 
414 
PathCart(3*I,IGeom) = z(I) 
415 
END DO 
416  
417 
XyzGeomF(IGeom,1,:)=x(:) 
418 
XyzGeomF(IGeom,2,:)=y(:) 
419 
XyzGeomF(IGeom,3,:)=z(:) 
420 
END IF 
421 
CASE ('MIXED') 
422 
GeomTmp=IntCoordF(IGeom,:) 
423 
Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom)) 
424 
CASE ('CART','HYBRID') 
425 
!!! CAUTION : PFL 29.AUG.2008 > 
426 
! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3) 
427 
! This should be made consistent !!!!!!!!!!!!!!! 
428 
GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 
429 
!!! < CAUTION : PFL 29.AUG.2008 
430 
PathCart(:,IGeom)=GeomTmp 
431 
CASE DEFAULT 
432 
WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L75.STOP" 
433 
STOP 
434 
END SELECT 
435  
436 
WRITE(Line,'(I5)') IGeom 
437 
FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".com" 
438 
WRITE(IOTMP2,'(1X,A)') TRIM(FileInp) 
439  
440 
OPEN(IOTMP,File=TRIM(FileInp)) 
441 
Current => Gauss_root 
442 
DO WHILE (ASSOCIATED(Current%next)) 
443 
WRITE(IOTMP,'(1X,A)') Trim(current%line) 
444 
WRITE(*,'(1X,A)') Trim(current%line) 
445 
Current => current%next 
446 
END DO 
447  
448 
WRITE(IOTMP,*) 
449  
450 
Current => Gauss_comment 
451 
DO WHILE (ASSOCIATED(Current%next)) 
452 
WRITE(IOTMP,'(1X,A)') Trim(current%line) 
453 
Current => current%next 
454 
END DO 
455 
! WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)"),Igeom,NgeomF,Iopt 
456  
457 
WRITE(IOTMP,*) 
458 
WRITE(IOTMP,*) Trim(Gauss_charge) 
459  
460 
DO I=1,Nat 
461 
If (renum) THEN 
462 
Iat=Order(I) 
463 
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)) & 
464 
,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom) & 
465 
,TRIM(Gauss_Paste(I)) 
466 
ELSE 
467 
Iat=OrderInv(I) 
468 
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)) & 
469 
,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom) & 
470 
, TRIM(Gauss_Paste(Iat)) 
471 
END IF 
472 
END DO 
473  
474 
WRITE(IOTMP,*) 
475  
476 
Current => Gauss_End 
477 
DO WHILE (ASSOCIATED(Current%next)) 
478 
WRITE(IOTMP,'(1X,A)') Trim(current%line) 
479 
Current => current%next 
480 
END DO 
481  
482 
WRITE(IOTMP,*) 
483 
CLOSE(IOTMP) 
484  
485 
END DO 
486  
487 
Close(IOTMP2) 
488  
489 
! We launch the parallel calculations 
490  
491 
Line=Trim(ProgExe) 
492 
IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(Line) 
493 
call system(Line) 
494 
IF (DEBUG) WRITE(*,*)'Back from RunCommand:' 
495  
496 
! We read energies and gradients from the log files 
497  
498 
DO IGeom=IGeom0,IGeomF 
499 
WRITE(Line,'(I5)') IGeom 
500 
FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".log" 
501 
OPEN(IOTMP,File=TRIM(FileInp)) 
502  
503 
! We first search for the forces 
504 
READ(IOTMP,'(A)') LINE 
505 
DO WHILE (INDEX(LINE,'Forces (Hartrees/Bohr)')==0) 
506 
! here, we have a problem, because there are so many ways to write E for Gaussian :) 
507 
! but we do not really need E except if we want to weight the different points 
508 
! on the path... that will be done later :) 
509 
! And I should use ConvGaussXmol ... 
510 
! PFL 3rd Jan 2008 > 
511 
! I have finally upgraded the energy reading phase so that it should be able to read 
512 
! many formats corresponding to HF, DFT, MP2, AM1, PM3, ONIOM with or w/o PCM 
513 
! 
514 

515 
! For MP2, energy is after the last = 
516 
IF ((Line(2:9)=="E2 = ")) THEN 
517 
Itmp1=Index(LINE,"=",BACK=.TRUE.)+1 
518 
READ(LINE(Itmp1:),*) e 
519 
END IF 
520 
! For other methods, it is after the first = 
521 
IF ((LINE(2:9)=='Energy= ').OR.(Line(2:9)=="SCF Done").OR.(Line(2:9)=="ONIOM: e") & 
522 
.OR.(Line(2:9)==" with al") & 
523 
.OR.(Line(2:31)=="Counterpoise: corrected energy")) THEN 
524 
Itmp1=Index(LINE,"=")+1 
525 
READ(LINE(Itmp1:),*) e 
526 
END IF 
527 
! < PFL 3 Jan 2008 
528 
READ(IOTMP,'(A)') LINE 
529 
END DO 
530 
READ(IOTMP,'(A)') LINE 
531 
READ(IOTMP,'(A)') LINE 
532 
DO I=1,Nat 
533 
Iat=I 
534 
IF (renum) Iat=Order(I) 
535 
READ(IOTMP,*) Itmp1,ITmp2, GradCart(3*Iat2:3*Iat) 
536 
END DO 
537  
538 
! Gaussian gives the Forces in ua/a0, so we convert it into the 
539 
! gradient in ua/Angstrom 
540 
gradCart=1.0d0*Unit*gradCart 
541  
542 
CLOSE(IOTMP) 
543 
Energies(IGeom)=E 
544 

545 
! We convert it into the right coordinate system 
546 
! This should be put into a subroutine (because it is also in egradient) 
547 
SELECT CASE (COORD) 
548 
CASE ("ZMAT") 
549 
ALLOCATE(GradTmp(3*Nat)) 
550 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int" 
551 
CALL CalcBMat_int (nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome) 
552 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
553  
554 
Grad(IGeom,1)=GradTmp(4) 
555 
Grad(IGeom,2)=GradTmp(7) 
556 
! We might have troubles whith angles that are not in the [0:pi] range because, 
557 
! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian... 
558 
! so that the derivative is not good, and a multiplicative factor should be added... 
559 
! 
560 
! This in fact should be taken care of in Bmat calculation... 
561 
! 
562 
IF ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN 
563 
Grad(IGeom,3)=1.0d0*GradTmp(8) 
564 
ELSE 
565 
Grad(IGeom,3)=GradTmp(8) 
566 
END IF 
567 
Idx=4 
568 
DO I=4,Nat 
569 
Grad(IGeom,Idx)=GradTmp(3*i2) 
570 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
571 
Grad(IGeom,Idx+2)=GradTmp(3*i) 
572 
Idx=Idx+3 
573 
END DO 
574 
DEALLOCATE(GradTmp) 
575 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
576 
! but we have to convert it into internal coordinates if COORD=Baker 
577 
CASE ('BAKER') 
578 
ALLOCATE(GradTmp(3*Nat)) 
579 
GradTmp=0.d0 
580 
! Below is to be corrected. 
581 
! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN' 
582 
!DO I=1, 3*Nat 
583 
! here GradTmp and Grad are gradients in Baker coordinates 
584 
! Not implemented IF ((PROG=="VASP").AND.(RunMode=="PARA")) 
585 
! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I) 
586 
!END DO 
587 
Grad(IGeom,:) = GradTmp 
588 
DEALLOCATE(GradTmp) 
589 
CASE ('MIXED') 
590 
ALLOCATE(GradTmp(3*Nat)) 
591 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed" 
592 
CALL CalcBMat_mixed(nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome) 
593 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
594 
DO I=1,Nat 
595 
WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I2:3*I) 
596 
END DO 
597 
SELECT CASE (NCART) 
598 
CASE (1) 
599 
Grad(IGeom,1:3)=GradTmp(1:3) 
600 
Grad(IGeom,4)=GradTmp(4) 
601 
Grad(IGeom,5)=GradTmp(7) 
602 
Grad(IGeom,6)=GradTmp(8) 
603 
Idx=7 
604 
IBeg=4 
605 
CASE(2) 
606 
Grad(IGeom,1:3)=GradTmp(1:3) 
607 
Grad(IGeom,4:6)=GradTmp(4:6) 
608 
Grad(IGeom,7)=GradTmp(7) 
609 
Grad(IGeom,8)=GradTmp(8) 
610 
Idx=9 
611 
IBeg=4 
612 
CASE DEFAULT 
613 
Idx=1 
614 
IBeg=1 
615 
END SELECT 
616 
DO I=IBeg,Nat 
617 
Grad(IGeom,Idx)=GradTmp(3*i2) 
618 
Grad(IGeom,Idx+1)=GradTmp(3*i1) 
619 
Grad(IGeom,Idx+2)=GradTmp(3*i) 
620 
Idx=Idx+3 
621 
END DO 
622 
DEALLOCATE(GradTmp) 
623 
CASE ("CART","HYBRID") 
624 
Grad(IGeom,:)=GradCart 
625 
CASE DEFAULT 
626 
WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP" 
627 
STOP 
628 
END SELECT 
629  
630 
IF (debug) THEN 
631 
Line='DBG Path, GradTmp' 
632 
Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat) 
633 
END IF 
634  
635 
END DO 
636  
637 
CASE DEFAULT 
638 
WRITE(*,*) "Prog=",TRIM(Prog)," not yet supported for Para mode." 
639 
STOP 
640 
END SELECT 
641  
642 
DEALLOCATE(GradCart,PathCart) 
643  
644 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
645 
! 
646 
! Serial executions 
647 
! 
648 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
649 
ELSE ! matches IF (RunMode=="PARA") THEN 
650 
! We will launch all calculations sequentially 
651 
ALLOCATE(GradTmp(NCoord)) 
652 
! We have the new path, we calculate its energy and gradient 
653 
IGeom0=2 
654 
IGeomF=NGeomF1 
655 
IF (OptReac.OR.CalcEReac) IGeom0=1 
656 
IF (OptProd.OR.CalcEProd) IGeomF=NGeomF 
657 
CalcEReac=.FALSE. 
658 
CalcEProd=.FALSE. 
659 
ALLOCATE(GeomCart(Nat,3)) 
660  
661 
DO IGeom=IGeom0,IGeomF 
662 
WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt 
663 
if (debug) WRITE(*,*) Line 
664 
SELECT CASE(Coord) 
665 
CASE ('ZMAT','MIXED') 
666 
GeomTmp=IntCoordF(IGeom,:) 
667 
CASE ('BAKER') 
668 
GeomTmp=IntCoordF(IGeom,:) 
669 
CASE ('CART','HYBRID') 
670 
!!! CAUTION : PFL 29.AUG.2008 > 
671 
! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3) 
672 
! This should be made consistent !!!!!!!!!!!!!!! 
673 
GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 
674 
! GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
675 
!!! < CAUTION : PFL 29.AUG.2008 
676 
CASE DEFAULT 
677 
WRITE(*,*) "Coord=",TRIM(Coord)," not recognized in EgradPath. L299" 
678 
STOP 
679 
END SELECT 
680  
681 
!if (debug) WRITE(*,*) "Calling PrintGeom" 
682 
!Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat) 
683  
684 
! IF (debug) THEN 
685 
! WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD), ' Igeom=',IGeom 
686 
! WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 
687 
! WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/)) 
688 
! DO Iat=1,Nat 
689 
! WRITE(*,*) XyzGeomF(IGeom,1:3,Iat) 
690 
! END DO 
691 
! END IF 
692  
693 
IF (debug) THEN 
694 
WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD),' IGeom=',IGeom,' Iopt=',Iopt 
695 
WRITE(*,*) GeomTmp(1:min(NCoord,12)) 
696 
END IF 
697  
698 
GeomCart=0.d0 
699 
IF(IOpt .EQ. 1) Then 
700 
Print *, 'GeomTmp=' 
701 
WRITE(*,'(12(1X,F6.3))') GeomTmp 
702 
Print *, 'GeomCart=' 
703 
WRITE(*,'(12(1X,F6.3))') GeomCart 
704 
END IF 
705 

706 
! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT). 
707 
! GeomCart has INTENT(OUT) 
708 
Call EGrad(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart) !GeomTmp=IntCoordF(IGeom,:) 
709  
710 
! Egrad calls ConvertBakerInternal_cart to convert GeomTmp=IntCoordF(IGeom,:) into 
711 
! GeomCart (Cartesian Coordinates) so that energy and gradients can be calculated 
712 
! for each geometry. 
713 
XyzGeomF(IGeom,1,:)=GeomCart(:,1) 
714 
XyzGeomF(IGeom,2,:)=GeomCart(:,2) 
715 
XyzGeomF(IGeom,3,:)=GeomCart(:,3) 
716 

717 
Energies(IGeom)=E 
718  
719 
if (debug) THEN 
720 
Line='DBG Path, GradTmp' 
721 
Call PrintGeom(Line,Nat,NCoord,GRadTmp,Coord,6,Atome,Order,OrderInv,IndZmat) 
722 
END IF 
723 
Grad(IGeom,:)=GradTmp 
724 
END DO ! matches DO IGeom=IGeom0,IGeomF 
725 
DEALLOCATE(GradTmp,GeomCart) 
726 
END IF ! matches IF (RunMode=="PARA") THEN AFTER ELSE 
727  
728 
DEALLOCATE(GeomTmp) 
729 
DEALLOCATE(x,y,z) 
730 
DEALLOCATE(x_k,y_k,z_k) 
731 
if (debug) Call header('Exiting EgradPath') 
732 
RETURN 
733 
999 WRITE(*,*) "EgradPath : We should not be here !" 
734 
STOP 
735  
736 
END SUBROUTINE EGRADPATH 