Revision 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/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/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/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/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/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/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/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/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/Path.f90 (revision 7)
212 212
! We add Siesta (and start thinking of add CP2K)
213 213
!  as a new "energy engine" code
214 214
! TO DO to go to 1.5: finish cleaning...
215
!
215 216
! 2013 Feb Opt'n Path v1.48
216 217
! We add some keyword for more output for Geometry Optimizations
217 218
! GeomFile: name of the file to print out the geometries and their energy
218 219
! as Xyz file. (only format for now)
220
!
221
! 2013 Feb opt'n Path v1.49
222
! We add the possibility to analyze the geometries in terms
223
! of Bonds, Angles, dihedrals
224
! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path)
225
! This is basicaly the merging of AnaPath into Opt'n Path... at last !
226
! This is done by adding a Flag into the Path namelist:
227
! AnaGeom: logical, if true Path looks for the AnaList namelist
228
! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path
229
! will read it and print the values of the variable in a Path.dat file
230
! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. 
219 231

  
220 232
PROGRAM PathOpt
221 233

  
......
249 261
  LOGICAL :: FChkFrozen
250 262

  
251 263
  ! Energies for all points along the path at the previous iteration
252
  REAL(KREAL), ALLOCATABLE ::  EPathold(:) ! NGeomF, where it is deallocated: Prakash
264
  REAL(KREAL), ALLOCATABLE ::  EPathold(:) ! NGeomF
253 265
  ! Maximum step allowed for this geometry
254
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF, where it is deallocated: Prakash 
266
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF
255 267

  
256 268
! these are used to read temporary coordinates
257 269
  LOGICAL :: FFrozen,FCart
......
272 284

  
273 285
!  INTEGER(KINT), EXTERNAL :: Iargc
274 286
  INTEGER(KINT), EXTERNAL ::  ConvertNumAt
287
  INTEGER(KINT) :: IoS
275 288

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

  
332 345
  Namelist/cartlist/list
333 346
  Namelist/frozenlist/list
347
  Namelist/analist/nb
334 348

  
335 349

  
336 350
  Flag_Opt_Geom = .FALSE. ! Initialized.
......
368 382
        WRITE(*,*) "Input: string that indicates the type of the input geometries."
369 383
        WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp"
370 384
        WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
371
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac, Test, Siesta or Ext"
372
        WRITE(*,*) "      In case of a Gaussian or Mopac calculations, input must be set to Cart."
385
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac, TurboMole, Siesta, Test or Ext"
386
        WRITE(*,*) "   *   In case of a Gaussian or Mopac calculations, input must be set to Cart."
373 387
        WRITE(*,*) "      One example of a gaussian/mopac input should be added at the end of the"
374
        WRITE(*,*) "      input file.See example file Test_G03.path or Test_Mopac.path"
375
        WRITE(*,*) "      In the case of a VASP calculation, if input is set to Cart, then"
376
        WRITE(*,*) "    the preamble of a VASP calculation should be added at the end of"
377
        WRITE(*,*) "    the input file. See example file Test_VASP_cart.path"
378
        WRITE(*,*) "    In the case of a VASP calculation, one should also give value of the "
388
        WRITE(*,*) "     input file. See example file Test_G03.path or Test_Mopac.path"
389
        WRITE(*,*) "   *   In the case of a VASP calculation, if input is set to Cart, then"
390
        WRITE(*,*) "     the preamble of a VASP calculation should be added at the end of"
391
        WRITE(*,*) "     the input file. See example file Test_VASP_cart.path"
392
        WRITE(*,*) "       In the case of a VASP calculation, one should also give value of the "
379 393
        WRITE(*,*) "    Runmode variable"
394
        WRITE(*,*) "   *   In the case of a SIESTA calculation, an example of a Siesta input file"
395
        WRITE(*,*) "     should be added at the end of the input file. See Test_Siesta.path"
380 396
        WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy"
381 397
        WRITE(*,*) "       and gradient of the whole path or not. If one wants to use VASP,"
382 398
        WRITE(*,*) "       Runmode must be set to PARA."
......
413 429
        WRITE(*,*) "Smax: Maximum length of a step during path optimization"
414 430
        WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
415 431
        WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken"
416
        WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"
432
        WRITE(*,*) "FTan: When moving the path, this gives the proportion of the displacement tangent to the path"
417 433
        WRITE(*,*) "      that is kept. FTan=1. corresponds to the full displacement, "
418
        WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path."
434
        WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path. Default: Ftan=0."
419 435
        WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization."
420 436
        WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path."
421 437
        WRITE(*,*) "         This option indicates the first step where spline interpolation is used."
......
451 467
        WRITE(*,*) "      If True, a &cartlist namelist containing the list of cart atoms must be given."
452 468
        WRITE(*,*) "      By default, only frozen atoms are described in cartesian coordinates."
453 469
        WRITE(*,*) ""
470
        WRITE(*,*) "AnaGEom: If TRUE, Opt'n Path will create a file .dat for geometries analysis."
471
        WRITE(*,*) "        If True, Opt'n Path will look for the AnaList namelist after the Path Namelist."
472
        WRITE(*,*) "AnaList: list of variables for geometry analysis. If present and if AnaGeom=T then"
473
        WRITE(*,*) "     Opt'n Path will read it and print the values of the variable in a .dat file"
474
        WRITE(*,*) "    If AnaGeom is T but Analist is not given, then Path.dat just contains the energy." 
475
        WRITE(*,*) "Analist contains the number of variables to monitor: Nb, and is followed by the description"
476
        WRITE(*,*) "of the variables among:"
477
        WRITE(*,*) "b(ond) At1 At2"
478
        WRITE(*,*) "a(ngle) At1 At2 At3"
479
        WRITE(*,*) "d(ihedral) At1 At2 At3 At4"
480
        WRITE(*,*) "c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added"
481
        WRITE(*,*) "at the end of the real atoms of the system"
482
        WRITE(*,*) ""
454 483
        WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
455 484
        WRITE(*,*) "        If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. "
456 485
        WRITE(*,*) "       It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]"
......
470 499
  END SELECT
471 500

  
472 501

  
473

  
474 502
  ! We initiliaze variables
475 503
  Pi=dacos(-1.0d0)
476 504
  Nat=1
......
507 535
  CalcEReac=.FALSE.
508 536
  CalcEProd=.FALSE.
509 537
  EReac=0.d0
510
  EProd=0.d0
538
  EProd=0.d0 
511 539
  OptGeom=-1
512 540
  GeomFile="EMPTY"
541
  AnaGeom=.FALSE.
542
  AnaFile="EMPTY"
543
  Nb=0
544
  NbCom=0
513 545
  PathOnly=.FALSE.
514 546
  Prog='EMPTY'
515 547
  ProgExe='EMPTY'
......
580 612
     END IF
581 613
  END IF
582 614

  
583
  If (COORD.EQ.'XYZ') THEN 
615
  If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN 
584 616
     COORD='CART'
585 617
  END IF
586 618

  
......
651 683
    WRITE(*,*) "Input has been set to the default: ",INPUT
652 684
 END IF
653 685

  
654

  
655

  
656 686
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
657 687

  
658 688
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
......
715 745
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
716 746
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
717 747

  
748
  if (AnaGeom) THEN
749
     REWIND(IOIN)
750
     READ(IOIN,AnaList,IOSTAT=IoS)
751
     IF (IOS==0) THEN ! we have a AnaList namelist 
752
        Call ReadAnaList
753
     END IF
754
     IF (AnaFile=="EMPTY") THEN
755
        AnaFile=TRIM(PathName) // '.dat'
756
     END IF
757
     OPEN(IODat,File=AnaFile)
758
  END IF
759

  
760
        
761

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

  
799
  AtName=""
800
  MassAt=1.
801

  
755 802
  ! We read the initial geometries
756 803
  Call Read_Geom(input)
757 804

  
......
761 808
     Atome(I)=ConvertNumAt(AtName(I))
762 809
  END DO
763 810

  
764
!*********** HERE *************
765
!* To DO: 
766
!* deplacer le test sur les frozen dans une sous routine
767
!* deplacer la lecture de sinput dans une sous routine
768
!* creer lecture input siesta
769
!* creer egrad siesta
770
!*
771
!*************************************
811
 If (AnaGeom) THEN
812
! We print the description of the varialbes in AnaFile
813
    Call PrintAnaList(IoDat)
814
 END IF
772 815

  
773 816
  ! PFL 23 Sep 2008 ->
774 817
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
......
1129 1172

  
1130 1173
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1131 1174
  Call Write_path(IOpt)
1175
! Shall we analyze the geometries ?
1176
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1177

  
1132 1178
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1133 1179

  
1134 1180
  ! We have the geometries, the first gradients... we will generate the first hessian matrices :-)
......
1803 1849
     GradOld=Grad
1804 1850
     EPathOld=Energies
1805 1851

  
1852

  
1853

  
1806 1854
     ! Shall we re-parameterize the path ?
1807 1855
     ! PFL 2007/Apr/10 ->
1808 1856
     ! We call PathCreate in all cases, it will take care of the 
......
1902 1950
        END DO ! matches DO IGeom=1,NGeomF
1903 1951

  
1904 1952
        Call Write_path(IOpt)
1953
! Shall we analyze the geometries ?
1954
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1905 1955

  
1906 1956
        ! We update the Hessian matrices
1907 1957
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
......
2048 2098

  
2049 2099
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2050 2100

  
2101
  Close(IOIN)
2102
  Close(IOOUT)
2103
  IF (AnaGeom) Close(IODAT)
2104

  
2051 2105
END PROGRAM PathOpt
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/Calc_baker.f90 (revision 7)
78 78
       real(KREAL) ::  angle_d,ca,sa
79 79
     END FUNCTION ANGLE_D
80 80

  
81
 
82 81

  
83 82

  
83

  
84 84
     SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
85 85
       !     
86 86
       ! This subroutine uses the description of a list of Coordinates
......
118 118
  END INTERFACE
119 119

  
120 120
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121
!
121
  !
122 122
!!!!!!!! PFL Debug
123
!
123
  !
124 124
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
125 125
  INTEGER(KINT) :: I1,I2,I3,I4
126 126
  LOGICAL :: debugPFL
......
130 130

  
131 131
  debug=valid("Calc_baker")
132 132
  debugPFL=valid("BakerPFL")
133
  if (debug) WRITE(*,*) "============= Entering Cal_baker =============="
134 133

  
134
  if (debug) Call Header(" Entering Cal_baker ")
135

  
135 136
  IF (Nat.le.2) THEN
136 137
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
137 138
     RETURN
......
164 165
  z_refGeom(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat) 
165 166

  
166 167
!!!!!!!!!!!!!!!!!!!!
167
     !
168
     !   PFL Debug
169
     !
168
  !
169
  !   PFL Debug
170
  !
170 171
!!!!!!!!!!!!!!!!!!!!
171
     if (debugPFL) THEN
172
        WRITE(*,*) "DBG PFL Calc_BAKER - Calculating Zmat"
172
  if (debugPFL) THEN
173
     WRITE(*,*) "DBG PFL Calc_BAKER - Calculating Zmat"
173 174
     if (.not.ALLOCATED(indzmat)) THEN
174 175
        ALLOCATE(indzmat(Nat,5))
175 176
     END IF
......
182 183
     IF (.NOT.ALLOCATED(indzmat0)) THEN
183 184
        ALLOCATE(indzmat0(nat,5))
184 185
     END IF
185
! We construct a Zmat
186
          IF (Frozen(1).NE.0) THEN 
187
           Renum=.TRUE.
186
     ! We construct a Zmat
187
     IF (Frozen(1).NE.0) THEN 
188
        Renum=.TRUE.
188 189
!!! remplcaer indzmat par indzmat0 puis renumerotation
189
           call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat, &
190
        call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat, &
190 191
             x_refGeom,y_refGeom,z_refGeom, &
191
                r_cov,fact,frozen)
192
        ELSE
193
           !no zmat... we create our own 
194
           call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat, &
192
             r_cov,fact,frozen)
193
     ELSE
194
        !no zmat... we create our own 
195
        call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat, &
195 196
             x_refGeom,y_refGeom,z_refGeom, &
196
                r_cov,fact)
197
        END IF
197
             r_cov,fact)
198
     END IF
198 199

  
199
         DO I=1,Nat
200
           Order(IndZmat0(I,1))=I
201
           OrderInv(I)=IndZmat0(I,1)
200
     DO I=1,Nat
201
        Order(IndZmat0(I,1))=I
202
        OrderInv(I)=IndZmat0(I,1)
203
     END DO
204
     IndZmat(1,1)=Order(IndZmat0(1,1))
205
     IndZmat(2,1)=Order(IndZmat0(2,1))
206
     IndZmat(2,2)=Order(IndZmat0(2,2))
207
     IndZmat(3,1)=Order(IndZmat0(3,1))
208
     IndZmat(3,2)=Order(IndZmat0(3,2))
209
     IndZmat(3,3)=Order(IndZmat0(3,3))
210
     DO I=4,Nat
211
        DO J=1,4
212
           IndZmat(I,J)=Order(IndZmat0(I,J))
202 213
        END DO
203
        IndZmat(1,1)=Order(IndZmat0(1,1))
204
        IndZmat(2,1)=Order(IndZmat0(2,1))
205
        IndZmat(2,2)=Order(IndZmat0(2,2))
206
        IndZmat(3,1)=Order(IndZmat0(3,1))
207
        IndZmat(3,2)=Order(IndZmat0(3,2))
208
        IndZmat(3,3)=Order(IndZmat0(3,3))
209
        DO I=4,Nat
210
           DO J=1,4
211
              IndZmat(I,J)=Order(IndZmat0(I,J))
212
           END DO
213
        end do
214
     end do
214 215

  
215 216

  
216
        WRITE(*,*) "DBG PFL CalcBaker, IndZmat0"
217
        DO I=1,Nat
218
            WRITE(*,*) I, indZmat0(I,:)
219
        end do
217
     WRITE(*,*) "DBG PFL CalcBaker, IndZmat0"
218
     DO I=1,Nat
219
        WRITE(*,*) I, indZmat0(I,:)
220
     end do
220 221

  
221
        WRITE(*,*) "DBG PFL CalcBaker, IndZmat"
222
        DO I=1,Nat
223
            WRITE(*,*) I, indZmat(I,:)
224
        end do
222
     WRITE(*,*) "DBG PFL CalcBaker, IndZmat"
223
     DO I=1,Nat
224
        WRITE(*,*) I, indZmat(I,:)
225
     end do
225 226

  
226 227

  
227 228
     DO I=1,Nat
......
229 230
        I2=indzmat0(I,2)
230 231
        I3=indzmat0(I,3)
231 232
        I4=indzmat0(I,4)
232
! Adding BOND between at1 and at2
233
      WRITE(*,*)  "DBG Indzmat",I,I1,I2,I3,I4
234
         if (I2.NE.0) THEN
233
        ! Adding BOND between at1 and at2
234
        WRITE(*,*)  "DBG Indzmat",I,I1,I2,I3,I4
235
        if (I2.NE.0) THEN
235 236
           NbBonds=NbBonds+1
236
    ! values of the coordinate (bond) is calculated for the reference geometry
237
    ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.)
238
    ! are determined using all geometries. vecteur is defined in bib_oper.f
239
                  Call vecteur(I1,I2,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
240
!                  CurrentCoord%value=Norm2
241
                  WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Bond",NbBonds," Norm2,val_zmat",Norm2,val_zmat(I,1)
242
                  CurrentCoord%value=val_zmat(I,1)
243
                  CurrentCoord%SignDihedral = 0
244
                  CurrentCoord%At1=I1
245
                  CurrentCoord%At2=I2
246
                  CurrentCoord%Type="BOND"
247
                  IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
248
                     ALLOCATE(CurrentCoord%next)
249
                     NULLIFY(CurrentCoord%next%next)
250
                  END IF
251
                  CurrentCoord => CurrentCoord%next
237
           ! values of the coordinate (bond) is calculated for the reference geometry
238
           ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.)
239
           ! are determined using all geometries. vecteur is defined in bib_oper.f
240
           Call vecteur(I1,I2,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
241
           !                  CurrentCoord%value=Norm2
242
           WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Bond",NbBonds," Norm2,val_zmat",Norm2,val_zmat(I,1)
243
           CurrentCoord%value=val_zmat(I,1)
244
           CurrentCoord%SignDihedral = 0
245
           CurrentCoord%At1=I1
246
           CurrentCoord%At2=I2
247
           CurrentCoord%Type="BOND"
248
           IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
249
              ALLOCATE(CurrentCoord%next)
250
              NULLIFY(CurrentCoord%next%next)
251
           END IF
252
           CurrentCoord => CurrentCoord%next
252 253
        END IF !IE.NE.0
253
       if (I3.NE.0) THEN
254
        if (I3.NE.0) THEN
254 255

  
255
                       NbAngles=NbAngles+1
256
                       ! values of the coordinate (bond) is calculated for the reference geometry
257
                       ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
258
                       ! determined using all geometries.
259
                       Call vecteur(i2,I3,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
260
                       Call vecteur(i2,I1,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
261
!                       CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
262
                  WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Angle",NbAngles," angle,val_zmat", &
263
                       angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),val_zmat(I,2)
264
                       CurrentCoord%value=val_zmat(I,2)*Pi/180.
265
                       CurrentCoord%SignDihedral = 0
266
                       CurrentCoord%At1=I1
267
                       CurrentCoord%At2=I2
268
                       CurrentCoord%At3=I3
269
                       CurrentCoord%Type="ANGLE"
270
                       IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
271
                          ALLOCATE(CurrentCoord%next)
272
                          NULLIFY(CurrentCoord%next%next)
273
                       END IF
274
                       CurrentCoord => CurrentCoord%next
275
                    END IF ! matches IF (I3.NE.0)
276
         if (I4.NE.0) THEN
277
                       NbDihedrals=NbDihedrals+1
278
              Call vecteur(I2,I1,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
279
              Call vecteur(I2,I3,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
280
              Call vecteur(I3,I4,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
281
                       CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
282
                            vx4,vy4,vz4,norm4)
283
                       CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
284
                            vx5,vy5,vz5,norm5)
285
!                       CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
286
!                            vx2,vy2,vz2,norm2)*Pi/180.
287
                  WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Dihedral",NbDihedrals,    &
288
                       " angle_d,val_zmat", &
289
                       angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
290
                       vx2,vy2,vz2,norm2) &
291
                       ,val_zmat(I,3)
256
           NbAngles=NbAngles+1
257
           ! values of the coordinate (bond) is calculated for the reference geometry
258
           ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
259
           ! determined using all geometries.
260
           Call vecteur(i2,I3,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
261
           Call vecteur(i2,I1,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
262
           !                       CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
263
           WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Angle",NbAngles," angle,val_zmat", &
264
                angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),val_zmat(I,2)
265
           CurrentCoord%value=val_zmat(I,2)*Pi/180.
266
           CurrentCoord%SignDihedral = 0
267
           CurrentCoord%At1=I1
268
           CurrentCoord%At2=I2
269
           CurrentCoord%At3=I3
270
           CurrentCoord%Type="ANGLE"
271
           IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
272
              ALLOCATE(CurrentCoord%next)
273
              NULLIFY(CurrentCoord%next%next)
274
           END IF
275
           CurrentCoord => CurrentCoord%next
276
        END IF ! matches IF (I3.NE.0)
277
        if (I4.NE.0) THEN
278
           NbDihedrals=NbDihedrals+1
279
           Call vecteur(I2,I1,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
280
           Call vecteur(I2,I3,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
281
           Call vecteur(I3,I4,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff