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