Statistics
| Revision:

root / src / Extrapol_cart.f90 @ 5

History | View | Annotate | Download (10.7 kB)

1
SUBROUTINE Extrapol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef)
2

    
3
  ! This subroutine constructs the path, and if dist<>Infinity, it samples
4
  ! the path to obtain geometries.
5
  ! Basically, you call it twice: i) dist=infinity, it will calculate the length of the path
6
  ! ii) dist finite, it will give you the images you want along the path.
7
  !
8
  ! For now, it gives equidistant geometries
9
  !
10

    
11
  use Path_module
12
  use Io_module
13

    
14

    
15
  IMPLICIT NONE
16

    
17

    
18
  REAL(KREAL), INTENT(OUT) :: s
19
  ! X0(Nat),Y0(Nat),Z0(Nat): reference geometry.
20
  REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat)
21
  REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)
22
! Iopt: Number of the cycles for the optimization
23
  INTEGER(KINT), INTENT(IN) :: Iopt
24
  INTEGER(KINT) :: IdxGeom, I, J, K, Iat
25
! NSpline is the number of points along the interpolating path
26
  INTEGER(KINT) :: NSpline
27
! FileSpline: Filename to save the interpolating path coordinates
28
  CHARACTER(LCHARS) :: FileSpline,TmpChar
29
  REAL(KREAL) :: Rmsd, MRot(3, 3), ds, u
30

    
31
  REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
32
  REAL(KREAL), ALLOCATABLE :: DerXyz(:,:) ! Nat,3
33

    
34

    
35
  LOGICAL ::  debug, print, printspline
36
  LOGICAL, EXTERNAL :: valid
37

    
38
  !We will calculate the length of the path, in MW coordinates...
39
  ! this is done is a stupid way: we interpolate the zmatrix values,
40
  ! convert them into cartesian, weight the cartesian
41
  ! and calculate the evolution of the distance ! 
42
  ! We have to follow the same procedure for every geometry
43
  ! so even for the first one, we have to convert it from zmat
44
  ! to cartesian !
45

    
46
  debug=valid("pathcreate")
47
  print=valid("printgeom")
48

    
49
  printspline=(valid("printspline").AND.(dist<=1e30))
50

    
51
  if (debug) Call Header("Entering Extrapol_cart")
52

    
53
! We want 100 points along the interpolating path
54
  NSpline=int(NMaxPtPath/100)
55

    
56
  if (printspline) THEN
57
     WRITE(TmpChar,'(I5)') Iopt
58
     FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar))
59
     OPEN(IOTMP,FILE=FileSpline)
60

    
61
  END IF
62

    
63

    
64
  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),DerXyz(Nat,3))
65

    
66
! PFL 2011 June 22
67
! There was an alignment procedure there.
68
! It has been moved to PathCreate before the 
69
! computation of the Spline coefficient.
70

    
71
     if (debug) THEN
72
        WRITE(*,*) "Extrapol_cart- L72 - geometries "
73
        DO I=1,NGeomI
74
           WRITE(*,*) "Geom  ",I
75
           DO J=1,Nat
76
              If (renum) THEN
77
                 Iat=Order(J)
78
                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,Iat)
79
              ELSE
80
                 Iat=OrderInv(J)
81
                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),XyzGeomI(I,1:3,J)
82
              END IF
83
           END DO
84
        END DO
85
     END IF
86

    
87

    
88
  XyzGeomF(1,:,:)=XyzGeomI(1,:,:)
89

    
90

    
91
  ! We initialize the first geometry
92

    
93
  XyzTmp=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/))
94

    
95

    
96

    
97
  s=0.d0
98
  SGeom(1)=0.d0
99

    
100
  if (printspline) THEN
101
     u=0.d0
102
     DO Iat=1,Nat
103
        ! We generate the interpolated coordinates
104
        if (Linear) THEN
105
           call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat))
106
           call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat))
107
           call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat))
108

    
109
        ELSE
110
           call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2))
111
           call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1))
112
           call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat))
113
        END IF
114
     END DO
115

    
116
     WRITE(IOTMP,*) Nat
117
     WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s
118
     DO Iat=1,Nat
119
        WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)),    &
120
                   (XyzTmp2(Iat,1:3))
121
     END DO
122
  END IF
123

    
124
 IdxGeom=1
125

    
126
  DO K=1,NMaxPtPath
127
     u=real(K)/NMaxPtPath*(NGeomI-1.)
128

    
129
     DO Iat=1,Nat
130
        ! We generate the interpolated coordinates
131
        if (Linear) THEN
132
           call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat))
133
           call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat))
134
           call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat))
135

    
136
        ELSE
137
           call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2))
138
           call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1))
139
           call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat))
140
        END IF
141
     END DO
142

    
143

    
144
! PFL 2011 June 22: We now align on the previous one
145
! as we do for internal coord interpolation
146

    
147
! We align this geometry with the previous one
148

    
149
  ! PFL 17/July/2006: only if we have more than 4 atoms.
150
!  IF (Nat.GT.4) THEN
151
! PFL 24 Nov 2008 ->
152
! If we have frozen atoms we align only those ones.
153
! PFL 8 Feb 2011 ->
154
! I add a flag to see if we should align or not.
155
! For small systems, it might be better to let the user align himself
156
      IF (Align) THEN
157
         if (NFroz.GT.0) THEN
158
            Call AlignPartial(Nat,                     &
159
                 xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3),       &
160
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
161
                 FrozAtoms,MRot)
162
         ELSE
163
            Call  CalcRmsd(Nat,                              &
164
                 xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3),       &
165
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
166
                 MRot,rmsd,.TRUE.,.TRUE.)
167
         END IF
168
! <- PFL 24 Nov 2008
169
      END IF
170
! -> PFL 8 Feb 2011
171
!  END IF
172

    
173

    
174

    
175
     ds=0.
176
     DO I=1,Nat
177
        DO J=1,3
178
           ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2
179
           XYZTmp(I,J)=XyzTMP2(I,J)
180
        ENDDO
181
     ENDDO
182
     s=s+sqrt(ds)
183

    
184
     !         if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist
185

    
186
  if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN
187
     WRITE(IOTMP,*) Nat
188
     WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s
189
     DO Iat=1,Nat
190
        WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)),    &
191
                   (XyzTmp2(Iat,1:3))
192
     END DO
193
  END IF
194

    
195

    
196
     if (s>=dist) THEN
197

    
198
        if (debug) THEN
199
           WRITE(*,*) "DBG Interpol_cart",s
200
           DO i=1,Nat
201
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
202
                   (XyzTmp2(I,J),J=1,3)
203
           END DO
204
        END IF
205

    
206
        s=s-dist
207
        IdxGeom=IdxGeom+1
208
        SGeom(IdxGeom)=s+dist*(IdxGeom-1)
209

    
210
  ! We align this geometry with the original one
211
  ! PFL 17/July/2006: only if we have more than 4 atoms.
212
!  IF (Nat.GT.4) THEN
213
! PFL 24 Nov 2008 ->
214
! If we have frozen atoms we align only those ones.
215
! PFL 8 Feb 2011 ->
216
! I add a flag to see if we should align or not.
217
! For small systems, it might be better to let the user align himself
218
      IF (Align) THEN
219
         if (NFroz.GT.0) THEN
220
            Call AlignPartial(Nat,x0,y0,z0,                     &
221
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
222
                 FrozAtoms,MRot)
223
         ELSE
224
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
225
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
226
                 MRot,rmsd,.TRUE.,.TRUE.)
227
         END IF
228
! <- PFL 24 Nov 2008
229
      END IF
230
! -> PFL 8 Feb 2011
231
!  END IF
232

    
233
        XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
234
        XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/))
235

    
236
        if (print) THEN
237
           WRITE(IOOUT,'(1X,I5)') Nat
238
           WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
239

    
240
           DO i=1,Nat
241
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
242
                   (XyzTmp2(I,J),J=1,3)
243
           END DO
244

    
245
        END IF
246
     END IF  ! s>= dist
247
  ENDDO  ! K
248

    
249

    
250
  if (s>=0.9*dist) THEN
251
     s=s-dist
252
     XyzTmp2=Reshape(XyzGeomI(NGeomI,:,:),(/Nat,3/),ORDER=(/2,1/))
253

    
254
  if (printspline) THEN
255
     WRITE(IOTMP,*) Nat
256
     WRITE(IOTMP,*) "Debug Spline - Coord=Cart,s=",s
257
     DO Iat=1,Nat
258
        WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)),    &
259
                   (XyzTmp2(Iat,1:3))
260
     END DO
261
  END IF
262

    
263

    
264
  ! We align this geometry with the original one
265
  ! PFL 17/July/2006: only if we have more than 4 atoms.
266
!  IF (Nat.GT.4) THEN
267
! PFL 24 Nov 2008 ->
268
! If we have frozen atoms we align only those ones.
269
! PFL 8 Feb 2011 ->
270
! I add a flag to see if we should align or not.
271
! For small systems, it might be better to let the user align himself
272
      IF (Align) THEN
273
         if (NFroz.GT.0) THEN
274
            Call AlignPartial(Nat,x0,y0,z0,                     &
275
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
276
                 FrozAtoms,MRot)
277
         ELSE
278
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
279
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
280
                 MRot,rmsd,.TRUE.,.TRUE.)
281
         END IF
282
! <- PFL 24 Nov 2008
283
      END IF
284
! -> PFL 8 Feb 2011
285
!  END IF
286

    
287
     IdxGeom=IdxGeom+1
288
     SGeom(IdxGeom)=s+dist*(IdxGeom-1)
289

    
290
    IF (IdxGeom.GT.NGeomF) THEN
291
        WRITE(IOOUT,*) "!!! ERROR in Extrapol_cart !!!!"
292
        WRITE(IOOUT,*) "Too many structures. Increase NMaxPath"
293
        WRITE(*,*) "** PathCreate ***"
294
        WRITE(*,*) "Distribution of points along the path is wrong."
295
        WRITE(*,*) "Increase value of NMaxPtPath in the input file"
296
        WRITE(*,*) "Present value is:", NMaxPtPath
297
        STOP
298
     END IF
299

    
300

    
301
  ! We align this geometry with the original one
302
  ! PFL 17/July/2006: only if we have more than 4 atoms.
303
!  IF (Nat.GT.4) THEN
304
! PFL 24 Nov 2008 ->
305
! If we have frozen atoms we align only those ones.
306
! PFL 8 Feb 2011 ->
307
! I add a flag to see if we should align or not.
308
! For small systems, it might be better to let the user align himself
309
      IF (Align) THEN
310
         if (NFroz.GT.0) THEN
311
            Call AlignPartial(Nat,x0,y0,z0,                     &
312
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
313
                 FrozAtoms,MRot)
314
         ELSE
315
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
316
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
317
                 MRot,rmsd,.TRUE.,.TRUE.)
318
         END IF
319
! <- PFL 24 Nov 2008
320
      END IF
321
! -> PFL 8 Feb 2011
322
!  END IF
323

    
324
     XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
325
     XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/))
326

    
327
     if (print) THEN
328
        WRITE(IOOUT,'(1X,I5)') Nat
329
        WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
330

    
331
        DO i=1,Nat
332
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),     &
333
                (XyzTmp2(I,J),J=1,3)
334
        END DO
335
     END IF
336
  END IF
337

    
338
  DEALLOCATE(XyzTmp,XyzTmp2)
339

    
340
  if (debug) WRITE(*,*) 's final =',s
341

    
342
  if (printspline) CLOSE(IOTMP)
343

    
344
  if (debug) Call Header("Extrapol_cart over")
345

    
346
END SUBROUTINE EXTRAPOL_CART