Révision 7 src/Path.f90
Path.f90 (revision 7)  

212  212 
! We add Siesta (and start thinking of add CP2K) 
213  213 
! as a new "energy engine" code 
214  214 
! TO DO to go to 1.5: finish cleaning... 
215 
! 

215  216 
! 2013 Feb Opt'n Path v1.48 
216  217 
! We add some keyword for more output for Geometry Optimizations 
217  218 
! GeomFile: name of the file to print out the geometries and their energy 
218  219 
! as Xyz file. (only format for now) 
220 
! 

221 
! 2013 Feb opt'n Path v1.49 

222 
! We add the possibility to analyze the geometries in terms 

223 
! of Bonds, Angles, dihedrals 

224 
! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path) 

225 
! This is basicaly the merging of AnaPath into Opt'n Path... at last ! 

226 
! This is done by adding a Flag into the Path namelist: 

227 
! AnaGeom: logical, if true Path looks for the AnaList namelist 

228 
! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path 

229 
! will read it and print the values of the variable in a Path.dat file 

230 
! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. 

219  231  
220  232 
PROGRAM PathOpt 
221  233  
...  ...  
249  261 
LOGICAL :: FChkFrozen 
250  262  
251  263 
! Energies for all points along the path at the previous iteration 
252 
REAL(KREAL), ALLOCATABLE :: EPathold(:) ! NGeomF, where it is deallocated: Prakash


264 
REAL(KREAL), ALLOCATABLE :: EPathold(:) ! NGeomF 

253  265 
! Maximum step allowed for this geometry 
254 
REAL(KREAL), ALLOCATABLE :: MaxStep(:) ! NGeomF, where it is deallocated: Prakash


266 
REAL(KREAL), ALLOCATABLE :: MaxStep(:) ! NGeomF 

255  267  
256  268 
! these are used to read temporary coordinates 
257  269 
LOGICAL :: FFrozen,FCart 
...  ...  
272  284  
273  285 
! INTEGER(KINT), EXTERNAL :: Iargc 
274  286 
INTEGER(KINT), EXTERNAL :: ConvertNumAt 
287 
INTEGER(KINT) :: IoS 

275  288  
276  289 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
277  290 
! 
...  ...  
327  340 
OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, & 
328  341 
FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, & 
329  342 
NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, & 
330 
NGintMax, Align, CalcEReac,CalcEProd, GeomFile 

343 
NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile


331  344  
332  345 
Namelist/cartlist/list 
333  346 
Namelist/frozenlist/list 
347 
Namelist/analist/nb 

334  348  
335  349  
336  350 
Flag_Opt_Geom = .FALSE. ! Initialized. 
...  ...  
368  382 
WRITE(*,*) "Input: string that indicates the type of the input geometries." 
369  383 
WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp" 
370  384 
WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations." 
371 
WRITE(*,*) " Accepted values are: Gaussian, Vasp, Mopac, Test, Siesta or Ext"


372 
WRITE(*,*) " In case of a Gaussian or Mopac calculations, input must be set to Cart." 

385 
WRITE(*,*) " Accepted values are: Gaussian, Vasp, Mopac, TurboMole, Siesta, Test or Ext"


386 
WRITE(*,*) " * In case of a Gaussian or Mopac calculations, input must be set to Cart."


373  387 
WRITE(*,*) " One example of a gaussian/mopac input should be added at the end of the" 
374 
WRITE(*,*) " input file.See example file Test_G03.path or Test_Mopac.path"


375 
WRITE(*,*) " In the case of a VASP calculation, if input is set to Cart, then" 

376 
WRITE(*,*) " the preamble of a VASP calculation should be added at the end of" 

377 
WRITE(*,*) " the input file. See example file Test_VASP_cart.path" 

378 
WRITE(*,*) " In the case of a VASP calculation, one should also give value of the " 

388 
WRITE(*,*) " input file. See example file Test_G03.path or Test_Mopac.path"


389 
WRITE(*,*) " * In the case of a VASP calculation, if input is set to Cart, then"


390 
WRITE(*,*) " the preamble of a VASP calculation should be added at the end of"


391 
WRITE(*,*) " the input file. See example file Test_VASP_cart.path"


392 
WRITE(*,*) " In the case of a VASP calculation, one should also give value of the "


379  393 
WRITE(*,*) " Runmode variable" 
394 
WRITE(*,*) " * In the case of a SIESTA calculation, an example of a Siesta input file" 

395 
WRITE(*,*) " should be added at the end of the input file. See Test_Siesta.path" 

380  396 
WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy" 
381  397 
WRITE(*,*) " and gradient of the whole path or not. If one wants to use VASP," 
382  398 
WRITE(*,*) " Runmode must be set to PARA." 
...  ...  
413  429 
WRITE(*,*) "Smax: Maximum length of a step during path optimization" 
414  430 
WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary" 
415  431 
WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken" 
416 
WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"


432 
WRITE(*,*) "FTan: When moving the path, this gives the proportion of the displacement tangent to the path"


417  433 
WRITE(*,*) " that is kept. FTan=1. corresponds to the full displacement, " 
418 
WRITE(*,*) " whereas FTan=0. gives a displacement orthogonal to the path." 

434 
WRITE(*,*) " whereas FTan=0. gives a displacement orthogonal to the path. Default: Ftan=0."


419  435 
WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization." 
420  436 
WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path." 
421  437 
WRITE(*,*) " This option indicates the first step where spline interpolation is used." 
...  ...  
451  467 
WRITE(*,*) " If True, a &cartlist namelist containing the list of cart atoms must be given." 
452  468 
WRITE(*,*) " By default, only frozen atoms are described in cartesian coordinates." 
453  469 
WRITE(*,*) "" 
470 
WRITE(*,*) "AnaGEom: If TRUE, Opt'n Path will create a file .dat for geometries analysis." 

471 
WRITE(*,*) " If True, Opt'n Path will look for the AnaList namelist after the Path Namelist." 

472 
WRITE(*,*) "AnaList: list of variables for geometry analysis. If present and if AnaGeom=T then" 

473 
WRITE(*,*) " Opt'n Path will read it and print the values of the variable in a .dat file" 

474 
WRITE(*,*) " If AnaGeom is T but Analist is not given, then Path.dat just contains the energy." 

475 
WRITE(*,*) "Analist contains the number of variables to monitor: Nb, and is followed by the description" 

476 
WRITE(*,*) "of the variables among:" 

477 
WRITE(*,*) "b(ond) At1 At2" 

478 
WRITE(*,*) "a(ngle) At1 At2 At3" 

479 
WRITE(*,*) "d(ihedral) At1 At2 At3 At4" 

480 
WRITE(*,*) "c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added" 

481 
WRITE(*,*) "at the end of the real atoms of the system" 

482 
WRITE(*,*) "" 

454  483 
WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry." 
455  484 
WRITE(*,*) " If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. " 
456  485 
WRITE(*,*) " It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]" 
...  ...  
470  499 
END SELECT 
471  500  
472  501  
473  
474  502 
! We initiliaze variables 
475  503 
Pi=dacos(1.0d0) 
476  504 
Nat=1 
...  ...  
507  535 
CalcEReac=.FALSE. 
508  536 
CalcEProd=.FALSE. 
509  537 
EReac=0.d0 
510 
EProd=0.d0 

538 
EProd=0.d0


511  539 
OptGeom=1 
512  540 
GeomFile="EMPTY" 
541 
AnaGeom=.FALSE. 

542 
AnaFile="EMPTY" 

543 
Nb=0 

544 
NbCom=0 

513  545 
PathOnly=.FALSE. 
514  546 
Prog='EMPTY' 
515  547 
ProgExe='EMPTY' 
...  ...  
580  612 
END IF 
581  613 
END IF 
582  614  
583 
If (COORD.EQ.'XYZ') THEN


615 
If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN


584  616 
COORD='CART' 
585  617 
END IF 
586  618  
...  ...  
651  683 
WRITE(*,*) "Input has been set to the default: ",INPUT 
652  684 
END IF 
653  685  
654  
655  
656  686 
WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP")) 
657  687  
658  688 
IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN 
...  ...  
715  745 
if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart 
716  746 
if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen 
717  747  
748 
if (AnaGeom) THEN 

749 
REWIND(IOIN) 

750 
READ(IOIN,AnaList,IOSTAT=IoS) 

751 
IF (IOS==0) THEN ! we have a AnaList namelist 

752 
Call ReadAnaList 

753 
END IF 

754 
IF (AnaFile=="EMPTY") THEN 

755 
AnaFile=TRIM(PathName) // '.dat' 

756 
END IF 

757 
OPEN(IODat,File=AnaFile) 

758 
END IF 

759  
760 


761  
718  762 
if (NMaxPtPath.EQ.1) then 
719  763 
NMaxPtPath=min(50*NGeomF,2000) 
720  764 
NMaxPtPath=max(20*NGeomF,NMaxPtPath) 
...  ...  
752  796 
ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat)) 
753  797 
ALLOCATE(MassAt(NAt),OrderInv(Nat)) 
754  798  
799 
AtName="" 

800 
MassAt=1. 

801  
755  802 
! We read the initial geometries 
756  803 
Call Read_Geom(input) 
757  804  
...  ...  
761  808 
Atome(I)=ConvertNumAt(AtName(I)) 
762  809 
END DO 
763  810  
764 
!*********** HERE ************* 

765 
!* To DO: 

766 
!* deplacer le test sur les frozen dans une sous routine 

767 
!* deplacer la lecture de sinput dans une sous routine 

768 
!* creer lecture input siesta 

769 
!* creer egrad siesta 

770 
!* 

771 
!************************************* 

811 
If (AnaGeom) THEN 

812 
! We print the description of the varialbes in AnaFile 

813 
Call PrintAnaList(IoDat) 

814 
END IF 

772  815  
773  816 
! PFL 23 Sep 2008 > 
774  817 
! We check that frozen atoms are really frozen, ie that their coordinates do not change 
...  ...  
1129  1172  
1130  1173 
if (debug) WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt 
1131  1174 
Call Write_path(IOpt) 
1175 
! Shall we analyze the geometries ? 

1176 
IF (AnaGeom) Call AnaPath(Iopt,IoDat) 

1177  
1132  1178 
if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp) Call Write_vasp(poscar) 
1133  1179  
1134  1180 
! We have the geometries, the first gradients... we will generate the first hessian matrices :) 
...  ...  
1803  1849 
GradOld=Grad 
1804  1850 
EPathOld=Energies 
1805  1851  
1852  
1853  
1806  1854 
! Shall we reparameterize the path ? 
1807  1855 
! PFL 2007/Apr/10 > 
1808  1856 
! We call PathCreate in all cases, it will take care of the 
...  ...  
1902  1950 
END DO ! matches DO IGeom=1,NGeomF 
1903  1951  
1904  1952 
Call Write_path(IOpt) 
1953 
! Shall we analyze the geometries ? 

1954 
IF (AnaGeom) Call AnaPath(Iopt,IoDat) 

1905  1955  
1906  1956 
! We update the Hessian matrices 
1907  1957 
IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord 
...  ...  
2048  2098  
2049  2099 
WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt 
2050  2100  
2101 
Close(IOIN) 

2102 
Close(IOOUT) 

2103 
IF (AnaGeom) Close(IODAT) 

2104  
2051  2105 
END PROGRAM PathOpt 
Formats disponibles : Unified diff