Statistics
| Revision:

root / src / CalcDist_cart.f90 @ 7

History | View | Annotate | Download (12.3 kB)

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

    
3
! This subroutine computes the curvilinear distance of the images on the path
4
! It is adapted from Extrapol_cart 
5

    
6
  use Path_module
7
  use Io_module
8

    
9

    
10
  IMPLICIT NONE
11

    
12

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

    
27
  REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
28
  REAL(KREAL), ALLOCATABLE :: DerXyz(:,:) ! Nat,3
29

    
30

    
31
  LOGICAL ::  debug, print, printspline
32
  LOGICAL, EXTERNAL :: valid
33

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

    
42
  debug=valid("pathcreate")
43
  print=valid("printgeom")
44

    
45
  printspline=(valid("printspline").AND.(dist<=1e30))
46

    
47
  if (debug) Call Header("Entering CalcDist_cart")
48

    
49
! We want 100 points along the interpolating path
50
  NSpline=int(NMaxPtPath/100)
51

    
52
  if (printspline) THEN
53
     WRITE(TmpChar,'(I5)') Iopt
54
     FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar))
55
     OPEN(IOTMP,FILE=FileSpline)
56

    
57
  END IF
58

    
59

    
60
  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),DerXyz(Nat,3))
61

    
62
!  XyzTmp2=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/))
63

    
64
  if (debug) THEN
65
     WRITE(*,*) "DBG Extrapol_cart Initial geometries"
66
     DO IGeom=1,NGeomI
67
        WRITE(*,*) 'XyzGeomI, IGeom=',IGeom
68
        DO I=1,Nat
69
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),      &
70
                (XyzGeomI(IGeom,J,I),J=1,3)
71
        END DO
72
     END DO
73
  END IF
74

    
75

    
76
! In order to mesure only the relevant distance between two points
77
! we align all geometries on the original one
78

    
79
   DO IGeom=1,NGeomI
80

    
81
      XyzTmp2=Reshape(XyzGeomI(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/))
82
  ! We align this geometry with the original one
83
  ! PFL 17/July/2006: only if we have more than 4 atoms.
84
!  IF (Nat.GT.4) THEN
85
! PFL 24 Nov 2008 ->
86
! If we have frozen atoms we align only those ones.
87
! PFL 8 Feb 2011 ->
88
! I add a flag to see if we should align or not.
89
! For small systems, it might be better to let the user align himself
90
      IF (Align) THEN
91
         if (NFroz.GT.0) THEN
92
            Call AlignPartial(Nat,x0,y0,z0,                     &
93
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
94
                 FrozAtoms,MRot)
95
         ELSE
96
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
97
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
98
                 MRot,rmsd,.TRUE.,.TRUE.)
99
         END IF
100
! <- PFL 24 Nov 2008
101
      END IF
102
! -> PFL 8 Feb 2011
103
!  END IF
104

    
105
  XyzGeomI(IGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
106
 END DO
107

    
108
   if (print) THEN
109
      WRITE(*,*) "Aligned geometries"
110
      DO J=1, NGeomI
111
         WRITE(IOOUT,*) Nat
112
         WRITE(IOOUT,*) " Aligned geometry ",J,"/",NGeomI
113
           DO i=1,Nat
114
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
115
                   XyzGeomI(J,1:3,I)
116
           END DO
117
        END DO
118
     END IF
119

    
120
  XyzGeomF(1,:,:)=XyzGeomI(1,:,:)
121

    
122

    
123
  ! We initialize the first geometry
124

    
125
  XyzTmp=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/))
126

    
127
  ! We align this geometry with the original one
128
  ! PFL 17/July/2006: only if we have more than 4 atoms.
129
!  IF (Nat.GT.4) THEN
130
! PFL 24 Nov 2008 ->
131
! If we have frozen atoms we align only those ones.
132
! PFL 8 Feb 2011 ->
133
! I add a flag to see if we should align or not.
134
! For small systems, it might be better to let the user align himself
135
      IF (Align) THEN
136
         if (NFroz.GT.0) THEN
137
            Call AlignPartial(Nat,x0,y0,z0,                     &
138
                 xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3),       &
139
                 FrozAtoms,MRot)
140
         ELSE
141
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
142
                 xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3),       &
143
                 MRot,rmsd,.TRUE.,.TRUE.)
144
         END IF
145
! <- PFL 24 Nov 2008
146
      END IF
147
! -> PFL 8 Feb 2011
148
!  END IF
149

    
150

    
151

    
152
  s=0.d0
153
  SGeom(1)=0.d0
154

    
155
  if (printspline) THEN
156
     u=0.d0
157
     DO Iat=1,Nat
158
        ! We generate the interpolated coordinates
159
        if (Linear) THEN
160
           call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat))
161
           call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat))
162
           call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat))
163

    
164
        ELSE
165
           call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2))
166
           call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1))
167
           call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat))
168
        END IF
169
     END DO
170

    
171
     WRITE(IOTMP,*) Nat
172
     WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s
173
     DO Iat=1,Nat
174
        WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)),    &
175
                   (XyzTmp2(Iat,1:3))
176
     END DO
177
  END IF
178

    
179
 IdxGeom=1
180

    
181
  DO K=1,NMaxPtPath
182
     u=real(K)/NMaxPtPath*(NGeomI-1.)
183

    
184
     DO Iat=1,Nat
185
        ! We generate the interpolated coordinates
186
        if (Linear) THEN
187
           call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat))
188
           call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat))
189
           call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat))
190

    
191
        ELSE
192
           call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2))
193
           call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1))
194
           call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat))
195
        END IF
196
     END DO
197

    
198

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

    
222

    
223

    
224
     ds=0.
225
     DO I=1,Nat
226
        DO J=1,3
227
           ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2
228
           XYZTmp(I,J)=XyzTMP2(I,J)
229
        ENDDO
230
     ENDDO
231
     s=s+sqrt(ds)
232

    
233
     !         if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist
234

    
235
  if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN
236
     WRITE(IOTMP,*) Nat
237
     WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s
238
     DO Iat=1,Nat
239
        WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)),    &
240
                   (XyzTmp2(Iat,1:3))
241
     END DO
242
  END IF
243

    
244

    
245
     if (s>=dist) THEN
246

    
247
        if (debug) THEN
248
           WRITE(*,*) "DBG Interpol_cart",s
249
           DO i=1,Nat
250
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
251
                   (XyzTmp2(I,J),J=1,3)
252
           END DO
253
        END IF
254

    
255
        s=s-dist
256
        IdxGeom=IdxGeom+1
257
        SGeom(IdxGeom)=s+dist*(IdxGeom-1)
258

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

    
282
        XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
283
        XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/))
284

    
285
        if (print) THEN
286
           WRITE(IOOUT,'(1X,I5)') Nat
287
           WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
288

    
289
           DO i=1,Nat
290
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
291
                   (XyzTmp2(I,J),J=1,3)
292
           END DO
293

    
294
        END IF
295
     END IF  ! s>= dist
296
  ENDDO  ! K
297

    
298

    
299
  if (s>=0.9*dist) THEN
300
     s=s-dist
301
     XyzTmp2=Reshape(XyzGeomI(NGeomI,:,:),(/Nat,3/),ORDER=(/2,1/))
302

    
303
  if (printspline) THEN
304
     WRITE(IOTMP,*) Nat
305
     WRITE(IOTMP,*) "Debug Spline - Coord=Cart,s=",s
306
     DO Iat=1,Nat
307
        WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)),    &
308
                   (XyzTmp2(Iat,1:3))
309
     END DO
310
  END IF
311

    
312

    
313
  ! We align this geometry with the original one
314
  ! PFL 17/July/2006: only if we have more than 4 atoms.
315
!  IF (Nat.GT.4) THEN
316
! PFL 24 Nov 2008 ->
317
! If we have frozen atoms we align only those ones.
318
! PFL 8 Feb 2011 ->
319
! I add a flag to see if we should align or not.
320
! For small systems, it might be better to let the user align himself
321
      IF (Align) THEN
322
         if (NFroz.GT.0) THEN
323
            Call AlignPartial(Nat,x0,y0,z0,                     &
324
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
325
                 FrozAtoms,MRot)
326
         ELSE
327
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
328
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
329
                 MRot,rmsd,.TRUE.,.TRUE.)
330
         END IF
331
! <- PFL 24 Nov 2008
332
      END IF
333
! -> PFL 8 Feb 2011
334
!  END IF
335

    
336
     IdxGeom=IdxGeom+1
337
     SGeom(IdxGeom)=s+dist*(IdxGeom-1)
338

    
339
    IF (IdxGeom.GT.NGeomF) THEN
340
        WRITE(IOOUT,*) "!!! ERROR in Extrapol_cart !!!!"
341
        WRITE(IOOUT,*) "Too many structures. Increase NMaxPath"
342
        WRITE(*,*) "** PathCreate ***"
343
        WRITE(*,*) "Distribution of points along the path is wrong."
344
        WRITE(*,*) "Increase value of NMaxPtPath in the input file"
345
        WRITE(*,*) "Present value is:", NMaxPtPath
346
        STOP
347
     END IF
348

    
349

    
350
  ! We align this geometry with the original one
351
  ! PFL 17/July/2006: only if we have more than 4 atoms.
352
!  IF (Nat.GT.4) THEN
353
! PFL 24 Nov 2008 ->
354
! If we have frozen atoms we align only those ones.
355
! PFL 8 Feb 2011 ->
356
! I add a flag to see if we should align or not.
357
! For small systems, it might be better to let the user align himself
358
      IF (Align) THEN
359
         if (NFroz.GT.0) THEN
360
            Call AlignPartial(Nat,x0,y0,z0,                     &
361
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
362
                 FrozAtoms,MRot)
363
         ELSE
364
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
365
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
366
                 MRot,rmsd,.TRUE.,.TRUE.)
367
         END IF
368
! <- PFL 24 Nov 2008
369
      END IF
370
! -> PFL 8 Feb 2011
371
!  END IF
372

    
373
     XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
374
     XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/))
375

    
376
     if (print) THEN
377
        WRITE(IOOUT,'(1X,I5)') Nat
378
        WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
379

    
380
        DO i=1,Nat
381
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),     &
382
                (XyzTmp2(I,J),J=1,3)
383
        END DO
384
     END IF
385
  END IF
386

    
387
  DEALLOCATE(XyzTmp,XyzTmp2)
388

    
389
  if (debug) WRITE(*,*) 's final =',s
390

    
391
  if (printspline) CLOSE(IOTMP)
392

    
393
  if (debug) Call Header("Extrapol_cart over")
394

    
395
END SUBROUTINE EXTRAPOL_CART