Revision 8

src/Egrad.f90 (revision 8)
4 4
  ! a molecule with Geometry Geom (may be in internal coordinates),
5 5
  ! using for now, either Gaussian or Ext, more general later.
6 6

  
7
  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
8
       prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, renum, &
7
  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Atome,massat, &
8
       prog,NCart,XyzGeomF, XyzGeomI, renum, &
9 9
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
10
       , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
11
      , Order,OrderInv, XPrimitiveF
10
       , BprimT,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
11
      , OrderInv, XPrimitiveF
12 12

  
13 13
  ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
14 14
  ! allocated in Path.f90
src/Hupdate_all.f90 (revision 8)
63 63
       character(*), intent(in) :: string
64 64
     end function valid
65 65

  
66

  
67
    SUBROUTINE die(routine, msg, file, line, unit)
68

  
69
      Use VarTypes
70
      Use io_module
71

  
72
      implicit none
73
!--------------------------------------------------------------- Input Variables
74
      character(len=*), intent(in)           :: routine, msg
75
      character(len=*), intent(in), optional :: file
76
      integer(KINT), intent(in), optional      :: line, unit
77

  
78
    END SUBROUTINE die
79

  
80
    SUBROUTINE Warning(routine, msg, file, line, unit)
81

  
82
      Use VarTypes
83
      Use io_module
84

  
85
      implicit none
86

  
87
      character(len=*), intent(in)           :: routine, msg
88
      character(len=*), intent(in), optional :: file
89
      integer(KINT), intent(in), optional      :: line, unit
90

  
91
    END SUBROUTINE Warning
92

  
93

  
66 94
  end interface
67 95

  
68 96
  ! ======================================================================
......
85 113
           CASE ("DFP")
86 114
              Call Hinvup_DFP(n,ds,dgrad,hess)
87 115
           CASE DEFAULT
88
              WRITE(*,*) Trim(HUpdate) // " not known: using BFGS"
116
              Call Warning("Hupdate_all Hinv=T",Trim(HUpdate) // " not known: using BFGS" &
117
                   ,Unit=IoOut)
89 118
              Call Hinvup_BFGS(n,ds,dgrad,hess)
90 119
         END SELECT
91 120
     CASE (.FALSE.)
......
101 130
! we use the fact that DFP is the dual of BFGS
102 131
              Call Hinvup_DFP(n,dgrad,ds,hess)
103 132
           CASE DEFAULT
104
              WRITE(*,*) Trim(HUpdate) // " not known: using Bofill"
133
              Call Warning("Hupdate_all Hinv=F",Trim(HUpdate) // " not known: using BoFill" &
134
                   ,Unit=IoOut)
105 135
              Call Hupdate_Bofill(n,ds,dgrad,hess)
106 136
           END SELECT
107 137
    END SELECT
src/Path_module.f90 (revision 8)
12 12
  INTEGER(KINT), PARAMETER :: MaxFroz=100
13 13
  REAL(KREAL), PARAMETER :: a0=0.529177249d0
14 14
  REAL(KREAL), PARAMETER :: Unit=1.d0/a0
15
  REAL(KREAL), PARAMETER :: au2kcal=627.509608d0
16 15
  REAL(KREAL) :: Pi
17 16

  
18 17
 ! Frozen contains the indices of frozen atoms
......
72 71
  LOGICAL :: AnaGeom
73 72
! Name of the file to analyse geometries
74 73
  CHARACTER(LCHARS) :: AnaFile
74
! Name of the Gplot file to see the path evolution
75
  CHARACTER(LCHARS) :: GplotFile
75 76
! Nb: number of variables to monitor, including Centers of Mass
76 77
  INTEGER(KINT) :: Nb
77 78
! NbVar: number of geometrical variables to monitor
src/Step_DIIS_all.f90 (revision 8)
141 141
            Vfree(:,I)=Vfree(:,I)-Norm*Tangent
142 142
     END DO
143 143

  
144
         Idx=0.
144
         Idx=0
145 145
         ! Schmidt orthogonalization of the Vfree vectors
146 146
         DO I=1,NCoord
147 147
            ! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable.
src/Makefile (revision 8)
253 253
      SearchInput.f90 \
254 254
      CleanString.f90 \
255 255
      InString.f90 \
256
      NoComment.f90 \
257
      NoString.f90 \
256 258
      Egrad.f90  \
257 259
      EgradPath.f90 \
258 260
      egrad_gaussian.f90 \
src/ReadInput_gaussian.f90 (revision 8)
16 16
  END INTERFACE
17 17

  
18 18

  
19
  CHARACTER(132) ::  Line,Line2
19
  CHARACTER(132) ::  Line
20 20
  INTEGER(KINT) :: LineL, Idx, Iat
21 21
  INTEGER(KINT) :: I
22
  REAL(KREAL) :: Xtmp,YTmp,ZTmp  
23
  LOGICAL :: Tchk
24 22
  
25 23
  LOGICAL :: Debug
26 24

  
src/Step_GEDIIS_All.f90 (revision 8)
90 90
            Vfree(:,I)=Vfree(:,I)-Norm*Tangent
91 91
     END DO
92 92

  
93
         Idx=0.
93
         Idx=0
94 94
         ! Schmidt orthogonalization of the Vfree vectors
95 95
         DO I=1,NCoord
96 96
            ! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable.
src/Path.f90 (revision 8)
331 331
     END subroutine Hinvup_BFGS_new
332 332

  
333 333

  
334
    FUNCTION InString(Line,String,Case,Clean,Back)  Result(Pos)
335

  
336
      Use VarTypes
337

  
338
      implicit none
339
! Input
340
      CHARACTER(*), INTENT(IN) :: Line
341
      CHARACTER(*), INTENT(IN) :: String
342
      LOGICAL, OPTIONAL, INTENT(IN) :: Case
343
      LOGICAL, OPTIONAL, INTENT(IN) :: Back
344
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean
345

  
346
! Output
347
! the position of String in Line (the first one) unless Back is present
348
      INTEGER(KINT) :: Pos
349
    END FUNCTION InString
350

  
351
    SUBROUTINE die(routine, msg, file, line, unit)
352

  
353
      Use VarTypes
354
      Use io_module
355

  
356
      implicit none
357
      character(len=*), intent(in)           :: routine, msg
358
      character(len=*), intent(in), optional :: file
359
      integer(KINT), intent(in), optional      :: line, unit
360

  
361
    END SUBROUTINE die
362

  
363
    SUBROUTINE Warning(routine, msg, file, line, unit)
364

  
365
      Use VarTypes
366
      Use io_module
367

  
368
      implicit none
369

  
370
      character(len=*), intent(in)           :: routine, msg
371
      character(len=*), intent(in), optional :: file
372
      integer(KINT), intent(in), optional      :: line, unit
373

  
374
    END SUBROUTINE Warning
375

  
376

  
334 377
  END INTERFACE
335 378

  
336 379
  Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, &
......
340 383
       OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, &
341 384
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
342 385
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
343
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile
386
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile,GplotFile
344 387

  
345 388
  Namelist/cartlist/list
346 389
  Namelist/frozenlist/list
......
355 398
  !
356 399
  WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013"
357 400

  
358
  ! We read the files names
401
  ! We read the files  names
359 402
!  SELECT CASE (iargc())
360 403
  SELECT CASE (command_argument_count())
361 404
  CASE (2)
......
540 583
  GeomFile="EMPTY"
541 584
  AnaGeom=.FALSE.
542 585
  AnaFile="EMPTY"
586
  GplotFile="EMPTY"
543 587
  Nb=0
544 588
  NbCom=0
545 589
  PathOnly=.FALSE.
......
638 682
     if (CalcName=="EMPTY") CalcName="Path"
639 683
     if (ISuffix=="EMPTY")  ISuffix="com"
640 684
     if (OSuffix=="EMPTY")  OSuffix="log"
685
     UnitProg="au"
686
     ConvE=au2kcal
641 687
    CASE ("G03")
642 688
     Prog='GAUSSIAN'
643 689
     If (ProgExe=="EMPTY") ProgExe="g03"
644 690
     if (CalcName=="EMPTY") CalcName="Path"
645 691
     if (ISuffix=="EMPTY")  ISuffix="com"
646 692
     if (OSuffix=="EMPTY")  OSuffix="log"
693
     UnitProg="au"
694
     ConvE=au2kcal
647 695
    CASE ("EXT")
648 696
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
649 697
     if (CalcName=="EMPTY") CalcName="Ext"
650 698
     if (ISuffix=="EMPTY")  ISuffix="in"
651 699
     if (OSuffix=="EMPTY")  OSuffix="out"
700
     UnitProg="au"
701
     ConvE=au2kcal
652 702
    CASE ('MOPAC')
653 703
     If (ProgExe=="EMPTY") ProgExe="mopac"
654 704
     if (CalcName=="EMPTY") CalcName="Path"
655 705
     if (ISuffix=="EMPTY")  ISuffix="mop"
656 706
     if (OSuffix=="EMPTY")  OSuffix="out"
707
     UnitProg="au"
708
     ConvE=au2kcal
657 709
    CASE ("SIESTA")
658 710
     If (ProgExe=="EMPTY") ProgExe="siesta"
659 711
     if (CalcName=="EMPTY") CalcName="Path"
660 712
     if (ISuffix=="EMPTY")  ISuffix="fdf"
661 713
     if (OSuffix=="EMPTY")  OSuffix="out"
714
     UnitProg="eV"
715
     ConvE=eV2kcal
662 716
    CASE ("VASP")
663 717
     If (ProgExe=="EMPTY") ProgExe="VASP"
718
     UnitProg="eV"
719
     ConvE=eV2kcal
664 720
    CASE ("TM","TURBOMOLE")
665 721
     Prog="TURBOMOLE"
666 722
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
723
     UnitProg="au"
724
     ConvE=au2kcal
667 725
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
668 726
     ProgExe=""
727
     UnitProg="au"
728
     ConvE=au2kcal
669 729
    CASE DEFAULT
670 730
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
671 731
     STOP
......
750 810
     READ(IOIN,AnaList,IOSTAT=IoS)
751 811
     IF (IOS==0) THEN ! we have a AnaList namelist 
752 812
        Call ReadAnaList
813
     ELSE
814
! No AnaList namelist, we will print only the energy
815
        FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)'
753 816
     END IF
754 817
     IF (AnaFile=="EMPTY") THEN
755 818
        AnaFile=TRIM(PathName) // '.dat'
756 819
     END IF
757 820
     OPEN(IODat,File=AnaFile)
821
     IF (GplotFile=="EMPTY") THEN
822
        GplotFile=TRIM(PathName) // '.gplot'
823
     END IF
824
     
825
     OPEN(IODat,File=AnaFile)
826
     If (OptGeom<=0) THEN
827
        OPEN(IOGplot,File=GPlotFile)
828
        WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot -persist"
829
        WRITE(IoGplot,'(A)') " set pointsize 2"
830
        WRITE(IoGplot,'(A)') " xcol=1"
831
        WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3
832
     END IF
758 833
  END IF
759 834

  
760 835
        
......
799 874
  AtName=""
800 875
  MassAt=1.
801 876

  
877
! As we have used rewind many times, the position in the input file
878
! is blurry (at least !) we thus have to go back to the Geometries description
879
! (TO DO: use a better Input parser !)
880

  
881
  REWIND(IOIN)
882
! We search the '&path' namelist
883
  READ(IOIN,'(A)',Iostat=Ios) Line
884
  Fini=.FALSE.
885
  DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0)
886
     if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line)
887
     READ(IOIN,'(A)',Iostat=Ios) Line
888
  END DO
889
  if (debug) WRITE(*,*) "Path found:",TRIM(LINE)
890
  ! We are on the &path line, we move to the end
891
  Call NoComment(Line)
892
  DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
893
     if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
894
     READ(IOIN,'(A)',Iostat=Ios) Line
895
     Call NoComment(Line)
896
     Call noString(Line)
897
  END DO
898
! We are at then end of &Path / namelist
899
! We scan for other Namelists
900
  READ(IOIN,'(A)',Iostat=Ios) Line
901
  if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
902
  DO WHILE  ((IoS==0).AND.(Index(Line,'&')/=0))
903
     Call NoComment(Line)
904
     if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
905
     Idx=InString(Line,'Analist')
906
     DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
907
        if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
908
        READ(IOIN,'(A)',Iostat=Ios) Line
909
        Call NoComment(Line)
910
     END DO
911
     If (Idx/=0) THEN
912
! we have just read &Analist namelist, we have to skip the variables description
913
        DO I=1,Nb
914
           READ(IOIN,'(A)',Iostat=Ios) Line
915
           if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE)
916
        END DO
917
     END IF
918
     READ(IOIN,'(A)',Iostat=Ios) Line
919
     if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE)
920
  END DO
921
  BACKSPACE(IOIN)
922

  
802 923
  ! We read the initial geometries
803 924
  Call Read_Geom(input)
804 925

  
......
809 930
  END DO
810 931

  
811 932
 If (AnaGeom) THEN
933
    If (NbVar>0) THEN
812 934
! We print the description of the varialbes in AnaFile
813
    Call PrintAnaList(IoDat)
935
       Call PrintAnaList(IoDat)
936
       If (OptGeom<=0) THEN
937
          WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
938
       ELSE
939
          WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
940
       END IF
941
    ELSE
942
       If (OptGeom<=0) THEN
943
          WRITE(IoDat,'(A)') "#  s    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"
944
       ELSE
945
          WRITE(IoDat,'(A)') "#  Iter    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"       
946
       END IF
947
    END IF
814 948
 END IF
815 949

  
816 950
  ! PFL 23 Sep 2008 ->
......
879 1013
              Dist=Norm1
880 1014
           END IF
881 1015
        END DO
882
        IAt2=Frozen(2)
883 1016
        ListAtFroz(Iat2)=.TRUE.
884 1017
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
885 1018
        x2(1)=x1(Iat2)-x1(Iat1)
......
888 1021
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
889 1022
        Dist=-1.
890 1023
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
1024
        Iat3=1
1025
        DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat))
1026
           Iat3=Iat3+1
1027
        END DO
891 1028
        DO I=1,Nat
892 1029
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
893 1030
           x2(2)=x1(Iat2)-x1(I)
......
909 1046
     DO I=2,NGeomI
910 1047
        ! First, we align the Ith geometry on the first one, using only
911 1048
        ! the positions of the "frozen" atoms. (see above: it might be that
912
        ! ListAtFroz contains non frozen atoms usd to align the geometries
1049
        ! ListAtFroz contains non frozen atoms used to align the geometries
913 1050
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
914 1051
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
915 1052
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
916 1053
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
917 1054

  
918

  
919 1055
        FChkFrozen=.FALSE.
920 1056
        DO J=1,NFroz
921 1057
           IAt=Frozen(J)
......
2101 2237
  Close(IOIN)
2102 2238
  Close(IOOUT)
2103 2239
  IF (AnaGeom) Close(IODAT)
2240
  IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot)
2104 2241

  
2105 2242
END PROGRAM PathOpt
src/egrad_ext.f90 (revision 8)
4 4
  ! a molecule, using an external code
5 5

  
6 6

  
7
  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, dzdc, indzmat,Nom,Atome, massat, unit,ProgExe
7
  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord,ProgExe
8 8
  use Io_module
9 9

  
10 10
  !
src/CalcCnct.f90 (revision 8)
31 31
SUBROUTINE CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)
32 32

  
33 33
  use Path_module, only : NMaxL, max_Z,Nom,Prog, KINT, KREAL
34
  
35
  IMPLICIT NONE
34 36

  
35 37
  integer(KINT) :: na, atome(Na)
36 38
  real(KREAL) ::  x(Na),y(Na),z(Na),fact
......
44 46

  
45 47
  ! Internals
46 48
  logical(KREAL) ::debug
49
  REAL(KREAL) :: DistTh
50
  INTEGER(KINT) :: I,j,iat
47 51

  
48 52
  INTERFACE
49 53
     function valid(string) result (isValid)
......
53 57
  END INTERFACE
54 58

  
55 59
  debug=valid('calccnct')
56
  if (debug) WRITE(*,*) "======================= Entering CalcCnct ======================="
60
  if (debug)  Call Header (" Entering CalcCnct ")
57 61

  
58 62
  DO i=1,na
59 63
     Liaisons(i,0)=0
......
113 117
        LIAISONS(i,0)=Nbli
114 118
     END DO
115 119
  END IF
116
  if (debug) WRITE(*,*) "======================= CalcCnct  over ======================="
120
  if (debug) Call Header (" CalcCnct over ")
117 121
END SUBROUTINE CalcCnct
src/Io_module.f90 (revision 8)
8 8
  SAVE
9 9

  
10 10
  INTEGER(KINT) :: IOIN=11, IOOUT=12, IOCART=14
11
  INTEGER(KINT) :: IOGEOM=15, IODAT=16
11
  INTEGER(KINT) :: IOGEOM=15, IODAT=16,IoGplot=17
12 12
  INTEGER(KINT), PARAMETER :: IOTMP=21,IOTMP2=22, IOTMP3=23
13 13
  INTEGER(KINT), PARAMETER :: IOERR=19
14 14
  CHARACTER(SCHARS) :: RunMode
......
19 19
!
20 20
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21 21

  
22
  REAL(KREAL), PARAMETER :: ev2au= 0.036749324445d0
22
  REAL(KREAL), PARAMETER :: au2eV=27.21183d0
23
  REAL(KREAL), PARAMETER :: ev2au= 1.d0/au2eV
24
  REAL(KREAL), PARAMETER :: au2kcal=627.509608d0
25
  REAL(KREAL), PARAMETER :: eV2kcal=23.06035d0
23 26

  
24 27
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
25 28
!
......
106 109
! OSuffix: suffix for the output file.
107 110
  CHARACTER(LCHARS) :: OSuffix
108 111

  
112
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
113
!
114
! For printing energies
115
!
116
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117
! Name of the unit used by internally by the 'engine' program
118
  CHARACTER(SCHARS) :: UnitProg="au"
119
! Conversion factor from energy program to kcal/mol
120
  REAL(KREAL) :: ConvE
109 121

  
110 122
END MODULE IO_MODULE
src/AnaPath.f90 (revision 8)
8 8

  
9 9

  
10 10
  use Io_module
11
  use Path_module, only : Nat,  NGeomF, XyzGeomF, AtName,SGeom,MassAt, &
12
       IReparam,NbVar, FormAna, Energies, au2kcal,XyzGeomI,Order,OrderInv, &
11
  use Path_module, only : Nat,  NGeomF, XyzGeomF, SGeom,MassAt, &
12
       IReparam,NbVar, FormAna, Energies,Order, &
13 13
       Renum, PrintGeomFactor, &
14
       Align, NFroz,  FrozAtoms
14
       Align, NFroz,  FrozAtoms, AnaFile,NbVar
15 15

  
16 16

  
17 17
  IMPLICIT NONE
18 18

  
19
! Input
20
! Iteration number
19
  ! Input
20
  ! Iteration number
21 21
  INTEGER(KINT), INTENT(IN) :: IOpt
22
! Unit file to print
22
  ! Unit file to print
23 23
  INTEGER(KINT), INTENT(IN) :: FileUnit
24 24

  
25 25

  
......
39 39
     END FUNCTION TEST_ZMAT
40 40

  
41 41

  
42
    SUBROUTINE die(routine, msg, file, line, unit)
42
     SUBROUTINE die(routine, msg, file, line, unit)
43 43

  
44
      Use VarTypes
45
      Use io_module
44
       Use VarTypes
45
       Use io_module
46 46

  
47
      implicit none
47
       implicit none
48 48

  
49
      character(len=*), intent(in)           :: routine, msg
50
      character(len=*), intent(in), optional :: file
51
      integer(KINT), intent(in), optional      :: line, unit
49
       character(len=*), intent(in)           :: routine, msg
50
       character(len=*), intent(in), optional :: file
51
       integer(KINT), intent(in), optional      :: line, unit
52 52

  
53
    END SUBROUTINE die
53
     END SUBROUTINE die
54 54

  
55
    SUBROUTINE Warning(routine, msg, file, line, unit)
55
     SUBROUTINE Warning(routine, msg, file, line, unit)
56 56

  
57
      Use VarTypes
58
      Use io_module
57
       Use VarTypes
58
       Use io_module
59 59

  
60
      implicit none
60
       implicit none
61 61

  
62
      character(len=*), intent(in)           :: routine, msg
63
      character(len=*), intent(in), optional :: file
64
      integer(KINT), intent(in), optional      :: line, unit
62
       character(len=*), intent(in)           :: routine, msg
63
       character(len=*), intent(in), optional :: file
64
       integer(KINT), intent(in), optional      :: line, unit
65 65

  
66
    END SUBROUTINE Warning
66
     END SUBROUTINE Warning
67 67

  
68 68
  END INTERFACE
69 69

  
......
74 74
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
75 75
  REAL(KREAL) ::  ds
76 76
  REAL(KREAL) :: MRot(3,3), Rmsd  
77
  CHARACTER(SCHARS) :: Plot="plot"
77 78

  
78 79
  debug=valid("AnaPath")
79 80

  
80 81
  if (debug) Call header("Entering AnaPath")
81 82

  
82
  if (debug) WRITE(*,*) "DBG AnaPaht PrintGeomFactor:",PrintGeomFactor
83
  !  if (debug) WRITE(*,*) "DBG AnaPaht PrintGeomFactor:",PrintGeomFactor
83 84
  ALLOCATE(GeomCart(Nat,3))
84
  ALLOCATE(Values(NbVar))
85
  If (NbVar>0) THEN
86
     ALLOCATE(Values(NbVar))
87
  END IF
85 88

  
86 89
  CalcDist=.FALSE.
87 90

  
......
94 97
  END DO
95 98

  
96 99

  
97
! Do we know the curvilinear values ?
100
  ! Do we know the curvilinear values ?
98 101
  IF (MOD(IOpt,IReparam)/=0) THEN
99 102
     CalcDist=.TRUE.
100 103
     SGeom=0.
......
116 119
     END DO
117 120
  END IF
118 121

  
119
  Call AnalyzeGeom(GeomCart,Values)
120
  WRITE(FileUnit,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
121
  if (debug) WRITE(*,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
122
  If (NbVar>0) THEN
123
     Call AnalyzeGeom(GeomCart,Values)
124
     WRITE(FileUnit,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
125
     if (debug) WRITE(*,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)
126
  ELSE
127
     WRITE(FileUnit,FormAna) SGeom(1),0.d0,Energies(1)
128
     if (debug) WRITE(*,FormAna) SGeom(1),0.d0,Energies(1)
129
  END IF
122 130

  
123 131
  DO I=2,NGeomF
124 132
     DO J=1,Nat
......
128 136
     END DO
129 137

  
130 138
     If (CalcDist) THEN
131
! PFL 2013 Feb 26
132
! For now, I do _NOT_ align the geometries
139
        ! PFL 2013 Feb 26
140
        ! For now, I do _NOT_ align the geometries
133 141
        if (debug) THEN
134 142
           WRITE(*,*) "AnaPath  - GeomCart avant align - I=",I
135 143
           DO K=1,Nat
......
138 146
        END IF
139 147

  
140 148
        IF ((Nat.GE.4).OR.Align) THEN
141
           WRITE(*,*) "DBG AnaPath: Aglin 2 geoms"
149
           if (debug) WRITE(*,*) "DBG AnaPath: Aglin 2 geoms"
142 150

  
143
! If we have frozen atoms we align only those ones.
151
           ! If we have frozen atoms we align only those ones.
144 152
           if (NFroz.GT.0) THEN
145 153
              Call AlignPartial(Nat,xRef,yRef,zRef,                 &
146 154
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),       &
......
150 158
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),   &
151 159
                   MRot,rmsd,.TRUE.,.TRUE.)
152 160
           END IF
153
              ! <- PFL 24 Nov 2008
154
           
161
           ! <- PFL 24 Nov 2008
162

  
155 163
        END IF
156 164

  
157 165
        if (debug) WRITE(*,*) "Mass:",MassAt(1:Nat)
158 166
        ds=0.
159 167
        DO J=1,Nat
160 168
           Iat=J
161
!        IF (renum) Iat=Order(J)
169
           !        IF (renum) Iat=Order(J)
162 170
           if (debug) WRITE(*,*) "Mass(J):",J,MassAt(Iat)
163
            
171

  
164 172
           do K=1,3
165 173
              ds=ds+MassAt(Iat)*(GeomCartPrec(J,K)-GeomCart(J,K))**2
166 174
           END DO
167 175
           if (debug) WRITE(*,*) "DBG DS:",J,ds
168 176
        END DO
169
        ds=sqrt(ds)
177
        ds=sqrt(ds) 
170 178
        if (debug) THEN
171 179
           WRITE(*,*) "AnaPath  - GeomCartPrec - I=",I
172 180
           DO K=1,Nat
......
180 188
        xRef=GeomCart(:,1)
181 189
        yRef=GeomCart(:,2)
182 190
        zRef=GeomCart(:,3)
183
        
191

  
184 192
     END IF
185 193

  
186 194
     if (debug) THEN
......
190 198
        END DO
191 199
     END IF
192 200

  
193
     Call AnalyzeGeom(GeomCart,Values)
194
     WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*au2kcal,Energies(i)
201
     If (NbVar>0) THEN
202
        Call AnalyzeGeom(GeomCart,Values)
203
        WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*ConvE,Energies(i)
204
        IF (Debug) WRITE(*,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*ConvE,Energies(i)
205
     ELSE
206
        WRITE(FileUnit,FormAna) SGeom(i),(Energies(i)-Energies(1))*ConvE,Energies(i)
207
        If (debug) WRITE(*,FormAna) SGeom(i),(Energies(i)-Energies(1))*ConvE,Energies(i)
208
     END IF
195 209
  END DO
196 210

  
197 211
  WRITE(FileUnit,*) ""
198 212
  WRITE(FileUnit,*) ""
199 213

  
200
  DeAllocate(GeomCart,Values)
214
  DeAllocate(GeomCart)
215
  If (NbVar>0) DeAllocate(Values)
201 216
  If (CalcDist) Deallocate(GeomCartPrec,xRef,yRef,zRef)
202 217

  
218
  Write(IoGplot,'(A,I5,A)') TRIM(Plot) // ' "' // TRIM(AnaFile) // '"  i ',Iopt," u xcol:Ecol w lp "
219
  Plot="replot"
220

  
203 221
  if (debug) Call header("AnaPath over")
204 222

  
205 223
END SUBROUTINE AnaPath
src/egrad_gaussian.f90 (revision 8)
3 3
  ! This routines calculates the energy and the gradient of 
4 4
  ! a molecule, using Gaussian 
5 5

  
6
  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, dzdc, indzmat,Nom,Atome, massat, unit,ProgExe
6
  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord,unit,ProgExe
7 7
  use Io_module
8 8

  
9 9
  IMPLICIT NONE
src/Read_geom.f90 (revision 8)
8 8

  
9 9

  
10 10
  CHARACTER(32), INTENT(IN) :: input    
11
  CHARACTER(132) :: LineTmp,Line
12 11

  
13
  INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na
14
  INTEGER(KINT) :: NbType, NbTypeUser
15
  INTEGER(KINT) :: I, J, Iat,NAtP,Idx,JStart
16
  INTEGER(KINT) :: IGeom
17
  REAL(KREAL) :: Xt,Yt,Zt
18
  LOGICAL :: Debug, TChk
19
  CHARACTER(4), ALLOCATABLE :: FFFt(:,:) ! 3,Na
20
  LOGICAL, ALLOCATABLE :: FFFl(:,:) ! 3,Na
12
  LOGICAL :: Debug
21 13

  
22
!!! for VASP inputs
23
! latx_loc used to check that the lattice parameters are the same for all geometries
24
  REAL(KREAL) :: lata_loc(3), latb_loc(3), latc_loc(3)
25
! VaspPar_loc is the local value of the vasp parametre
26
  REAL(KREAL) :: VaspPar_loc
27
! Vdir_loc is used to convert the read geometry to Cartesian.
28
! V_direct is set by the first geometry.
29
  CHARACTER(LCHARS) :: Vdir_loc
30 14

  
31

  
32 15
  INTERFACE
33 16
     function valid(string) result (isValid)
34 17
       CHARACTER(*), intent(in) :: string
src/Calc_zmat_constr_frag.f90 (revision 8)
9 9
  !     we first construct the zmat using only the frozen atoms, and then
10 10
  !     we add the other atoms on top of the first ones...
11 11
  !     
12
  Use Path_module, only : max_Z, NMaxL, Nom,MaxFroz
12
  Use Path_module, only : max_Z, NMaxL, Nom
13 13
  Use Io_module
14 14

  
15 15
  IMPLICIT NONE
src/ReadAnaList.f90 (revision 8)
31 31

  
32 32
  END INTERFACE
33 33

  
34
  LOGICAL :: Debug,FirstVar,FirstCom
34
  LOGICAL :: Debug
35 35
  INTEGER(KINT) :: At1,At2,At3,At4
36 36
  INTEGER(KINT) :: I,Idx,J
37 37

  
......
50 50

  
51 51
  NbCom=0
52 52
  NbVar=0
53
  FormAna='(1X'
53
  FormAna='(1X,F8.3'
54 54
  DO I=1, Nb
55 55
     READ(IOIN,'(A)') Line
56 56
     Line=AdjustL(Line)
......
110 110
           if (debug) THEN
111 111
              WRITE(*,'("# d ",4(I3))') At1,At2,At3,At4
112 112
           END IF
113
           FormAna=TRIM(FormAna) //',1X,F7.2'
113 114
           Allocate(CurVar%Next)
114 115
           CurVar => CurVar%Next
115 116
           Nullify(CurVar%Next)
......
129 130
              WRITE(*,'("# c ",I4,20(I3))') At1, &
130 131
                   (CurBary%ListAtoms(j),j=1,At1)
131 132
           END IF
132
           FormAna=Trim(FormAna) //',1X,F7.2'
133

  
134 133
           Allocate(CurBary%Next)
135 134
           CurBary => CurBary%Next
136 135
           Nullify(CurBary%Next)
src/Calc_mixed_frag.f90 (revision 8)
9 9
  !     we first construct the zmat using only the frozen atoms, and then
10 10
  !     we add the other atoms on top of the first ones...
11 11
  !     
12
  Use Path_module, only : max_Z, NMaxL, Nom,MaxFroz
12
  Use Path_module, only : max_Z, NMaxL, Nom
13 13
  Use Io_module
14 14

  
15 15
  IMPLICIT NONE
......
1011 1011

  
1012 1012
! This subroutine comes for zmat_g92
1013 1013

  
1014
  Use Path_module, only : max_Z, NMaxL, Nom,MaxFroz, Pi
1014
  Use Path_module, only :   Nom
1015 1015
  Use Io_module
1016 1016

  
1017 1017
  IMPLICIT NONE
src/Calc_zmat.f90 (revision 8)
2 2
           r_cov,  fact)
3 3

  
4 4
      use Path_module, only : Pi,max_Z, NMaxL, Nom, KINT, KREAL
5
      use Io_module, only : IoIn, IoOut
5
      use Io_module, only :  IoOut
6 6

  
7 7
      IMPLICIT NONE
8 8

  
src/EgradPath.f90 (revision 8)
6 6
  use Path_module, only : ProgExe, NCoord, PROG, Nat, NGeomF, Coord, IntCoordF, &
7 7
       IndZmat, XyzGeomF, Atome, Order, OrderInv, Latr, FrozAtoms, &
8 8
       Grad, Energies, renum, dzdc, massat, Pi, NCART, OptReac, &
9
       OptProd, IntCoordI, GeomOld_all, AtName, Unit, CalcEReac, CalcEProd
9
       OptProd,  AtName, Unit, CalcEReac, CalcEProd
10 10
  ! GeomOld_all(NGeomF,NCoord), BTransInv (3*Nat,NCoord): used for Baker case
11 11
  use Io_module
12 12

  
src/egrad_LEPS.f90 (revision 8)
14 14
       REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3)
15 15

  
16 16

  
17
!     Bohr --> Angstr
18
      real(KREAL), parameter :: BOHR   = 0.52917726D+00
19
!
20 17
! Parameters to define the surface
21 18
      INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/-1,9,-45,45,-9,1/)
22 19
      INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/-3,-2,-1,1,2,3/)
23 20
      REAL(KREAL), PARAMETER :: hh=0.001d0
24 21
      
25 22
! Variables
26
       INTEGER(KINT) :: i,j,iat,jat
23
       INTEGER(KINT) :: i,iat,jat
27 24
       REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)
28
       REAL(KREAL) :: xp,yp
29

  
30
       CHARACTER(132) :: Line
25
 
31 26
       LOGICAL :: Debug
32 27
  
33
      real(KREAL), parameter :: zero = 0.0_KREAL,   one = 1.0_KREAL
34 28
      REAL(KREAL), external :: ELEPS_xyz
35 29

  
36 30

  
......
94 88
      function ELEPS_xyz(natoms,Xyz)
95 89

  
96 90
        use Path_module, only : order
91
        use Io_module, only : au2ev
97 92

  
98 93
   IMPLICIT NONE
99 94
     integer, parameter :: KINT = kind(1)
......
113 108

  
114 109
      REAL(KREAL) :: Qbc,QAc,Qab,Jab,Jbc,Jac
115 110

  
116
      Real(KREAL), Parameter ::  autoA=0.52917715d0,autoeV=27.21183d0
117

  
118 111
      rAB=0.
119 112
      rAC=0.
120 113
      rBC=0.
......
140 133
! V(x,y) = Qab(x)+Qbc(y)+Qac(rac)-sqrt(Jab(x)**2+Jbc(y)**2+Jac(x)**2-Jab(x)*Jbc(y)-Jbc(y)*Jac(rac)-Jab(x)*Jac(rac))
141 134

  
142 135

  
143
      ELEPS_xyz=(Qab+Qbc+Qac-sqrt(Jab**2+Jbc**2+Jac**2-Jab*Jbc-Jbc*Jac-Jab*Jac))/autoeV
136
      ELEPS_xyz=(Qab+Qbc+Qac-sqrt(Jab**2+Jbc**2+Jac**2-Jab*Jbc-Jbc*Jac-Jab*Jac))/au2eV
144 137

  
145 138
      return
146 139
    end function ELEPS_xyz
src/egrad_test_2D.f90 (revision 8)
126 126

  
127 127
      function E2D(natoms,Xyz)
128 128

  
129
        use Path_module, only : order
130 129
        use Io_module, only : IOTMP
131 130

  
132 131
   IMPLICIT NONE
......
140 139
      REAL(KREAL), save :: Data(501,501)
141 140

  
142 141
      INTEGER(KINT) :: i,j
143
      REAL(KREAL) :: dCH,dNH,dCN
144 142

  
145
!      real(KREAL) :: autoA,autoeV
146
      real(KREAL) :: autoeV
147
!      parameter (autoA=0.52917715d0)
148
      parameter (autoeV=27.21183d0)
149

  
150

  
151
      real(KREAL) :: r1, r2
152

  
153 143
      LOGICAL :: Debug
154 144
      LOGICAL, SAVE :: First=.TRUE.
155 145
 
src/egrad_chamfre.f90 (revision 8)
11 11
! WAS xa,ya,za,v,dvdr)
12 12

  
13 13

  
14
  use Path_module, only : Lat_a,Lat_b,Lat_c,order,orderinv,renum
14
  use Path_module, only : Lat_a,Lat_b,order,renum
15 15
  use Io_module
16 16

  
17 17

  
src/PrintAnaList.f90 (revision 8)
40 40

  
41 41
  If (Debug) Call Header("Entering PrintAnaList")
42 42

  
43
  WRITE(FileUnit,'(A)') "# Centers of mass used "
44

  
43 45
  CurBary => Bary
44 46

  
45 47
  DO I=1,NbCom
......
53 55
     if (associated(CurBary%Next)) CurBary=> CurBary%Next
54 56
  END DO
55 57

  
58
  WRITE(FileUnit,'(A)') "# Variables monitored "
59

  
56 60
  CurVar => GeomList
57 61
  DO I=1,NbVar
58 62
     SELECT CASE (CurVar%Type)
src/Opt_Geom.f90 (revision 8)
3 3
!so  It has been designed to be almost independent of the rest of the code.
4 4
! It is now an (almost) officiel option.
5 5

  
6
 SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method)
6
SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method)
7 7

  
8 8
  use VarTypes
9 9
  use Path_module, only : maxcyc,sthresh,gthresh,Hinv,Smax,Nom,Atome,OrderInv,indzmat,&
10
                                            XyzGeomF,IntCoordF,UMat,ScanCoord,Coordinate,NPrim,&
11
                                            BTransInv,BTransInv_local,XyzGeomI,Xprimitive_t,ScanCoord,&
12
                                            Coordinate,Pi,UMat_local,NCart, DynMaxStep,FPrintGeom,GeomFile, AtName,Renum,Order
13
  use Io_module, only : IoGeom
10
       XyzGeomF,IntCoordF,UMat,ScanCoord,Coordinate,NPrim,&
11
       BTransInv,BTransInv_local,XyzGeomI,Xprimitive_t,ScanCoord,&
12
       Coordinate,Pi,UMat_local,NCart, DynMaxStep,FPrintGeom,GeomFile, AtName,Renum,Order, &
13
       FormAna,NbVar, PrintGeomFactor,AnaGeom
14 14

  
15
  use Io_module
16

  
15 17
  IMPLICIT NONE
16 18

  
17 19
  INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom
......
30 32

  
31 33
     subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old)
32 34

  
33
  ! This routines calculates the energy E and the gradient Grad of 
34
  ! a molecule with Geometry Geom (may be in internal coordinates),
35
  ! using for now, either Gaussian or Ext, more general later.
35
       ! This routines calculates the energy E and the gradient Grad of 
36
       ! a molecule with Geometry Geom (may be in internal coordinates),
37
       ! using for now, either Gaussian or Ext, more general later.
36 38

  
37
  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
38
       prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &
39
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
40
       , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
41
      , Order,OrderInv, XPrimitiveF
39
       use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
40
            prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &
41
            GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
42
            , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
43
            , Order,OrderInv, XPrimitiveF
42 44

  
43
  ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
44
  ! allocated in Path.f90
45
       ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
46
       ! allocated in Path.f90
45 47

  
46
  use Io_module
48
       use Io_module
47 49

  
48
  ! Energy (calculated if F300K=.F., else estimated)
49
  REAL(KREAL), INTENT (OUT) :: E
50
  ! NCoord: Number of the degrees of freedom
51
  ! IGeom: index of the geometry.
52
  INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
53
  ! Geometry at which gradient is calculated (cf Factual also):
54
  REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
55
  ! Gradient calculated at Geom geometry:
56
  REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
57
  ! Cartesian geometry corresponding to (Internal Geometry) Geom:
58
  REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)
50
       ! Energy (calculated if F300K=.F., else estimated)
51
       REAL(KREAL), INTENT (OUT) :: E
52
       ! NCoord: Number of the degrees of freedom
53
       ! IGeom: index of the geometry.
54
       INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
55
       ! Geometry at which gradient is calculated (cf Factual also):
56
       REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
57
       ! Gradient calculated at Geom geometry:
58
       REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
59
       ! Cartesian geometry corresponding to (Internal Geometry) Geom:
60
       REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)
59 61
!!! Optional, just for geometry optimization with Baker coordinates
60
  REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
61
  REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
62
! FOptGeom is a flag indicating if we are doing a geom optimization
63
! it can be omitted so that we use a local flag: Flag_Opt_Geom
64
  LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
65
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
66
  LOGICAL  :: Flag_Opt_Geom
62
       REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
63
       REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
64
       ! FOptGeom is a flag indicating if we are doing a geom optimization
65
       ! it can be omitted so that we use a local flag: Flag_Opt_Geom
66
       LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
67
       ! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
68
       LOGICAL  :: Flag_Opt_Geom
67 69

  
68
 END subroutine Egrad
70
     END subroutine Egrad
69 71

  
70 72
  END INTERFACE
71 73

  
......
74 76
  LOGICAL :: Fini
75 77
  LOGICAL, SAVE :: FRST=.TRUE.
76 78

  
77
! Variables
79
  ! Variables
78 80
  INTEGER(KINT) :: IOpt, I,J,Iat, IBEG
79 81
  REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6
80 82
  REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)
......
89 91
  REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)
90 92
  REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis
91 93
  REAL(KREAL) :: NormStep, FactStep, HP, MaxStep
92
  REAL(KREAL) :: Eold  
94
  REAL(KREAL) :: Eold, Eini
95
  ! Values contains the values for the geometry analizes
96
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
93 97

  
94 98
  debug=valid('OptGeom')
95 99

  
96
    E=0.
97
    Eold=0.
98
    MaxStep=SMax
100
  E=0.
101
  Eold=0.
102
  Eini=0.
103
  MaxStep=SMax
99 104

  
100
    if (debug) Call Header("Entering Geom Optimization")
105
  if (debug) Call Header("Entering Geom Optimization")
101 106

  
102
     ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
103
     ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
104
     ALLOCATE(GeomRef(NCoord))
105
     ALLOCATE(Hess_local(NCoord,NCoord))
106
     ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
107
     ALLOCATE(NewGeom(NCoord))
108
     ALLOCATE(NewGradTmp(NCoord))
109
     ALLOCATE(Hess_local_inv(NCoord,NCoord))
110
        
111
     IOpt=0
112
     ZeroVec=0.d0
113
     
114
     ! Initialize the Hessian
115
     Hess_local=0.
107
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
108
  ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
109
  ALLOCATE(GeomRef(NCoord))
110
  ALLOCATE(Hess_local(NCoord,NCoord))
111
  ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))
112
  ALLOCATE(NewGeom(NCoord))
113
  ALLOCATE(NewGradTmp(NCoord))
114
  ALLOCATE(Hess_local_inv(NCoord,NCoord))
116 115

  
117
     SELECT CASE (COORD)
118
        CASE ('ZMAT')
119
        ! We use the fact that we know that approximate good hessian values
120
        ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
121
        Hess_local(1,1)=0.5d0
122
        Hess_local(2,2)=0.5d0
123
        Hess_local(3,3)=0.1d0
124
        DO J=1,Nat-3
125
           Hess_local(3*J+1,3*J+1)=0.5d0
126
           Hess_local(3*J+2,3*J+2)=0.1d0
127
           Hess_local(3*J+3,3*J+3)=0.02d0
116
  if (NbVar>0) THEN
117
     ALLOCATE(Values(NbVar))
118
  END IF
119
  FormAna(5:8)=" I5 "
120
  IOpt=0
121
  ZeroVec=0.d0
122

  
123
  ! Initialize the Hessian
124
  Hess_local=0.
125

  
126
  SELECT CASE (COORD)
127
 CASE ('ZMAT')
128
     ! We use the fact that we know that approximate good hessian values
129
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
130
     Hess_local(1,1)=0.5d0
131
     Hess_local(2,2)=0.5d0
132
     Hess_local(3,3)=0.1d0
133
     DO J=1,Nat-3
134
        Hess_local(3*J+1,3*J+1)=0.5d0
135
        Hess_local(3*J+2,3*J+2)=0.1d0
136
        Hess_local(3*J+3,3*J+3)=0.02d0
137
     END DO
138
     IF (HInv) THEN
139
        DO I=1,NCoord
140
           Hess_local(I,I)=1.d0/Hess_local(I,I)
128 141
        END DO
129
        IF (HInv) THEN
130
           DO I=1,NCoord
131
              Hess_local(I,I)=1.d0/Hess_local(I,I)
132
           END DO
133
        END IF
142
     END IF
134 143

  
135
        IF (Step_method .EQ. "RFO") Then
136
                   Hess_local=0.
137
           DO I=1,NCoord
138
              Hess_local(I,I)=0.5d0
139
           END DO
140
        END IF
144
     IF (Step_method .EQ. "RFO") Then
145
        Hess_local=0.
146
        DO I=1,NCoord
147
           Hess_local(I,I)=0.5d0
148
        END DO
149
     END IF
141 150

  
142
     CASE ('BAKER')
143
            ! UMat(NPrim,3*Nat-6)
144
        BTransInv_local = BTransInv
145
        UMat_local = UMat
146
         ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord))
147
         Hprim=0.d0
148
         ScanCoord=>Coordinate
149
         I=0
150
         DO WHILE (Associated(ScanCoord%next))
151
            I=I+1
152
            SELECT CASE (ScanCoord%Type)
153
                     CASE ('BOND')
154
                          Hprim(I,I) = 0.5d0
155
            CASE ('ANGLE')
156
                          Hprim(I,I) = 0.2d0
157
            CASE ('DIHEDRAL')
158
                          Hprim(I,I) = 0.1d0
159
               END SELECT
160
            ScanCoord => ScanCoord%next
161
         END DO
162
               
163
         ! Hprim U:
164
         HprimU=0.d0
165
         DO I=1, NCoord
166
                  DO J=1,NPrim
167
               HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
168
            END DO
169
         END DO
151
  CASE ('BAKER')
152
     ! UMat(NPrim,3*Nat-6)
153
     BTransInv_local = BTransInv
154
     UMat_local = UMat
155
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord))
156
     Hprim=0.d0
157
     ScanCoord=>Coordinate
158
     I=0
159
     DO WHILE (Associated(ScanCoord%next))
160
        I=I+1
161
        SELECT CASE (ScanCoord%Type)
162
        CASE ('BOND')
163
           Hprim(I,I) = 0.5d0
164
        CASE ('ANGLE')
165
           Hprim(I,I) = 0.2d0
166
        CASE ('DIHEDRAL')
167
           Hprim(I,I) = 0.1d0
168
        END SELECT
169
        ScanCoord => ScanCoord%next
170
     END DO
170 171

  
171
         ! Hess = U^T Hprim U:
172
         Hess_local=0.d0
173
         DO I=1, NCoord
174
            DO J=1,NPrim
175
                     ! UMat^T is needed below.
176
                        Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I)
177
            END DO
178
         END DO
179
               
180
               !Print *, 'Hprim='
181
               DO I=1,NPrim
182
               !   WRITE(*,'(10(1X,F6.3))') Hprim(I,:)
183
            END DO
184
               !Print *, 'UMat='
185
               DO I=1,NPrim
186
                !  WRITE(*,'(3(1X,F6.3))') UMat(I,1:3)
187
            END DO
188
               !Print *, 'HprimU='
189
               DO I=1,NPrim
190
                !  WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3)
191
            END DO
192
               !Print *, 'Hess_local='
193
               DO I=1,NCoord
194
                !  WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
195
            END DO
196
               
197
            DEALLOCATE(Hprim,HprimU)
198
               
199
         IF (HInv) THEN
200
                     Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord)
201
                  Hess_local = Hess_local_inv
202
        END IF
203
              
204
              !Print *, 'Hess_local after inversion='
205
               DO I=1,NCoord
206
               !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
207
            END DO
208
                             
209
        IF (Step_method .EQ. "RFO") Then
210
                 Hess_local=0.
211
           DO I=1,NCoord
212
              Hess_local(I,I)=0.5d0
213
           END DO
214
        END IF
215
                 
216
        CASE ('MIXED','CART','HYBRID')
217
        DO J=1,NCoord
218
           Hess_local(J,J)=1.
172
     ! Hprim U:
173
     HprimU=0.d0
174
     DO I=1, NCoord
175
        DO J=1,NPrim
176
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
219 177
        END DO
220
        CASE DEFAULT
221
        WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"
222
        STOP
223
     END SELECT
178
     END DO
224 179

  
225
     ! Go to optimization
226
     GeomOld=0.d0 ! Internal coordinates
227
     GeomCart_old=0.d0 ! Cartesian coordinates
180
     ! Hess = U^T Hprim U:
181
     Hess_local=0.d0
182
     DO I=1, NCoord
183
        DO J=1,NPrim
184
           ! UMat^T is needed below.
185
           Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I)
186
        END DO
187
     END DO
228 188

  
229
     IF (FPrintGeom) THEN
230
        OPEN(IOGeom,File=GeomFile)
189
     !Print *, 'Hprim='
190
     DO I=1,NPrim
191
        !   WRITE(*,'(10(1X,F6.3))') Hprim(I,:)
192
     END DO
193
     !Print *, 'UMat='
194
     DO I=1,NPrim
195
        !  WRITE(*,'(3(1X,F6.3))') UMat(I,1:3)
196
     END DO
197
     !Print *, 'HprimU='
198
     DO I=1,NPrim
199
        !  WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3)
200
     END DO
201
     !Print *, 'Hess_local='
202
     DO I=1,NCoord
203
        !  WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
204
     END DO
205

  
206
     DEALLOCATE(Hprim,HprimU)
207

  
208
     IF (HInv) THEN
209
        Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord)
210
        Hess_local = Hess_local_inv
231 211
     END IF
232 212

  
233
     Fini=.FALSE.
234
     DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
235
        IOpt=IOpt+1
213
     !Print *, 'Hess_local after inversion='
214
     DO I=1,NCoord
215
        !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
216
     END DO
236 217

  
237
        ! Calculate the  energy and gradient
238
        IF (debug) THEN 
239
           WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
240
           WRITE(*,*) "Energy and Coordinates:"
241
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
242
           WRITE(*,*) "Geom Old:"
243
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
244
           WRITE(*,*) "Grad Old:"
245
           WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))
246
        END IF
218
     IF (Step_method .EQ. "RFO") Then
219
        Hess_local=0.
220
        DO I=1,NCoord
221
           Hess_local(I,I)=0.5d0
222
        END DO
223
     END IF
247 224

  
248
        !Call EGrad(E,Geom,GradTmp,NCoord)
249
        IF (COORD.EQ.'BAKER') THEN
250
           ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
251
           ! GeomCart has INTENT(OUT)
252
           ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable
253
           GeomRef=GeomOld
254
           Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
255
                GeomRef,GeomCart_old)
256
           GeomCart_old=GeomCart
225
  CASE ('MIXED','CART','HYBRID')
226
     DO J=1,NCoord
227
        Hess_local(J,J)=1.
228
     END DO
229
  CASE DEFAULT
230
     WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"
231
     STOP
232
  END SELECT
233

  
234
  ! Go to optimization
235
  GeomOld=0.d0 ! Internal coordinates
236
  GeomCart_old=0.d0 ! Cartesian coordinates
237

  
238
  IF (FPrintGeom) THEN
239
     OPEN(IOGeom,File=GeomFile)
240
  END IF
241

  
242
  Fini=.FALSE.
243
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
244
     IOpt=IOpt+1
245

  
246
     ! Calculate the  energy and gradient
247
     IF (debug) THEN 
248
        WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
249
        WRITE(*,*) "Energy and Coordinates:"
250
        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
251
        WRITE(*,*) "Geom Old:"
252
        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))
253
        WRITE(*,*) "Grad Old:"
254
        WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))
255
     END IF
256

  
257
     !Call EGrad(E,Geom,GradTmp,NCoord)
258
     IF (COORD.EQ.'BAKER') THEN
259
        ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
260
        ! GeomCart has INTENT(OUT)
261
        ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable
262
        GeomRef=GeomOld
263
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
264
             GeomRef,GeomCart_old)
265
        GeomCart_old=GeomCart
266
     ELSE
267
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
268
     END IF
269

  
270
     If (Iopt==1) EIni=E
271

  
272
     IF (debug) THEN 
273
        WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
274
        WRITE(*,*) "Energy and Coordinates:"
275
        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
276
        WRITE(*,*) "Geom Old:"
277
        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,12))
278
        WRITE(*,*) "Grad:"
279
        WRITE(*,'(F12.6,12(1X,F8.3))') GradTmp(1:min(NCoord,12))
280
     END IF
281

  
282
     IF (FPrintGeom) THEN
283
        WRITE(IoGeom,*) Nat
284
        WRITE(IoGeom,"(' Geom. ',I5,' Energy= ',F15.6)") IOpt, E
285

  
286
        DO I=1,Nat
287
           If (renum) THEN
288
              Iat=Order(I)
289
              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat-2:3*Iat)
290
           ELSE
291
              Iat=OrderInv(I)
292
              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I-2:3*I)
293
           END IF
294
        END DO
295
     END IF
296

  
297
     ! do we have to analyze geometries ?
298
     If (AnaGeom) THEN
299
        If (NbVar>0) THEN
300
           Call AnalyzeGeom(GeomCart,Values)
301
           WRITE(IoDat,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E
302
           if (debug) WRITE(*,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E
257 303
        ELSE
258
           Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
304
           WRITE(IoDat,FormAna) Iopt,(E-Eini)*ConvE,E
305
           if (debug) WRITE(*,FormAna) Iopt,(E-Eini)*ConvE,E
259 306
        END IF
260
              
261
        IF (debug) THEN 
262
           WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
263
           WRITE(*,*) "Energy and Coordinates:"
264
           WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))
265
           WRITE(*,*) "Geom Old:"
266
           WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,12))
267
           WRITE(*,*) "Grad:"
268
           WRITE(*,'(F12.6,12(1X,F8.3))') GradTmp(1:min(NCoord,12))
269
        END IF
307
     END IF
270 308

  
271
        IF (FPrintGeom) THEN
272
           WRITE(IoGeom,*) Nat
273
           WRITE(IoGeom,"(' Geom. ',I5,' Energy= ',F15.6)") IOpt, E
274 309

  
275
           DO I=1,Nat
276
              If (renum) THEN
277
                 Iat=Order(I)
278
                 WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat-2:3*Iat)
279
              ELSE
280
                 Iat=OrderInv(I)
281
                 WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I-2:3*I)
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff