Révision 7

src/ReadInput_siesta.F90 (revision 7)
1
 SUBROUTINE ReadInput_siesta
2

  
3
! This routine reads an input template for Siesta
4

  
5
  use VarTypes
6
  use Path_module
7
  use Io_module
8

  
9
  IMPLICIT NONE
10

  
11
  INTERFACE
12
     function valid(string) result (isValid)
13
       CHARACTER(*), intent(in) :: string
14
       logical                  :: isValid
15
     END function VALID
16
  END INTERFACE
17

  
18

  
19
  CHARACTER(132) ::  Line,LineUp,Line2
20
  INTEGER(KINT) :: LineL, Idx, Iat
21
  INTEGER(KINT) :: IoRead
22
  
23
  LOGICAL :: Debug
24
  LOGICAL :: FSpecies, FCoord
25

  
26

  
27
  Debug=Valid("readinput").OR.Valid("readinput_siesta")
28

  
29
 if (debug) Call Header("Entering ReadInput_Siesta")
30

  
31
  ! We read the Siesta input file
32
     ALLOCATE(Siesta_Input)
33
     NULLIFY(Siesta_Input%next)
34
     Current => Siesta_Input
35

  
36
     FSpecies=.FALSE.
37
     FCoord=.FALSE.
38

  
39
     READ(IOIN,'(A)',iostat=IoRead) Line
40
     DO WHILE (IoRead==0)
41
        Line=AdjustL(Line)
42
        LineL=len(Trim(Line))
43
        LineUp=Line
44
        Call UpCase(LineUp)
45

  
46
! We look for the NumberofAtoms
47
        Idx=INDEX(LineUp,"NUMBEROFATOMS")
48
        IF (Idx/=0) THEN
49
           Idx=Index(Line," ")
50
           Line2=Trim(AdjustL(Line(Idx+1:)))
51
           Read(Line2,*) IAt
52
           if (Iat/=Nat) THEN
53
              Call Die('ReadInput_siesta','Nat in FDF sample different from  Nat Path input', & 
54
                   __FILE__, __LINE__, IOOUT)
55
                   
56
           END IF
57
        END IF
58

  
59
! We look for the ChemicalSpeciesLabel block
60
        Idx=INDEX(LineUp,"CHEMICALSPECIESLABEL")
61
        IF (Idx/=0) THEN
62
           FSpecies=.TRUE.
63
        END IF
64

  
65
! We look for SystemLabel
66
        Idx=INDEX(LineUp,"SYSTEMLABEL")
67
        If (Idx/=0) THEN
68
           Idx=Index(Line," ")
69
           Siesta_Label=Trim(adjustl(Line(Idx+1:)))
70
        END IF
71

  
72
! We look for Coordinates
73

  
74
        IF ((LineL.NE.0).AND.(Idx.EQ.0)) THEN
75
           current%Line=TRIM(Line)
76
           ALLOCATE(current%next)
77
           Current => Current%next
78
           Nullify(Current%next)
79
        END IF
80
     END DO
81

  
82
     if (debug) Call Header("Exiting ReadInput_Siesta")
83

  
84
   END SUBROUTINE READINPUT_SIESTA
src/sinangle.f90 (revision 7)
1 1
     FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
2 2

  
3
       use Path_module, only :  Pi,KINT, KREAL
3
       use Path_module, only :  KREAL
4 4

  
5 5
       IMPLICIT NONE
6 6

  
src/PrintAnaList.f90 (revision 7)
1
 SUBROUTINE PrintAnaList(FileUnit)
2
! This routines prints the list of geometrical variables to monitor
3

  
4

  
5
  use VarTypes
6
  use Path_module
7
  use Io_module
8

  
9
  IMPLICIT NONE
10

  
11
  INTERFACE
12
     function valid(string) result (isValid)
13
       CHARACTER(*), intent(in) :: string
14
       logical                  :: isValid
15
     END function VALID
16

  
17
     SUBROUTINE die(routine, msg, file, line, unit)
18

  
19
       Use VarTypes
20
       Use io_module
21

  
22
       implicit none
23

  
24
       character(len=*), intent(in)           :: routine, msg
25
       character(len=*), intent(in), optional :: file
26
       integer(KINT), intent(in), optional      :: line, unit
27

  
28
     END SUBROUTINE die
29

  
30
  END INTERFACE
31
! Input
32
! Unit to print the description
33
  INTEGER(KINT) :: FileUnit
34

  
35
! Local
36
  LOGICAL :: Debug
37
  INTEGER(KINT) :: I,J,At1,At2,At3,At4
38

  
39
  Debug=Valid('PrintAnaList')
40

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

  
43
  CurBary => Bary
44

  
45
  DO I=1,NbCom
46
     WRITe(*,*) "Dbg PrintAnalist: I,NbCom",I,NbCom
47
     At1=CurBary%ListAtoms(0)
48
     if (debug)      WRITE(*,'("# c ",I4,20(A3,I3))') At1, &
49
                (AtName(CurBary%ListAtoms(j)),CurBary%ListAtoms(j),j=1,At1)
50

  
51
     WRITE(FileUnit,'("# c ",I4,20(A3,I3))') At1, &
52
                (AtName(CurBary%ListAtoms(j)),CurBary%ListAtoms(j),j=1,At1)
53
     if (associated(CurBary%Next)) CurBary=> CurBary%Next
54
  END DO
55

  
56
  CurVar => GeomList
57
  DO I=1,NbVar
58
     SELECT CASE (CurVar%Type)
59
        CASE ('BOND') 
60
! this is a bond
61
           At1=CurVar%At1
62
           AT2=CurVar%At2
63
           if (debug) WRITE(*,'("# b ",2(1X,A,I3))') TRIM(AtName(At1)),At1,AtName(At2),At2
64
           WRITE(FileUnit,'("# b ",2(1X,A,I3))') TRIM(AtName(At1)),At1,AtName(At2),At2
65
        CASE ('ANGLE') 
66
! this is an angle
67
           At1=CurVar%At1
68
           AT2=CurVar%At2
69
           At3=CurVar%At3
70
           if (debug) WRITE(*,'("# a ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3
71
           WRITE(FileUnit,'("# a ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3
72
        CASE ('DIHEDRAL') 
73
! this is a dihedral
74
           At1=CurVar%At1
75
           AT2=CurVar%At2
76
           At3=CurVar%At3
77
           At4=CurVar%At4
78
           if (debug) WRITE(*,'("# d ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3,AtName(At4),At4
79
           WRITE(FileUnit,'("# d ",4(A3,I3))') AtName(At1),At1,AtName(At2),At2,AtName(At3),At3,AtName(At4),At4
80
     END SELECT
81
     If (Associated(CurVar%Next)) CurVar => CurVar%next
82
  END DO
83

  
84
  If (Debug) Call Header("Exiting PrintAnaList")
85

  
86
END SUBROUTINE PrintAnaList
src/AnalyzeGeom.f90 (revision 7)
1
 SUBROUTINE AnalyzeGeom(GeomCart,Values)
2
! This routines read a list of geometrical variables to monitor
3
! This is inspired from Xyz2Path (that was inspired by Xyz2scan ...)
4

  
5

  
6
  use VarTypes
7
  use Path_module
8
  use Io_module
9

  
10
  IMPLICIT NONE
11

  
12

  
13
  INTERFACE
14
     function valid(string) result (isValid)
15
       CHARACTER(*), intent(in) :: string
16
       logical                  :: isValid
17
     END function VALID
18

  
19
     SUBROUTINE die(routine, msg, file, line, unit)
20

  
21
       Use VarTypes
22
       Use io_module
23

  
24
       implicit none
25

  
26
       character(len=*), intent(in)           :: routine, msg
27
       character(len=*), intent(in), optional :: file
28
       integer(KINT), intent(in), optional      :: line, unit
29

  
30
     END SUBROUTINE die
31

  
32

  
33
     SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
34

  
35
       Use VarTypes
36
       Use Io_module
37
       Use Path_module, only : pi
38

  
39
       IMPLICIT NONE
40

  
41
       Type (ListCoord), POINTER :: Coordinate
42
       INTEGER(KINT), INTENT(IN) :: Nat,NPrim
43
       REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat)
44
       REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
45
       REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim)
46

  
47
     END SUBROUTINE CALC_XPRIM
48

  
49
  END INTERFACE
50

  
51
! Input
52
  REAL(KREAL),INTENT(IN) :: GeomCart(Nat,3)
53
  REAL(KREAL), INTENT(OUT) :: Values(NbVar)
54

  
55
  LOGICAL :: Debug
56
  INTEGER(KINT) :: I,J,K,NatT
57
  REAL(KREAL), ALLOCATABLE :: GeoCartLoc(:,:) ! (Nat+NbCom,3)
58
  REAL(KREAL), ALLOCATABLE :: x(:),y(:),z(:) ! Nat+NbCom
59
  REAL(KREAL) :: COG(3),Weight
60

  
61

  
62

  
63
  Debug=Valid('AnaGeom')
64

  
65
  If (Debug) Call Header("Entering AnalyzeGeom")
66

  
67
  if (debug) THEN
68
     WRITE(*,*) "AnalyzeGeom  - GeomCart"
69
     DO K=1,Nat
70
        WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
71
     END DO
72
  END IF
73

  
74
  NAtt=Nat+NbCom
75
  ALLOCATE(GeoCartLoc(Natt,3),x(Natt),y(Natt),z(Natt))
76
  GeoCartLoc(1:Nat,:)=GeomCart(:,:)
77
  CurBary => Bary
78

  
79
  DO I=1, NbCom
80
     COG=0.
81
     Weight=0.
82
     DO j=1,CurBary%ListAtoms(0)
83
        DO k=1,3
84
           COG(k)=COG(k)+GeomCart(CurBary%ListAtoms(j),k)*CurBary%Weights(j)
85
        END DO
86
        Weight=Weight+CurBary%Weights(j)
87
     END DO
88
     COG=COG/Weight
89
     DO k=1,3
90
        GeoCartLoc(Nat+i,k)=COG(k)
91
     END DO
92
  END DO
93
     
94
  Values=0.
95

  
96
  if (debug) THEN
97
     WRITE(*,*) "AnalyzeGeom before Calc_Xprim - GeomCartLoc"
98
     DO K=1,Nat+NbCom
99
        WRITE(*,*) K,GeoCartLoc(K,1:3)
100
     END DO
101
  END IF
102

  
103
  x = GeoCartLoc(:,1)
104
  y = GeoCartLoc(:,2)
105
  z = GeoCartLoc(:,3)
106

  
107
   Call Calc_XPrim(Natt,x,y,z,GeomList,NbVar,Values)
108

  
109
   DeALLOCATE(GeoCartLoc,x,y,z)
110

  
111
   if (debug) THEN
112
      WRITE(*,*) 'AnalyzeGeom: NbVar,Values',NbVar,Values
113
   END IF
114

  
115
  If (Debug) Call Header("Exiting AnalizeGeom")
116

  
117
END SUBROUTINE AnalyzeGeom
src/Io_module.f90 (revision 7)
8 8
  SAVE
9 9

  
10 10
  INTEGER(KINT) :: IOIN=11, IOOUT=12, IOCART=14
11
  INTEGER(KINT) :: IOGEOM=15
11
  INTEGER(KINT) :: IOGEOM=15, IODAT=16
12 12
  INTEGER(KINT), PARAMETER :: IOTMP=21,IOTMP2=22, IOTMP3=23
13 13
  INTEGER(KINT), PARAMETER :: IOERR=19
14

  
15
  TYPE Input_Line
16
     CHARACTER(LCHARS) :: Line
17
     TYPE (Input_Line), POINTER :: Prev
18
     TYPE (Input_Line), POINTER :: Next
19
  END TYPE Input_Line
20

  
21 14
  CHARACTER(SCHARS) :: RunMode
22 15

  
23 16
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
src/VarTypes.f90 (revision 7)
13 13
     INTEGER(KINT) :: At1,At2,At3,At4
14 14
! Value of the Coordinate
15 15
     REAL(KREAL) :: Value
16
!! ! This is the value to use to multiply the value to print it in a human readable format
17
!!     REAL(KREAL) :: PrintFactor
16 18
     INTEGER(KINT) :: SignDihedral
17 19
     TYPE (ListCoord), POINTER :: Next
18 20
  END TYPE ListCoord
19 21

  
22
  Type Barycenter
23
     INTEGER(KINT), ALLOCATABLE :: ListAtoms(:) ! Nat max
24
! Weights is not used for now 
25
     REAL(KREAL), ALLOCATABLE :: Weights(:) ! Nat max
26
     TYPE (Barycenter), POINTER :: Next
27
  END type BARYCENTER
20 28

  
29

  
30
  TYPE Input_Line
31
     CHARACTER(LCHARS) :: Line
32
     TYPE (Input_Line), POINTER :: Prev
33
     TYPE (Input_Line), POINTER :: Next
34
  END TYPE Input_Line
35

  
36

  
21 37
END MODULE VARTYPES
src/Path_module.f90 (revision 7)
54 54
! two initial geometries
55 55
  REAL(KREAL) :: BoxTol=0.5
56 56

  
57

  
58 57
  LOGICAL :: Freq, MW, Bohr, Renum, Hinv
59 58
! OptReac: if TRUE the reactants will be optimized. Default is FALSE
60 59
  LOGICAL :: OptReac
......
69 68
! CalEProd: if TRUE the products energy will be computed. Default is FALSE
70 69
! Not compatible with RunMode=Para
71 70
  LOGICAL :: CalcEProd
71
! AnaGeom: if TRUE the geometries are analyzed
72
  LOGICAL :: AnaGeom
73
! Name of the file to analyse geometries
74
  CHARACTER(LCHARS) :: AnaFile
75
! Nb: number of variables to monitor, including Centers of Mass
76
  INTEGER(KINT) :: Nb
77
! NbVar: number of geometrical variables to monitor
78
  INTEGER(KINT) :: NbVar
79
! NbCom: number of center of mass to create
80
  INTEGER(KINT) :: NbCom
81
  TYPE(ListCoord), POINTER :: GeomList,CurVar
82
  TYPE(Barycenter), POINTER :: Bary,CurBary
83
! Format to print the values
84
  CHARACTER(VLCHARS) :: FormAna
85
! Factor to use to print the values
86
  REAL(KREAL), ALLOCATABLE :: PrintGeomFactor(:) ! NbVar
87
! How to update the Hessian
72 88
  LOGICAL :: IniHup,HupNeighbour
73 89

  
74 90
  REAL(KREAL) :: Fact,FTan
src/Extrapol_int.f90 (revision 7)
15 15
  ! Default value of FAlign=.TRUE.
16 16
  use Path_module, only : NMaxPtPath, IntCoordI, Pi, IndZmat, XyzGeomF, &
17 17
       IntCoordF, IntTangent, Renum, Nom, Order, MassAt, SGeom, XyzGeomI, &
18
       Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms
18
       Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms,Align
19 19
  ! IndZmat(Nat,5)
20 20

  
21 21
  use Io_module
......
127 127
  END DO
128 128

  
129 129
  ! We align this geometry with the original one
130
  ! PFL 17/July/2006: only if we have more than 4 atoms.
131
  IF (Nat.GE.4) THEN
130
  ! PFL 17/July/2006: only if we have more than 4 atoms. 
131
! PFL 2013 Feb 27 ... or if the user asks for it
132
  IF ((Nat.GE.4).OR.Align) THEN
132 133
     ! The rotation matrix MRot has INTENT(OUT) in subroutine rotation_matrix(...),
133 134
     ! which is called in the CalcRmsd(...).
134 135
     ! PFL 24 Nov 2008 ->
......
263 264
             XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
264 265
     END DO
265 266

  
266
     IF (Nat.GE.4) THEN
267
     IF ((Nat.GE.4).OR.Align) THEN
267 268
        ! PFL 24 Nov 2008 ->
268 269
        ! If we have frozen atoms we align only those ones.
269 270
        if (NFroz.GT.0) THEN
......
288 289
           XYZTmp(I,J)=XyzTMP2(I,J)
289 290
        END DO
290 291
     END DO
292
     
291 293

  
292 294

  
293 295
     s=s+sqrt(ds)
......
338 340
           WRITE(IOOUT,'(1X,I5)') Nat
339 341
           WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
340 342
           ! PFL 17/July/2006: only if we have more than 4 atoms.
341
           IF (Nat.GE.4) THEN
343
           IF ((Nat.GE.4).OR.Align) THEN
342 344
              ! PFL 24 Nov 2008 ->
343 345
              ! If we have frozen atoms we align only those ones.
344 346
              if (NFroz.GT.0) THEN
......
382 384

  
383 385
  if (s>=0.9*dist) THEN
384 386
     if (debug) WRITE(*,*) "DBG Extrapol_int L383: Adding last geom"
385
     write(*,*) "Extrapol_int u,xgeom(NGeomI),s,dist,s-dist",u,xgeom(NGeomI),s,dist,s-dist
387
     if (debug) write(*,*) "Extrapol_int u,xgeom(NGeomI),s,dist,s-dist",u,xgeom(NGeomI),s,dist,s-dist
386 388
!     u=xgeom(NGeomI)
387 389
     s=s-dist
388 390

  
......
476 478
        WRITE(IOOUT,'(1X,I5)') Nat
477 479
        WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
478 480
        ! PFL 17/July/2006: only if we have more than 4 atoms.
479
        IF (Nat.GE.4) THEN
481
        IF ((Nat.GE.4).OR.Align) THEN
480 482
           ! PFL 24 Nov 2008 ->
481 483
           ! If we have frozen atoms we align only those ones.
482 484
           if (NFroz.GT.0) THEN
src/AnaPath.f90 (revision 7)
1
!================================================================
2
!
3
! This subroutine analyzes the geometries along a path
4
! and prints it to FileUnit
5
!================================================================
6

  
7
SUBROUTINE AnaPath(IOpt,FileUnit)
8

  
9

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

  
16

  
17
  IMPLICIT NONE
18

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

  
25

  
26
  INTERFACE
27
     function valid(string) result (isValid)
28
       CHARACTER(*), intent(in) :: string
29
       logical                  :: isValid
30
     END function VALID
31

  
32
     FUNCTION test_zmat(na,ind_zmat)
33

  
34
       USE Path_module
35

  
36
       logical :: test_zmat
37
       integer(KINT) :: na
38
       integer(KINT) :: ind_zmat(Na,5)
39
     END FUNCTION TEST_ZMAT
40

  
41

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

  
44
      Use VarTypes
45
      Use io_module
46

  
47
      implicit none
48

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

  
53
    END SUBROUTINE die
54

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

  
57
      Use VarTypes
58
      Use io_module
59

  
60
      implicit none
61

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

  
66
    END SUBROUTINE Warning
67

  
68
  END INTERFACE
69

  
70
  LOGICAL :: Debug, CalcDist
71
  INTEGER(KINT) :: I,J,K,Iat
72
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:),GeomCartPrec(:,:) ! Nat,3
73
  REAL(KREAL), ALLOCATABLE :: xRef(:),yRef(:),zRef(:) ! Nat
74
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
75
  REAL(KREAL) ::  ds
76
  REAL(KREAL) :: MRot(3,3), Rmsd  
77

  
78
  debug=valid("AnaPath")
79

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

  
82
  if (debug) WRITE(*,*) "DBG AnaPaht PrintGeomFactor:",PrintGeomFactor
83
  ALLOCATE(GeomCart(Nat,3))
84
  ALLOCATE(Values(NbVar))
85

  
86
  CalcDist=.FALSE.
87

  
88
  WRITE(FileUnit,'("# Iteration ",I5)') Iopt
89

  
90
  DO J=1,Nat
91
     Iat=J
92
     IF (renum) Iat=Order(J)
93
     GeomCart(J,:)=XyzGeomF(1,:,Iat)
94
  END DO
95

  
96

  
97
! Do we know the curvilinear values ?
98
  IF (MOD(IOpt,IReparam)/=0) THEN
99
     CalcDist=.TRUE.
100
     SGeom=0.
101
     ALLOCATE(xRef(Nat),yRef(Nat),zref(Nat),GeomCartPrec(Nat,3))
102
     xRef=GeomCart(:,1)
103
     yRef=GeomCart(:,2)
104
     zRef=GeomCart(:,3)
105
     GeomCartPrec=GeomCart
106
  END IF
107

  
108
  if (debug) THEN
109
     WRITE(*,*) "AnaPath  - GeomCart - I=",1
110
     DO K=1,Nat
111
        WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
112
     END DO
113
     WRITE(*,*) "AnaPath  - XYzGeomF(I,:,:) - I=",1
114
     DO K=1,Nat
115
        WRITE(*,'(1X,I5,3(1X,F15.8))') K,XyzGeomF(1,:,K)
116
     END DO
117
  END IF
118

  
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

  
123
  DO I=2,NGeomF
124
     DO J=1,Nat
125
        Iat=J
126
        IF (renum) Iat=Order(J)
127
        GeomCart(J,:)=XyzGeomF(I,:,Iat)
128
     END DO
129

  
130
     If (CalcDist) THEN
131
! PFL 2013 Feb 26
132
! For now, I do _NOT_ align the geometries
133
        if (debug) THEN
134
           WRITE(*,*) "AnaPath  - GeomCart avant align - I=",I
135
           DO K=1,Nat
136
              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
137
           END DO
138
        END IF
139

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

  
143
! If we have frozen atoms we align only those ones.
144
           if (NFroz.GT.0) THEN
145
              Call AlignPartial(Nat,xRef,yRef,zRef,                 &
146
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),       &
147
                   FrozAtoms,MRot)
148
           ELSE 
149
              Call  CalcRmsd(Nat, xRef,yRef,zRef,               &
150
                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),   &
151
                   MRot,rmsd,.TRUE.,.TRUE.)
152
           END IF
153
              ! <- PFL 24 Nov 2008
154
           
155
        END IF
156

  
157
        if (debug) WRITE(*,*) "Mass:",MassAt(1:Nat)
158
        ds=0.
159
        DO J=1,Nat
160
           Iat=J
161
!        IF (renum) Iat=Order(J)
162
           if (debug) WRITE(*,*) "Mass(J):",J,MassAt(Iat)
163
            
164
           do K=1,3
165
              ds=ds+MassAt(Iat)*(GeomCartPrec(J,K)-GeomCart(J,K))**2
166
           END DO
167
           if (debug) WRITE(*,*) "DBG DS:",J,ds
168
        END DO
169
        ds=sqrt(ds)
170
        if (debug) THEN
171
           WRITE(*,*) "AnaPath  - GeomCartPrec - I=",I
172
           DO K=1,Nat
173
              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCartPrec(K,1:3)
174
           END DO
175
           WRITE(*,*) "DBG Anapaht, ds=",ds
176
        END IF
177

  
178
        SGeom(I)=SGeom(I-1)+ds
179
        GeomCartPrec=GeomCart
180
        xRef=GeomCart(:,1)
181
        yRef=GeomCart(:,2)
182
        zRef=GeomCart(:,3)
183
        
184
     END IF
185

  
186
     if (debug) THEN
187
        WRITE(*,*) "AnaPath  - GeomCart - I=",I
188
        DO K=1,Nat
189
           WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)
190
        END DO
191
     END IF
192

  
193
     Call AnalyzeGeom(GeomCart,Values)
194
     WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*au2kcal,Energies(i)
195
  END DO
196

  
197
  WRITE(FileUnit,*) ""
198
  WRITE(FileUnit,*) ""
199

  
200
  DeAllocate(GeomCart,Values)
201
  If (CalcDist) Deallocate(GeomCartPrec,xRef,yRef,zRef)
202

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

  
205
END SUBROUTINE AnaPath
206

  
src/Makefile (revision 7)
232 232
      ReadInput_vasp.f90 \
233 233
      ReadInput_mopac.f90 \
234 234
      ReadInput_siesta.f90 \
235
      ReadAnaList.f90 \
236
      PrintAnaList.f90 \
237
      AnaPath.f90 \
238
      AnalyzeGeom.f90 \
235 239
      Calc_zmat.f90 \
236 240
      Calc_zmat_frag.f90 \
237 241
      Calc_zmat_constr_frag.f90 \
src/Calc_Xprim.f90 (revision 7)
80 80

  
81 81
  debug=valid("Calc_Xprim")
82 82
  debugPFL=valid("BakerPFL")
83
  if (debug) WRITE(*,*) "============= Entering Cal_XPrim =============="
84 83

  
84
  if (debug) Call Header("Entering Cal_XPrim")
85 85

  
86

  
86 87
  IF (debug) THEN
87 88
     WRITE(*,*) "DBG Calc_Xprim, geom"
88 89
     DO I=1,Nat
89 90
        WRITE(*,'(1X,I5,3(1X,F15.3))') I,X(I),y(i),z(i)
90 91
     END DO
91
     WRITE(*,*) "XPrimRef"
92
     WRITE(*,'(15(1X,F10.6))') XPrimRef
92
     IF (Present(XPrimRef)) THEN
93
        WRITE(*,*) "XPrimRef"
94
        WRITE(*,'(15(1X,F10.6))') XPrimRef
95
     END IF
96
     WRITE(*,*) "NPrim:",NPrim
93 97
  END IF
94 98

  
95 99

  
96
!  WRITE(*,*) "Coordinate:",associated(Coordinate)
100
  WRITE(*,*) "Coordinate:",associated(Coordinate),coordinate%Type
97 101

  
98
     ScanCoord=>Coordinate
102
     ScanCoord => Coordinate
99 103

  
100
!  WRITE(*,*) "ScanCoord:",associated(ScanCoord),ScanCoord%Type
104
  WRITE(*,*) "ScanCoord:",associated(ScanCoord),ScanCoord%Type
101 105

  
102 106
!     WRITE(*,*) "coucou"
103 107
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
......
120 124
                vx4,vy4,vz4,norm4)
121 125
           Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
122 126
                vx5,vy5,vz5,norm5)
123

  
124
                 DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
125
                         vx2,vy2,vz2,norm2)
126
                 Xprimitive(I) = DiheTmp*Pi/180.
127
           
128
           DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
129
                vx2,vy2,vz2,norm2)
130
           Xprimitive(I) = DiheTmp*Pi/180.
127 131
! PFL 15th March 2008 ->
128 132
! I think that the test on changes less than Pi should be enough
129 133
!! We treat large dihedral angles differently as +180=-180 mathematically and physically
......
139 143
!!!! <- PFL 15th March 2008
140 144
! Another test... will all this be consistent ??? 
141 145
! We want the shortest path, so we check that the change in dihedral angles is less than Pi:
142
                  IF (Present(XPrimRef)) THEN
143
                     IF (abs(Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
144
                       if (debug) THEN 
145
                         WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I)
146
                         WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi
147
                       END IF
148
                        if ((Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
149
                           Xprimitive(I)=Xprimitive(I)-2.d0*Pi
150
                        ELSE
151
                           Xprimitive(I)=Xprimitive(I)+2.d0*Pi
152
                        END IF
153
                     END IF
154
                       if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi
155
                  END IF
156
        END SELECT
146
           IF (Present(XPrimRef)) THEN
147
              IF (abs(Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
148
                 if (debug) THEN 
149
                    WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I)
150
                    WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi
151
                 END IF
152
                 if ((Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
153
                    Xprimitive(I)=Xprimitive(I)-2.d0*Pi
154
                 ELSE
155
                    Xprimitive(I)=Xprimitive(I)+2.d0*Pi
156
                 END IF
157
              END IF
158
              if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi
159
           END IF
160
        CASE DEFAULT
161
           WRITE(*,*) "Scancoord%type unknown for I",I,scancoord%type
162
      END SELECT
157 163
        ScanCoord => ScanCoord%next
158 164
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
159 165

  
......
161 167

  
162 168
  IF (debug) THEN
163 169
     WRITE(*,*) "DBG Calc_Xprim Values"
170
     WRITE(*,*) "Found ",I," primitives"
164 171
     
165 172
     ScanCoord=>Coordinate
166 173
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
......
179 186
        ScanCoord => ScanCoord%next
180 187
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
181 188
  END IF
182
  if (debug) WRITE(*,*) "============= Cal_XPrim OVER =============="
189

  
190
  if (debug) Call Header(" Cal_XPrim OVER ")
191

  
183 192
END SUBROUTINE Calc_Xprim
src/bib_oper.f90 (revision 7)
104 104
                               v2x,v2y,v2z, &
105 105
                               v3x,v3y,v3z,norm3)
106 106

  
107
        use Path_module, only :  Pi,KINT, KREAL
107
        use Path_module, only :   KREAL
108 108

  
109 109
        real(KREAL) ::  v1x,v1y,v1z ! what do you do with norm1, norm2???
110 110
        real(KREAL) ::  v2x,v2y,v2z
src/WriteInput_siesta.f90 (revision 7)
1
 SUBROUTINE ReadInput_siesta
1
 SUBROUTINE WriteInput_siesta(GeomCart,FileUnit)
2 2

  
3
! This routine reads an input template for Siesta
3
! This routine writes an input for Siesta
4 4

  
5 5
  use VarTypes
6 6
  use Path_module
......
14 14
       logical                  :: isValid
15 15
     END function VALID
16 16

  
17
    FUNCTION SearchInput(Input,String,Line,Clean)  Result (Found)
18 17

  
19
      Use Vartypes
20
      use io_module
21
! Input
22
      TYPE (Input_line), POINTER, INTENT(IN) :: Input
23
      CHARACTER(*), INTENT(IN) :: String
24
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean
18
     SUBROUTINE WriteList(Input,Unit)
19
        ! This routine reads an input template for Siesta
25 20

  
26
! Output
27
      TYPE (Input_line), POINTER, INTENT(OUT) :: Line
21
       use VarTypes
22
       use Io_module
28 23

  
29
      LOGICAL :: Found
24
       IMPLICIT NONE
30 25

  
31
    END FUNCTION SearchInput
26
     ! Input
27
       TYPE (Input_line), POINTER, INTENT(IN) :: Input
28
       INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit
29
       
30
     END SUBROUTINE WriteList
32 31

  
33
    SUBROUTINE die(routine, msg, file, line, unit)
34

  
35
      Use VarTypes
36
      Use io_module
37

  
38
      implicit none
39
!--------------------------------------------------------------- Input Variables
40
      character(len=*), intent(in)           :: routine, msg
41
      character(len=*), intent(in), optional :: file
42
      integer(KINT), intent(in), optional      :: line, unit
43

  
44
    END SUBROUTINE die
45

  
46
    SUBROUTINE Warning(routine, msg, file, line, unit)
47

  
48
      Use VarTypes
49
      Use io_module
50

  
51
      implicit none
52
!--------------------------------------------------------------- Input Variables
53
      character(len=*), intent(in)           :: routine, msg
54
      character(len=*), intent(in), optional :: file
55
      integer(KINT), intent(in), optional      :: line, unit
56

  
57
    END SUBROUTINE Warning
58

  
59 32
  END INTERFACE
60 33

  
61

  
62
  CHARACTER(132) ::  Line,LineUp,Line2
63
  INTEGER(KINT) :: LineL, Idx
64
  INTEGER(KINT) :: ISpec, IAt,I
65
  INTEGER(KINT) :: IoRead, ITmp, JTmp
66
  REAL(KREAL) :: Xtmp, Ytmp, Ztmp
34
!Input
35
 ! Geometry in cartesian coordinates
36
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
37
! Unit to write to
38
  INTEGER(KINT), INTENT(IN) :: FileUnit
67 39
  
68 40
  LOGICAL :: Debug
69
  LOGICAL :: FSpecies, FCoord
41
  INTEGER(KINT) :: I,Iat
70 42

  
71
  TYPE(Input_Line), POINTER :: Search,Bla
43
  Debug=Valid("WriteInput").OR.Valid("WriteInput_siesta")
72 44

  
73
  Debug=Valid("readinput").OR.Valid("readinput_siesta")
45
  if (debug) Call Header("Entering WriteInput_Siesta")
74 46

  
75
 if (debug) Call Header("Entering ReadInput_Siesta")
76 47

  
77
  ! We read the Siesta input file
48
  Call WriteList(Siesta_input,Unit=FileUnit)
78 49

  
79
     ALLOCATE(Siesta_Input)
80
     NULLIFY(Siesta_Input%next)
81
     NULLIFY(Siesta_Input%prev)
82
     Current => Siesta_Input
50
  WRITE(FileUnit,*) 
83 51

  
84
     FSpecies=.FALSE.
85
     FCoord=.FALSE.
52
  WRITE(FileUnit,'(1X,A)')  '%block AtomicCoordinatesAndAtomicSpecies' 
86 53

  
87
     READ(IOIN,'(A)',iostat=IoRead) Line
54
  DO I=1,Nat
55
     If (renum) THEN
56
        Iat=Order(I)
57
        WRITE(FileUnit,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(Iat,:)/Siesta_Unit_Read, IdxSpecies(Iat),TRIM(Siesta_Paste(I))
58
     ELSE
59
        Iat=OrderInv(I)
60
        WRITE(FileUnit,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(I,:)/Siesta_Unit_Read, IdxSpecies(Iat), TRIM(Siesta_Paste(Iat))
61
     END IF
62
  END DO
88 63

  
89
     DO WHILE (IoRead==0)
90
        Line=AdjustL(Line)
91
        if (debug) WRITE(*,*) 'Line:', Line
92
        current%Line=TRIM(Line)
93
        ALLOCATE(current%next)
94
        Current%next%prev => Current
95
        Current => Current%next
96
        Nullify(Current%next)
64
  WRITE(FileUnit,'(1X,A)')  '%endblock AtomicCoordinatesAndAtomicSpecies'
65
  WRITE(FileUnit,*) 
97 66

  
98
        READ(IOIN,'(A)',iostat=IoRead) Line
99
     END DO
100 67

  
101
! We analyse the input
102 68

  
103
! We look for the NumberofAtoms
104
        If (SearchInput(Siesta_Input,"NUMBEROFATOMS",Search,Clean=".-_")) THEN
105
           Line=AdjustL(Search%line)
106
           Idx=Index(Line," ")
107
           Line2=Trim(AdjustL(Line(Idx+1:)))
108
           Read(Line2,*) IAt
109
           if (Iat/=Nat) THEN
110
              Call Die('ReadInput_siesta','Nat in FDF sample different from  Nat Path input', Unit=IOOUT)
111
           END IF
112
        ELSE
113
! There is no atom number defined !!!
114
           Call Warning('ReadInput_siesta','No NumberOfAtoms in FDF sample', Unit=IOOUT)
115
           WRITE(current%Line,'(1X,"NumberOfAtoms ",I5)') Nat
116
           ALLOCATE(current%next)
117
           Current%next%prev => Current
118
           Current => Current%next
119
           Nullify(Current%next)
120
         
121
        END IF
122

  
123
! We look for the NumberofSpecies
124
        IF (SearchInput(Siesta_Input,"NUMBEROFSPECIES",Search,Clean=".-_")) THEN
125
           Line=AdjustL(Search%line)
126
           Idx=Index(Line," ")
127
           Line2=Trim(AdjustL(Line(Idx+1:)))
128
           Read(Line2,*) Siesta_NbSpecies
129
        END IF
130

  
131
! We look for SystemLabel
132
        If (SearchInput(Siesta_Input,"SYSTEMLABEL",Search,Clean=".-_")) THEN
133
           Line=AdjustL(Search%line)
134
           Idx=Index(Line," ")
135
           Siesta_Label=Trim(adjustl(Line(Idx+1:)))
136
        ELSE
137
           Siesta_label='siesta'
138
        END IF
139

  
140
! We look for the ChemicalSpeciesLabel block
141
        IF (SearchInput(Siesta_Input,"CHEMICALSPECIESLABEL",Search,Clean=".-_")) THEN
142
           ALLOCATE(ListSpecies(Siesta_NbSpecies))
143
           ALLOCATE(Siesta_SpeciesName(Siesta_NbSpecies))
144
           ALLOCATE(IdxSpecies(NAt))
145
           DO I=1,Siesta_NbSpecies
146
              Search => Search%next
147
              Line=AdjustL(Search%Line)
148
              Read(Line,*) ITmp
149
              Idx=Index(Line,' ')
150
              Line=AdjustL(Line(Idx+1:))
151
              Read(Line,*) Ztmp
152
              ListSpecies(ITmp)=ZTmp
153
              Idx=Index(Line,' ')
154
              Line=AdjustL(Line(Idx+1:))
155
              Idx=Index(Line,' ')
156
              Siesta_SpeciesName(ITmp)=Line(1:Idx-1)
157
           END DO
158
           if (Debug) THEN
159
              WRITE(*,*) 'Found ',Siesta_NbSpecies,' species'
160
              DO I=1,Siesta_NbSpecies
161
                 WRITE(*,*) I, ListSpecies(I),TRIM(Siesta_speciesName(I))
162
              END DO
163
        END IF
164

  
165

  
166
! We look for the AtomicCoordinatesAndAtomicSpecies  block
167
        IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESANDATOMICSPECIES",Search,Clean=".-_")) THEN
168
           ALLOCATE(Siesta_Paste(Nat))
169
           Current=>Search
170
           DO I=1,NAt
171
              Search => Search%next
172
              Read(Search%line,*) XTmp,YTmp,ZTmp,IdxSpecies(I)
173
! We save everything but the x,y,z, and species description              
174
              Line=AdjustL(search%line)
175
! We skip x
176
              Idx=Index(Line,' ')
177
              Line=AdjustL(Line(Idx+1:))
178
! we skip y
179
              Idx=Index(Line,' ')
180
              Line=AdjustL(Line(Idx+1:))
181
! we skip z
182
              Idx=Index(Line,' ')
183
              Line=AdjustL(Line(Idx+1:))
184
! we skip species
185
              Idx=Index(Line,' ')
186
              Line=AdjustL(Line(Idx+1:))
187
              Siesta_Paste(I)=TRIM(Line)
188
           END DO
189
! We will now delete this block from our input sample as it will then be
190
! written directly by Opt'n Path
191
! Search%next point on %endblock
192
           search => search%next
193
           IF (ASSOCIATED(Search%next)) THEN
194
! we are not at the end of the input file
195
              if (ASSOCIATED(Current%Prev)) THEN
196
                 Search%next%prev => current%prev
197
              ELSE
198
! the coordinate block is at the begining of the input
199
                 siesta_input => Search%next
200
                 Nullify(Siesta_Input%Prev)
201
              END IF
202
           ELSE
203
! the coordinate block is the last part of the input
204
              nullify(current%prev%next)
205
           END IF
206
              DO I=1,Nat+2
207
                 bla=>current
208
                 current=> current%next
209
                 deallocate(bla)
210
              END DO
211
        ELSE
212
           IF (SearchInput(Siesta_Input,"ZMATRIX",Search,Clean=".-_")) THEN
213
              call Die('ReadInput_Siesta','For now, I need the full block'//&
214
                   'AtomicCoordinatesAndAtomicSpecies, ZMatrix not yet handled.')
215
           ELSE
216
              call Die('ReadInput_Siesta','For now, I need the full block' //&
217
              'AtomicCoordinatesAndAtomicSpecies.')
218
           END IF
219
        END IF
220

  
221
        IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN
222
           Line=Adjustl(Search%Line)
223
           Idx=Index(Line,' ')
224
           Line=AdjustL(Line(Idx+1:))
225
           Call UpCase(Line)
226
           SELECT CASE (Line) 
227
              CASE ('ANG','NOTSCALEDCARTESIANANG')
228
                 Siesta_Unit_Read=1.d0
229
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
230
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
231
                 Siesta_Unit_Read=a0
232
            END SELECT
233
         END IF
234

  
235
        IF (SearchInput(Siesta_Input,"ATOMICCOORFORMATOUT",Search,Clean=".-_")) THEN
236
           Line=Adjustl(Search%Line)
237
           Idx=Index(Line,' ')
238
           Line=AdjustL(Line(Idx+1:))
239
           Call UpCase(Line)
240
           SELECT CASE (Line) 
241
              CASE ('ANG','NOTSCALEDCARTESIANANG')
242
                 Siesta_Unit_Write=1.d0
243
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
244
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
245
                 Siesta_Unit_Write=Unit
246
            END SELECT
247
         END IF
248

  
249
         IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN
250
           Line=Adjustl(Search%Line)
251
! We discar the label
252
           Idx=Index(Line,' ')
253
           Line=AdjustL(Line(Idx+1:))
254
! we read the value
255
           Read(Line,*) Siesta_LatticeConstant
256
! We discard the value
257
           Idx=Index(Line,' ')
258
           Line=AdjustL(Line(Idx+1:))
259
! We read the unit
260
           Call UpCase(Line)
261
           SELECT CASE (Line) 
262
              CASE ('ANG')
263
                 Siesta_Lat_Unit=1.d0
264
              CASE ('BOHR')
265
                 Siesta_Lat_Unit=Unit
266
            END SELECT
267
! for now!
268
            Call Die('ReadInput_siesta:LatticeConstant',"For now, periodic calculations are NOT possible in Opt'n Path")
269

  
270
         END IF
271

  
272
         IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN
273

  
274
            Call Die('ReadInput_siesta:LatticeVectors',"For now, periodic calculations are NOT possible in Opt'n Path")
275

  
276
         END IF
277

  
278
         IF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN
279

  
280
            Call Die('ReadInput_siesta:LatticeParameters',"For now, periodic calculations are NOT possible in Opt'n Path")
281

  
282
         END IF
283

  
284
         IF (SearchInput(Siesta_Input,"SUPERCELL",Search,Clean=".-_")) THEN
285

  
286
            Call Die('ReadInput_siesta:SuperCell',"SuperCell  NOT possible in Opt'n Path")
287

  
288
         END IF
289

  
290
         Call Die('Readinput_siesta:fin',' Lecture finie')
291

  
292
     if (debug) Call Header("Exiting ReadInput_Siesta")
293

  
294
   END SUBROUTINE READINPUT_SIESTA
69
  if (debug) Call Header("Exiting WriteInput_Siesta")
70
  
71
END SUBROUTINE WRITEINPUT_SIESTA
src/Extrapol_mixed.f90 (revision 7)
122 122
     if (s>=dist) THEN
123 123
        s=s-dist
124 124
        IdxGeom=IdxGeom+1
125
        SGeom(IdxGeom)=s+IdxGeom*dist
125 126
        XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
126 127

  
127 128
        IntCoordF(IdxGeom,:)=IntCoordTmp
......
154 155

  
155 156
        END IF
156 157
     END IF
157
  ENDDO
158
  END DO
158 159

  
159 160

  
160 161
  if ((s>=0.1*dist).AND.(IdxGeom.EQ.NGeomF)) THEN
......
168 169
  IdxGeom=NGeomF
169 170

  
170 171
! We have to add the last geometry. We copy the last geom of Initial geometries.
171
     IntCoordF(IdxGeom,:)=IntCoordI(NGeomI,:)
172
     IntTangent(IdxGeom,:)=DerInt
172
  IntCoordF(IdxGeom,:)=IntCoordI(NGeomI,:)
173
  IntTangent(IdxGeom,:)=DerInt
173 174

  
174 175
! we convert it to cartesian geom
175 176
  call Mixed2Cart(Nat,IndZmat,IntCoordF(IdxGeom,1),XyzTmp2(1,1))
src/ReadAnaList.f90 (revision 7)
1
 SUBROUTINE ReadAnaList
2
! This routines read a list of geometrical variables to monitor
3
! This is inspired from Xyz2Path (that was inspired by Xyz2scan ...)
4

  
5

  
6
  use VarTypes
7
  use Path_module
8
  use Io_module
9

  
10
  IMPLICIT NONE
11

  
12

  
13
  INTERFACE
14
     function valid(string) result (isValid)
15
       CHARACTER(*), intent(in) :: string
16
       logical                  :: isValid
17
     END function VALID
18

  
19
     SUBROUTINE die(routine, msg, file, line, unit)
20

  
21
       Use VarTypes
22
       Use io_module
23

  
24
       implicit none
25

  
26
       character(len=*), intent(in)           :: routine, msg
27
       character(len=*), intent(in), optional :: file
28
       integer(KINT), intent(in), optional      :: line, unit
29

  
30
     END SUBROUTINE die
31

  
32
  END INTERFACE
33

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

  
38
  CHARACTER(LCHARS) :: Line
39

  
40
  Debug=Valid('ReadAnaList')
41
  If (Debug) Call Header("Entering ReadAnaList")
42

  
43
  Allocate(Bary)
44
  Nullify(Bary%Next)
45
  CurBary => Bary
46

  
47
  ALLOCATE(GeomList)
48
  Nullify(GeomList%Next)
49
  CurVar => GeomList
50

  
51
  NbCom=0
52
  NbVar=0
53
  FormAna='(1X'
54
  DO I=1, Nb
55
     READ(IOIN,'(A)') Line
56
     Line=AdjustL(Line)
57
     Call Upcase(Line)
58

  
59
     SELECT CASE (Line(1:1))
60
        CASE ('B') 
61
! this is a bond
62
           CurVar%Type="BOND"
63
           Idx=Index(Line," ")
64
           Line=Line(Idx+1:)
65
           Read(Line,*) At1, At2
66
           CurVar%At1=At1
67
           CurVar%At2=At2
68
           CurVar%Value=0.
69
!           CurVar%PrintFactor=1.
70
           Curvar%SignDihedral=1
71
           if (debug) THEN
72
              WRITE(*,'("# b ",2I3)') At1,At2
73
           END IF
74
           FormAna=TRIM(FormAna) //',1X,F7.3'
75
           Allocate(CurVar%Next)
76
           CurVar => CurVar%Next
77
           Nullify(CurVar%Next)
78
        CASE ('A') 
79
! this is a valence angle
80
           CurVar%Type="ANGLE"
81
           Idx=Index(Line," ")
82
           Line=Line(Idx+1:)
83
           Read(Line,*) At1, At2,At3
84
           CurVar%At1=At1
85
           CurVar%At2=At2
86
           CurVar%At3=At3
87
           CurVar%Value=0.
88
!           CurVar%PrintFactor=180./Pi
89
           Curvar%SignDihedral=1
90
           if (debug) THEN
91
              WRITE(*,'("# a ",4(I3))') At1,At2,At3
92
           END IF
93
           FormAna=TRIM(FormAna) //',1X,F7.2'
94
           Allocate(CurVar%Next)
95
           CurVar => CurVar%Next
96
           Nullify(CurVar%Next)
97
        CASE ('D') 
98
! this is a dihedral
99
           CurVar%Type="DIHEDRAL"
100
           Idx=Index(Line," ")
101
           Line=Line(Idx+1:)
102
           Read(Line,*) At1, At2,At3,At4
103
           CurVar%At1=At1
104
           CurVar%At2=At2
105
           CurVar%At3=At3
106
           CurVar%At4=At4
107
           CurVar%Value=0.
108
!           CurVar%PrintFactor=180./Pi
109
           Curvar%SignDihedral=1
110
           if (debug) THEN
111
              WRITE(*,'("# d ",4(I3))') At1,At2,At3,At4
112
           END IF
113
           Allocate(CurVar%Next)
114
           CurVar => CurVar%Next
115
           Nullify(CurVar%Next)
116
        CASE ('C') 
117
           NbCom=NbCom+1
118
           Idx=Index(Line," ")
119
           Line=Line(Idx+1:)
120
           READ(Line,*) At1
121
           Allocate(CurBary%ListAtoms(0:At1))
122
           Allocate(CurBary%Weights(1:At1))
123
           CurBary%Weights=1.
124
           CurBary%ListAtoms(0)=At1
125
           Idx=Index(Line," ")
126
           Line=Line(Idx+1:)
127
           READ(Line,*) CurBary%ListAtoms(1:At1)
128
           if (debug) THEN
129
              WRITE(*,'("# c ",I4,20(I3))') At1, &
130
                   (CurBary%ListAtoms(j),j=1,At1)
131
           END IF
132
           FormAna=Trim(FormAna) //',1X,F7.2'
133

  
134
           Allocate(CurBary%Next)
135
           CurBary => CurBary%Next
136
           Nullify(CurBary%Next)
137

  
138
        CASE Default
139
           Call Die('ReadAnaList','Variable not recognized: ' //Line, Unit=IOOUT)
140
     END SELECT
141

  
142
  END DO
143

  
144
  NbVar=Nb-NbCom
145

  
146
  CurVar => GeomList
147
  I=0
148

  
149
  Allocate(PrintGeomFactor(NbVar))
150
  DO WHILE (associated(CurVar%Next))
151
     I=I+1
152
     SELECT CASE (CurVar%TYPE)
153
        CASE ('ANGLE','DIHEDRAL')
154
           PrintGeomFactor(I)=180./pi
155
        CASE DEFAULT
156
           PrintGeomFactor(I)=1.
157
     END SELECT
158
     CurVar => CurVar%Next
159
  END DO
160
  
161
  FormAna=Trim(FormAna) // ',1X,F7.2,1X,F15.6)'
162

  
163

  
164
!  if (debug) WRITE(*,*) "NbVar, NbCom,Nb",NbVar,NbCom,Nb
165

  
166
  if (debug) Call Header(" Exiting ReadAnaList")
167

  
168
END SUBROUTINE ReadAnaList
src/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 
... Ce différentiel a été tronqué car il excède la taille maximale pouvant être affichée.

Formats disponibles : Unified diff