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 re-parameterize 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