Statistiques
| Révision :

root / src / Extrapol_int.f90 @ 7

Historique | Voir | Annoter | Télécharger (16,67 ko)

1
SUBROUTINE Extrapol_int(iopt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
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
  ! Default value of FRot=.TRUE.
12
  ! Default value of NMaxPtPath = 1000
13
  ! IntCoordI(:,:) ! (NGeomI,3*Nat-6)
14
  ! XyzGeomF(:,:,:) ! (NGeomF,3,Nat)
15
  ! Default value of FAlign=.TRUE.
16
  use Path_module, only : NMaxPtPath, IntCoordI, Pi, IndZmat, XyzGeomF, &
17
       IntCoordF, IntTangent, Renum, Nom, Order, MassAt, SGeom, XyzGeomI, &
18
       Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms,Align
19
  ! IndZmat(Nat,5)
20

    
21
  use Io_module
22

    
23

    
24
  IMPLICIT NONE
25

    
26

    
27
  REAL(KREAL), INTENT(OUT) :: s
28
  ! X0(Nat),Y0(Nat),Z0(Nat): reference geometry.
29
  REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat)
30
  REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)
31
  ! Iopt: Number of the cycles for the optimization
32
  INTEGER(KINT), INTENT(IN) :: Iopt
33
  REAL(KREAL), INTENT(OUT) :: XGeomF(NGeomF)
34

    
35
  INTEGER(KINT) :: IdxGeom, I, J, K, Idx
36
  REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v
37
  REAL(KREAL) :: a_val, d
38

    
39
  REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
40
  REAL(KREAL), ALLOCATABLE :: ValZmat(:,:),DerInt(:) ! (Nat,3)
41

    
42

    
43
  LOGICAL ::  debug, print,printspline
44
  LOGICAL, EXTERNAL :: valid
45

    
46
  INTEGER(KINT) :: NSpline
47
  CHARACTER(LCHARS) :: FileSpline,TmpChar
48

    
49

    
50
  !  LOGICAL :: FAlign=.FALSE., FRot=.FALSE.
51

    
52
  ! We will calculate the length of the path, in MW coordinates...
53
  ! this is done is a stupid way: we interpolate the zmatrix values,
54
  ! convert them into cartesian, weight the cartesian
55
  ! and calculate the evolution of the distance ! 
56
  ! We have to follow the same procedure for every geometry
57
  ! so even for the first one, we have to convert it from zmat
58
  ! to cartesian !
59

    
60
  debug=valid("pathcreate")
61
  print=valid("printgeom")
62
  printspline=(valid("printspline").AND.(dist<=1e30))
63

    
64
  if (debug) Call Header("Entering Extrapol_int")
65

    
66
! We want 100 points along the interpolating path
67
  NSpline=int(NMaxPtPath/100)
68
  if (printspline) THEN
69
     WRITE(TmpChar,'(I5)') Iopt
70
     FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar))
71
     OPEN(IOTMP,FILE=FileSpline)
72

    
73
  END IF
74

    
75
  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),valzmat(Nat,3),DerInt(NCoord))
76
  IdxGeom=1
77

    
78
  XYZTmp2(1,1)=0.
79
  XYZTmp2(1,2)=0.
80
  XYZTmp2(1,3)=0.
81
  XYZTmp2(2,2)=0.
82
  XYZTmp2(2,3)=0.
83
  XYZTmp2(3,3)=0.
84

    
85
  valzmat=0.
86
  valzmat(2,1)=IntCoordI(1,1)
87
  valzmat(3,1)=IntCoordI(1,2)
88
  valzmat(3,2)=IntCoordI(1,3)*180./Pi
89
  Idx=4
90
  DO I=4,Nat
91
     ValZmat(I,1)=IntCoordI(1,Idx)
92
     Idx=Idx+1
93
     DO J=2,3
94
        ValZmat(I,J)=IntCoordI(1,Idx)*180./Pi
95
        Idx=Idx+1
96
     END DO
97
  END DO
98

    
99
  XyzTmp2(2,1)=valzmat(2,1)
100
  d=valzmat(3,1)
101
  a_val=valzmat(3,2)/180.*Pi
102
  !              write(*,*) "aval,pi",a_val,valzmat(3,2),pi
103
  if (Nat.GE.3) THEN
104
     if (IndZmat(3,2).EQ.1)  THEN
105
        XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)
106
     ELSE
107
        XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)
108
     ENDIF
109
     XyzTmp2(3,2)=d*sin(a_val)
110
  ENDIF
111
  !              i=1
112
  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
113
  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)
114
  !              i=2
115
  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
116
  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)
117
  !              i=3
118
  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
119
  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)
120

    
121
  DO i=4,Nat
122
     call ConvertZmat_cart(i,IndZmat,valzmat,                &
123
          XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
124
     !                  WRITE(*,*) 'TOTOZma:',i,IndZmat(I,1),            &
125
     !                        (IndZmat(I,J+1),valzmat(I,J),J=1,3)
126
     !                   WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
127
  END DO
128

    
129
  ! We align this geometry with the original one
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
133
     ! The rotation matrix MRot has INTENT(OUT) in subroutine rotation_matrix(...),
134
     ! which is called in the CalcRmsd(...).
135
     ! PFL 24 Nov 2008 ->
136
     ! If we have frozen atoms we align only those ones.
137
     if (NFroz.GT.0) THEN
138
        Call AlignPartial(Nat,x0,y0,z0,                     &
139
             xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
140
             FrozAtoms,MRot)
141
     ELSE
142
        Call  CalcRmsd(Nat, x0,y0,z0,                              &
143
             xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
144
             MRot,rmsd,.TRUE.,.TRUE.)
145
     END IF
146
     ! <- PFL 24 Nov 2008
147
  END IF
148

    
149
  XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
150
  IntCoordF(IdxGeom,:)=IntCoordI(1,:)
151

    
152
  ! We calculate the first derivatives
153
  u=0.d0
154
  DO Idx=1,NCoord
155
     call splintder(u,v,DerInt(iDX),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
156
  END DO
157
  IntTangent(IdxGeom,:)=DerInt
158

    
159
  if (print.AND.(Dist.LE.1e20)) THEN
160
     WRITE(IOOUT,'(1X,I5)') Nat
161
     WRITE(IOOUT,*) "# Cartesian Coordinates for geom",IdxGeom
162
     DO i=1,Nat
163
        If (Renum) THEN
164
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),      &
165
                (XyzTmp2(Order(I),J),J=1,3)
166
        ELSE
167
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),      &
168
                (XyzTmp2(I,J),J=1,3)
169
        END IF
170
     END DO
171
  END IF
172

    
173
  ! We initialize the first geometry
174
  ! PFL 26 Nov 2008 ->
175
  ! Now, I copy XyzTmp2 into XyzTmp because XyzTmp2 has been aligned
176
  ! with the reference geometry
177

    
178
  XyzTmp=XyzTmp2
179

    
180
  !   XyzTmp(1,1)=0.
181
  !   XyzTmp(1,2)=0.
182
  !   XyzTmp(1,3)=0.
183
  !   XyzTmp(2,2)=0.
184
  !   XyzTmp(2,3)=0.     
185
  !   XyzTmp(3,3)=0.
186

    
187
  !   XyzTmp(2,1)=valzmat(2,1)
188
  !   d=valzmat(3,1)
189
  !   a_val=valzmat(3,2)/180.*Pi
190

    
191
  !   if (Nat.GE.3) THEN
192
  !      if (IndZmat(3,2).EQ.1)  THEN
193
  !         XyzTmp(3,1)=XYzTmp(1,1)+d*cos(a_val)
194
  !      ELSE
195
  !         XyzTmp(3,1)=XyzTmp(2,1)-d*cos(a_val)
196
  !      ENDIF
197
  !      XyzTmp(3,2)=d*sin(a_val)
198
  !   ENDIF
199

    
200
  !   DO i=4,Nat
201
  !      call ConvertZmat_cart(i,IndZmat,valzmat,        &
202
  !           XYzTMP(1,1), XYzTMP(1,2),XYzTMP(1,3))
203
  !   END DO
204

    
205
  ! <- PFL 26 Nov 2008
206

    
207
  s=0.
208
  if (printspline) THEN
209
     SELECT CASE (Nat)
210
     CASE(2)
211
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)
212
     CASE (3)
213
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)
214
     CASE(4:)
215
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)
216
     END SELECT
217
  END IF
218
  valzmat=0.d0
219

    
220
  DO K=1,NMaxPtPath
221
     u=real(K)/NMaxPtPath*(NGeomI-1.)
222

    
223
     XYZTmp2(1,1)=0.
224
     XYZTmp2(1,2)=0.
225
     XYZTmp2(1,3)=0.
226
     XYZTmp2(2,2)=0.
227
     XYZTmp2(2,3)=0.
228
     XYZTmp2(3,3)=0.
229

    
230
     ! We generate the interpolated Zmatrix
231
     call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))
232
     valzmat(2,1)=v
233
     call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))
234
     valzmat(3,1)=v
235
     call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))
236
     valzmat(3,2)=v*180./Pi
237
     Idx=4
238
     DO I=4,Nat
239
        call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
240
        valzmat(I,1)=v
241
        Idx=Idx+1
242
        DO J=2,3
243
           call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
244
           valzmat(I,J)=v*180./pi
245
           Idx=Idx+1
246
        END DO
247
     END DO ! matches DO I=4,Nat
248

    
249
     ! We convert it into Cartesian coordinates
250
     XyzTmp2(2,1)=valzmat(2,1)
251
     d=valzmat(3,1)
252
     a_val=valzmat(3,2)/180.*Pi
253
     if (Nat.GE.3) THEN
254
        if (IndZmat(3,2).EQ.1)  THEN
255
           XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)
256
        ELSE
257
           XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)
258
        ENDIF
259
        XyzTmp2(3,2)=d*sin(a_val)
260
     ENDIF
261

    
262
     DO I=4,Nat
263
        call ConvertZmat_cart(I,IndZmat,valzmat,       &
264
             XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
265
     END DO
266

    
267
     IF ((Nat.GE.4).OR.Align) THEN
268
        ! PFL 24 Nov 2008 ->
269
        ! If we have frozen atoms we align only those ones.
270
        if (NFroz.GT.0) THEN
271
           Call AlignPartial(Nat,x0,y0,z0,                     &
272
                xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
273
                FrozAtoms,MRot)
274
        ELSE
275
           ! PFL 26 Nov 2008!
276
           ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here
277
           ! we align on the previous geom... what is the best ? Is there a difference ?
278
           Call  CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &
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

    
285
     ds=0.
286
     DO I=1,Nat
287
        DO J=1,3
288
           ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2
289
           XYZTmp(I,J)=XyzTMP2(I,J)
290
        END DO
291
     END DO
292
     
293

    
294

    
295
     s=s+sqrt(ds)
296

    
297
     if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN
298
        SELECT CASE (Nat)
299
        CASE(2)
300
           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)
301
        CASE (3)
302
           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)
303
        CASE(4:)
304
           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)
305
        END SELECT
306
     END IF
307

    
308
     !         if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist
309
     if (s>=dist) THEN
310
        if (debug) THEN
311
           WRITE(*,*) "DBG Zmat",s
312
           WRITE(*,'(1X,I5)') IndZmat(1,1)
313
           WRITE(*,'(1X,I5,I5,1X,F10.4)') IndZmat(2,1),IndZmat(2,2),valzmat(2,1)
314
           WRITE(*,'(1X,I5,2(I5,1X,F10.4))') IndZmat(3,1),IndZmat(3,2),valzmat(3,1)  &
315
                ,IndZmat(3,3),valzmat(3,2)
316
           DO I=4,Nat
317
              WRITE(*,'(1X,I5,3(I5,1X,F10.4))') IndZmat(I,1),IndZmat(I,2),valzmat(I,1)  &
318
                   ,IndZmat(I,3),valzmat(I,2),IndZmat(I,4),valzmat(I,3)
319
           END DO
320
        END IF ! matches if (debug) THEN
321

    
322
        s=s-dist
323
        IdxGeom=IdxGeom+1
324
        SGeom(IdxGeom)=s+IdxGeom*dist
325
        XgeomF(IdxGeom)=u
326
        XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
327
        IntCoordF(IdxGeom,1)=valzmat(2,1)
328
        IntCoordF(IdxGeom,2)=valzmat(3,1)
329
        IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi
330
        Idx=4
331
        DO I=4,Nat
332
           IntCoordF(IdxGeom,Idx)=valzmat(I,1)
333
           IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi
334
           IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi
335
           Idx=Idx+3
336
        END DO
337
        IntTangent(IdxGeom,:)=DerInt
338

    
339
        if (print) THEN
340
           WRITE(IOOUT,'(1X,I5)') Nat
341
           WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
342
           ! PFL 17/July/2006: only if we have more than 4 atoms.
343
           IF ((Nat.GE.4).OR.Align) THEN
344
              ! PFL 24 Nov 2008 ->
345
              ! If we have frozen atoms we align only those ones.
346
              if (NFroz.GT.0) THEN
347
                 Call AlignPartial(Nat,x0,y0,z0,                     &
348
                      xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
349
                      FrozAtoms,MRot)
350
              ELSE
351
                 Call  CalcRmsd(Nat, x0,y0,z0,                               &
352
                      xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),   &
353
                      MRot,rmsd,.TRUE.,.TRUE.)
354
              END IF
355
              ! <- PFL 24 Nov 2008
356

    
357
           END IF
358

    
359
           DO I=1,Nat
360
              If (Renum) THEN
361
                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),   &
362
                      (XyzTmp2(Order(I),J),J=1,3)
363
              ELSE
364
                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),    &
365
                      (XyzTmp2(I,J),J=1,3)
366
              END IF
367
           END DO
368

    
369
        END IF ! matches if (print) THEN
370
     END IF ! matches if (s>=dist) THEN
371
  END DO ! matches DO K=1,NMaxPtPath, Is it correct??
372

    
373

    
374
  if (printspline) THEN
375
     SELECT CASE (Nat)
376
     CASE(2)
377
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)
378
     CASE (3)
379
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)
380
     CASE(4:)
381
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)
382
     END SELECT
383
  END IF
384

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

    
391

    
392
!      ! We generate the interpolated Zmatrix
393
!      call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))
394
!      valzmat(2,1)=v
395
!      call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))
396
!      valzmat(3,1)=v
397
!      call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))
398
!      valzmat(3,2)=v*180./Pi
399
!      Idx=4
400
!      DO I=4,Nat
401
!         call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
402
!         valzmat(I,1)=v
403
!         Idx=Idx+1
404
!         DO J=2,3
405
!            call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
406
!            valzmat(I,J)=v*180./pi
407
!            Idx=Idx+1
408
!         END DO
409
!      END DO ! matches DO I=4,Nat
410

    
411
!      ! We convert it into Cartesian coordinates
412
!      XyzTmp2(2,1)=valzmat(2,1)
413
!      d=valzmat(3,1)
414
!      a_val=valzmat(3,2)/180.*Pi
415
!      if (Nat.GE.3) THEN
416
!         if (IndZmat(3,2).EQ.1)  THEN
417
!            XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)
418
!         ELSE
419
!            XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)
420
!         ENDIF
421
!         XyzTmp2(3,2)=d*sin(a_val)
422
!      ENDIF
423

    
424
!      DO I=4,Nat
425
!         call ConvertZmat_cart(I,IndZmat,valzmat,       &
426
!              XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
427
!      END DO
428

    
429
!      IF (Nat.GE.4) THEN
430
!         ! PFL 24 Nov 2008 ->
431
!         ! If we have frozen atoms we align only those ones.
432
!         if (NFroz.GT.0) THEN
433
!            Call AlignPartial(Nat,x0,y0,z0,                     &
434
!                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
435
!                 FrozAtoms,MRot)
436
!         ELSE
437
!            ! PFL 26 Nov 2008!
438
!            ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here
439
!            ! we align on the previous geom... what is the best ? Is there a difference ?
440
!            Call  CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &
441
!                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
442
!                 MRot,rmsd,.TRUE.,.TRUE.)
443
!         END IF
444
!         ! <- PFL 24 Nov 2008
445
!      END IF
446

    
447
     IdxGeom=IdxGeom+1
448

    
449

    
450

    
451
     IF (IdxGeom.GT.NGeomF) THEN
452
        WRITE(IOOUT,*) "!!! ERROR in Extrapol_int !!!!"
453
        WRITE(IOOUT,*) "Too many structures. Increase NMaxPath"
454
        WRITE(*,*) "** PathCreate ***"
455
        WRITE(*,*) "Distribution of points along the path is wrong."
456
        WRITE(*,*) "Increase value of NMaxPtPath in the input file"
457
        WRITE(*,*) "Present value is:", NMaxPtPath
458
        STOP
459
     END IF
460

    
461
     SGeom(IdxGeom)=s+IdxGeom*dist
462
     XgeomF(IdxGeom)=min(u,NGeomI-1.d0)
463
     XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
464
     !  XyzGeomF(IdxGeom,:,:)=XyzTmp2(:,:)
465
     IntCoordF(IdxGeom,1)=valzmat(2,1)
466
     IntCoordF(IdxGeom,2)=valzmat(3,1)
467
     IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi
468
     Idx=4
469
     DO I=4,Nat
470
        IntCoordF(IdxGeom,Idx)=valzmat(I,1)
471
        IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi
472
        IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi
473
        Idx=Idx+3
474
     END DO
475
     IntTangent(IdxGeom,:)=DerInt     
476

    
477
     if (print) THEN
478
        WRITE(IOOUT,'(1X,I5)') Nat
479
        WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
480
        ! PFL 17/July/2006: only if we have more than 4 atoms.
481
        IF ((Nat.GE.4).OR.Align) THEN
482
           ! PFL 24 Nov 2008 ->
483
           ! If we have frozen atoms we align only those ones.
484
           if (NFroz.GT.0) THEN
485
              Call AlignPartial(Nat,x0,y0,z0,                     &
486
                   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
487
                   FrozAtoms,MRot)
488
           ELSE
489
              Call  CalcRmsd(Nat, x0,y0,z0,                               &
490
                   xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),   &
491
                   MRot,rmsd,.TRUE.,.TRUE.)
492
           END IF
493
           ! <- PFL 24 Nov 2008
494
        END IF
495

    
496
        DO I=1,Nat
497
           If (Renum) THEN
498
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
499
                   (XyzTmp2(Order(I),J),J=1,3)
500
           ELSE
501
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),     &
502
                   (XyzTmp2(I,J),J=1,3)
503
           END IF
504
        END DO ! matches DO I=1,Nat
505
     END IF ! matches if (print) THEN
506
  END IF ! matches if (s>=0.9*dist) THEN
507

    
508
  if (debug) WRITE(*,*) 's final =',s
509

    
510
  DEALLOCATE(XyzTmp,XyzTmp2,valzmat)
511

    
512
  if (printspline) CLOSE(IOTMP)
513

    
514
  if (debug) Call Header("Extrapol_int Over")
515

    
516

    
517
END SUBROUTINE EXTRAPOL_INT