Statistics
| Revision:

root / src / PathCreate.f90 @ 5

History | View | Annotate | Download (39.5 kB)

1
!================================================================
2
!
3
! This subroutine create a path using starting geometries
4
! it is 'original' in the sense that the path is generated
5
! by interpolating internal coordinates instead of cartesian.
6
! Comes from the Extrapol.f subroutine of Cart suite.
7
! Rewritten in F90 to be more flexible
8
!
9
!================================================================
10

    
11
SUBROUTINE PathCreate(IOpt)
12

    
13
  use Io_module
14
  use Path_module, only : Nat, NGeomI, NCoord, NGeomF, IGeomRef, NMaxL, IReparam, IReparamT, &
15
       Coord, Frozen, Cart, NCart, NFroz, XyzGeomI, atome, r_cov, fact, &
16
       IndZmat, Renum, Order, OrderInv, IntCoordI, IntCoordF,Pi,  Nom, &
17
       ISpline, IntTangent,  NPrim,Xprimitive_t, OptGeom,  &
18
       UMatF, UMat_local, XyzTangent, Linear, Align, FrozAtoms,AtName
19
  ! BMat_BakerT (3*Nat,NCoord), XyzGeomI(NGeomI,3,Nat), IGeomRef=-1 (default value)
20
  ! IntCoordI(NGeomI,NCoord)
21

    
22
  IMPLICIT NONE
23

    
24
  REAL(KREAL), PARAMETER :: Inf=1e32
25
  INTEGER(KINT), INTENT(IN) :: IOpt
26

    
27
  REAL(KREAL), ALLOCATABLE :: val_zmat(:,:), val_zmatTmp(:,:) ! (3,Nat)
28
  REAL(KREAL), ALLOCATABLE :: X0(:), Y0(:), Z0(:) ! Nat
29
  REAL(KREAL), ALLOCATABLE ::  XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
30
  !  REAL(KREAL), ALLOCATABLE ::  XyzTmp3(:,:) ! (3,Nat)
31
  REAL(KREAL), ALLOCATABLE ::  IntCoordTmp(:) ! (NCoord)
32
  INTEGER(KINT), ALLOCATABLE :: IndZmat0(:,:) !Nat,5
33
  REAL(KREAL), ALLOCATABLE :: Coef(:,:) ! (NGeomI,NCoord)
34
  ! NPrim=Number of primitive coordinates.
35
  REAL(KREAL), ALLOCATABLE :: XGeom(:) ! NGeomI
36
  REAL(KREAL), ALLOCATABLE :: XGeomF(:) ! NGeomF
37
  INTEGER(KINT) :: NGeomS
38

    
39
  REAL(KREAL) :: Rmsd, MRot(3, 3), Delta
40
  REAL(KREAL), PARAMETER ::  TwoPi=360.d0
41
  REAL(KREAL)  ::s, dist, a_val, d
42
  INTEGER(KINT) :: I, J, igeom, iat
43
  INTEGER(KINT) :: Itmp, Jat
44
  INTEGER(KINT) :: Idx,IBeg
45
  CHARACTER(LCHARS) :: Title
46
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
47
  LOGICAL ::  debug, print
48
  LOGICAL, SAVE :: First=.TRUE.
49
  ! Frozen contains the indices of frozen atoms
50
  LOGICAL, ALLOCATABLE :: FCart(:) ! Nat
51
  INTEGER, ALLOCATABLE :: AtCart(:) !Nat
52
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:),LiaisonsIni(:,:) ! (Nat,0:NMaxL)
53
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !nat
54
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !nat,nat
55
  INTEGER(KINT) :: NbFrag,NFragRef
56

    
57

    
58
  INTERFACE
59
     function valid(string) result (isValid)
60
       CHARACTER(*), intent(in) :: string
61
       logical                  :: isValid
62
     END function VALID
63

    
64
     FUNCTION test_zmat(na,ind_zmat)
65

    
66
       USE Path_module
67

    
68
       logical :: test_zmat
69
       integer(KINT) :: na
70
       integer(KINT) :: ind_zmat(Na,5)
71
     END FUNCTION TEST_ZMAT
72

    
73

    
74
    SUBROUTINE die(routine, msg, file, line, unit)
75

    
76
      Use VarTypes
77
      Use io_module
78

    
79
      implicit none
80

    
81
      character(len=*), intent(in)           :: routine, msg
82
      character(len=*), intent(in), optional :: file
83
      integer(KINT), intent(in), optional      :: line, unit
84

    
85
    END SUBROUTINE die
86

    
87
    SUBROUTINE Warning(routine, msg, file, line, unit)
88

    
89
      Use VarTypes
90
      Use io_module
91

    
92
      implicit none
93

    
94
      character(len=*), intent(in)           :: routine, msg
95
      character(len=*), intent(in), optional :: file
96
      integer(KINT), intent(in), optional      :: line, unit
97

    
98
    END SUBROUTINE Warning
99

    
100
  END INTERFACE
101

    
102

    
103

    
104
  debug=valid("pathcreate")
105
  print=valid("printgeom")
106
  if (debug) Call header("Entering PathCreate")
107
  if (debug.AND.First) WRITE(*,*) "=     First call          ="
108

    
109
  if (debug) WRITE(*,*) "Iopt,IReparam=",Iopt,IReparam
110

    
111
  ALLOCATE(Coef(NGeomI,NCoord),val_zmat(Nat,3),val_zmatTMP(Nat,3))
112
  ALLOCATE(X0(Nat),Y0(Nat),Z0(Nat),Indzmat0(Nat,5))
113
  ALLOCATE(X(Nat),Y(Nat),Z(Nat),Xgeom(NGeomI),XgeomF(NGeomF))
114
  ! ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),XyzTmp3(3,Nat))
115
  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3))
116
  ALLOCATE(IntCoordTmp(NCoord))
117

    
118
  Do I=1,NGeomI
119
     XGeom(I)=FLoat(I)-1.d0
120
     if (Print) THEN
121
        WRITE(*,*) "PathCreate - L121 - Initial geometries "
122
        DO J=1,Nat
123
           WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,J)
124
        END DO
125
     END IF
126
  END DO
127

    
128
  ! First iteration of the optimization:
129
  IF (First) THEN ! matches at L484
130
     if (debug) WRITE(*,*) "first iteration in PathCreate.90, L93"
131

    
132
     ! This is the first call, so we might have to generate a Zmat
133
     First=.FALSE.
134

    
135
     IF ( ((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD.EQ."MIXED") &
136
          .OR.(COORD.EQ."BAKER")).AND.(IGeomRef.LE.0)) THEN
137

    
138
        ! User has not defined the Reference geometry, we have to choose it ourselves!
139
        ! This is done by counting the number of fragments of each geometry. The 
140
        ! reference geometry is the one with the least fragments. When there are 
141
        ! frozen or cart atoms, they are discarded to count the fragments.
142

    
143
        ALLOCATE(Liaisons(Nat,0:NMaxL))
144
        ALLOCATE(Fragment(nat),NbAtFrag(nat),FragAt(nat,nat),FCart(Nat))
145
        FCart=.FALSE.
146
        ! FCart(Nat) was under IF below earlier but is being accessed after this loop.
147
        IF ((Frozen(1).NE.0).OR.(Cart(1).NE.0)) THEN
148
           ALLOCATE(AtCart(Nat),LiaisonsIni(Nat,0:NMaxL))
149
           FCart=.FALSE.
150
           NCart=0
151
           Idx=1
152
           DO WHILE (Frozen(Idx).NE.0)
153
              IF (Frozen(Idx).LT.0) THEN
154
                 DO I=Frozen(Idx-1),abs(Frozen(Idx))
155
                    IF (.NOT.FCart(I)) THEN
156
                       NCart=NCart+1
157
                       AtCart(NCart)=I
158
                       FCart(I)=.TRUE.
159
                    END IF
160
                 END DO
161
              ELSE IF (.NOT.FCart(Frozen(Idx))) THEN
162
                 FCart(Frozen(Idx))=.TRUE.
163
                 NCart=NCart+1
164
                 AtCart(NCart)=Frozen(Idx)
165
              END IF
166
              Idx=Idx+1
167
           END DO ! matches DO WHILE (Frozen(Idx).NE.0)
168
           NFroz=NCart
169
           Idx=1
170
           DO WHILE (Cart(Idx).NE.0)
171
              IF (Cart(Idx).LT.0) THEN
172
                 DO I=Cart(Idx-1),abs(Cart(Idx))
173
                    IF (.NOT.FCart(I)) THEN
174
                       NCart=NCart+1
175
                       AtCart(NCart)=I
176
                       FCart(I)=.TRUE.
177
                    END IF
178
                 END DO
179
              ELSEIF (.NOT.FCart(Cart(Idx))) THEN
180
                 FCart(Cart(Idx))=.TRUE.
181
                 NCart=NCart+1
182
                 AtCart(NCart)=Cart(Idx)
183
              END IF
184
              Idx=Idx+1
185
           END DO !matches DO WHILE (Cart(Idx).NE.0)
186
           IF (debug) THEN
187
              WRITE(*,'(1X,A,I4,A,I4,A)') "DBG PathCreate - GeomRef :: Found ", &
188
                   NFroz," frozen atoms, and a total of ",NCart, &
189
                   " atoms described in cartesian coordinates"
190
              WRITE(*,*) "AtCart:",AtCart(1:NCart)
191
           END IF
192
        END IF ! IF ((Frozen(1).NE.0).OR.(Cart(1).NE.0)) THEN
193

    
194
        NFragRef=2*Nat
195
        DO IGeom=1,NGeomI
196
           x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
197
           y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
198
           z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
199
           ! we calculate the connectivities
200
           Call CalcCnct(nat,atome,x,y,z,LIAISONS,r_cov,fact)
201
           ! if needed, we get rid of links between cart or frozen atoms
202
           IF (NCart.GE.2) THEN
203
              LiaisonsIni=Liaisons
204
              DO I=1,NCart
205
                 Iat=AtCart(I)
206
                 if (debug) WRITE(*,'(20(1X,I3))') I,Iat,  &
207
                      (LiaisonsIni(Iat,Itmp),ITmp=0,LiaisonsIni(Iat,0))
208
                 DO J=1,LiaisonsIni(Iat,0)
209
                    Jat=LiaisonsIni(Iat,J)
210
                    IF (FCart(Jat)) THEN
211
                       Call Annul(Liaisons,Iat,Jat)
212
                       Call Annul(Liaisons,Jat,Iat)
213
                       Call Annul(LiaisonsIni,Jat,Iat)
214
                       Liaisons(Iat,0)=Liaisons(Iat,0)-1
215
                       Liaisons(Jat,0)=Liaisons(Jat,0)-1
216
                       LiaisonsIni(Jat,0)=LiaisonsIni(Jat,0)-1
217
                    END IF
218
                 END DO
219
              END DO
220
           END IF ! matches IF (NCart.GE.2) THEN
221
           ! we now calculate the number of fragments.
222
           Call Decomp_frag(nat,liaisons,.NOT.Fcart,nbfrag,Fragment,NbAtFrag,FragAt)
223
           IF (debug) THEN
224
              WRITE(*,*) 'Debug PathCreat Line190'
225
              WRITE(*,*) 'NbFrag, NbFragRef=',NbFrag,NFragRef
226
           END IF
227

    
228
           IF (NbFrag.LT.NFragRef) THEN
229
              NFragRef=NbFrag
230
              ! The reference geometry, IGeomRef, is determined based on the least number
231
              ! of fragments (here if (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD
232
              !.EQ."MIXED")).AND.(IGeomRef.LE.0))), otherwise it still has the value of 
233
              ! IGeomRef=-1.
234
              IGeomRef=IGeom
235
           END IF
236
        END DO ! matches DO IGeom=1,NGeomI
237
        DEALLOCATE(FCart)
238
     END IF ! matches IF (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD.EQ."MIXED")
239
     ! .OR.(COORD.EQ."BAKER")).AND.(IGeomRef.LE.0)).
240

    
241
     IF ((COORD.EQ."CART").AND.(IGeomREF.LE.0)) THEN
242
        IGeomRef=1
243
        CALL Warning('PathCreate L209','IGeomRef<=0',UNIT=IOOUT)
244
     END IF
245

    
246
     if (debug) WRITE(*,*) "DBG PathCreate : IGeomRef= ",IGeomRef
247

    
248
     ! we now compute the internal coordinates for this geometry !
249
     ! The reference geometry, IGeomRef, is determined based on the least number
250
     ! of fragments if (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD
251
     !.EQ."MIXED")).AND.(IGeomRef.LE.0)), otherwise it has the value of 
252
     ! IGeomRef=-1.
253
     x(1:Nat) = XyzGeomI(IGeomRef,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat)
254
     y(1:Nat) = XyzGeomI(IGeomRef,2,1:Nat)
255
     z(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat)
256

    
257
     IndZmat=0
258
     IndZmat0=0
259

    
260
     SELECT CASE(COORD)
261
     CASE ('BAKER')
262
        ! atome(Nat), r_cov(0:Max_Z)
263
        ! Frozen(Nat) contains the indices of frozen atoms
264
        Call Calc_baker(atome,r_cov,fact,frozen,IGeomRef)
265
        if (debug) THEN
266
           WRITE(*,*) "UMat_local after Calc_baker"
267
           DO I=1,3*Nat-6
268
              WRITE(*,'(50(1X,F12.8))') UMat_local(:,I)
269
           END DO
270
        END IF
271
        DO IGeom=1,NGeomF
272
           UMatF(IGeom,:,:)=UMat_Local(:,:)
273
        END DO
274

    
275
        ALLOCATE(Xprimitive_t(NPrim))
276
        ! x, y, z corresponds the reference geometry.
277
        DO I=1,Nat
278
           Order(I)=I
279
           OrderInv(I)=I
280
        END DO
281
        x0=x
282
        y0=y
283
        z0=z
284
     CASE ('ZMAT','HYBRID')
285
        ! IN case we are using Hybrid or zmat coordinate system, we have to 
286
        ! generate the internal coordinates with a Z-Matrix
287
        ! IN case we have frozen atoms, we generate a constrained Zmat
288
        IF (Frozen(1).NE.0) THEN 
289
           Renum=.True.
290
           call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact,frozen)
291
        ELSE
292
           !no zmat... we create our own
293
           call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact)
294
           if (valid('old_zmat')) THEN 
295
              call Calc_zmat(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact)
296
              WRITE(*,*) "DBG PathCreate Calc_Zmat.... STOP"
297
              STOP
298
           END IF
299
        END IF
300

    
301
        if (debug) WRITE(*,*) 'Back from Calc_zmat'
302
        IF (.NOT.test_zmat(nat,IndZmat0)) THEN
303
           WRITE(IOOUT,*) "I cannot generate a valid Zmat... help"
304
           STOP
305
        END IF
306

    
307
        ! The IndZmat0 generated by Calc_zmat refers to the order of the input x,y,z
308
        ! However, we need it with the new order so that we can use the zmat 
309
        ! coordinates to generate cartesian coordinates.
310
        ! So, we have to reorder it...
311
        DO I=1,Nat
312
           Order(IndZmat0(I,1))=I
313
           OrderInv(I)=IndZmat0(I,1)
314
           ! I let X0 being the Cartesian coordinates of the reference geometry. Before it was
315
           ! those of the First geometry (which was also the reference one). This might change!
316
           X0(i)=X(IndZmat0(i,1))
317
           Y0(i)=Y(IndZmat0(i,1))
318
           Z0(i)=Z(IndZmat0(i,1))
319
        END DO
320
        IndZmat(1,1)=Order(IndZmat0(1,1))
321
        IndZmat(2,1)=Order(IndZmat0(2,1))
322
        IndZmat(2,2)=Order(IndZmat0(2,2))
323
        IndZmat(3,1)=Order(IndZmat0(3,1))
324
        IndZmat(3,2)=Order(IndZmat0(3,2))
325
        IndZmat(3,3)=Order(IndZmat0(3,3))
326
        DO I=4,Nat
327
           DO J=1,4
328
              IndZmat(I,J)=Order(IndZmat0(I,J))
329
           END DO
330
        end do
331
        IF (DEBUG) THEN
332
           WRITE(IOOUT,*) "DBG PathCreate : IndZmat"
333
           DO I=1,Nat
334
              WRITE(IOOUT,*) (IndZmat(I,J),J=1,5)
335
           END DO
336
           WRITE(IOOUT,*) "DBG PathCreate : IndZmat0"
337
           DO I=1,Nat
338
              WRITE(IOOUT,*) (IndZmat0(I,J),J=1,5)
339
           END DO
340

    
341
        END IF
342

    
343
        ! It is the first call, we have to create all internal coordinates
344
        if (debug) WrITE(*,*) "DBG PathCreate L302: First Call, we generate internal coords"
345
        DO igeom=1,NgeomI
346
           x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
347
           y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
348
           z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
349

    
350
           if (debug) WRITE(*,*) "DBG PathCreate L308: we generate internal coords for geom",IGeom
351

    
352
           Call ConvXyzZmat(Nat,x,y,z,IndZmat0,val_zmatTmp)
353
           IntCoordI(IGeom,1)=val_zmatTmp(2,1)
354
           IntCoordI(IGeom,2)=val_zmatTmp(3,1)
355
           IntCoordI(IGeom,3)=val_zmatTmp(3,2)*Pi/180.
356
           Idx=4
357
           DO i=4,nat
358
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,1)
359
              Idx=Idx+1
360
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180.
361
              Idx=Idx+1
362

    
363
              ! Some tests to check that the dihedral angles are similar... this is
364
              ! important if they are close to Pi.
365
              Delta=0.
366

    
367
              if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN
368
                 if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN
369
                    Delta=-1.0d0*TwoPi
370
                 ELSE
371
                    Delta=TwoPi
372
                 END IF
373
                 if (debug) THEN
374
                    WRITE(*,*) Nom(atome(i)),        &
375
                         abs(val_zmatTMP(i,3)-val_zmat(i,3)), &
376
                         val_zmatTMP(i,3),val_zmat(i,3),Delta, &
377
                         val_zmatTMP(i,3)-Delta
378
                 END IF
379
              END IF
380
              IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi
381
              Idx=Idx+1
382
           END DO
383
        END DO  ! End Loop on Igeom
384

    
385
     CASE ('MIXED')
386
        ! IN case we are using mixed coordinate system, we have to 
387
        ! we generate the internal coordinates with a Z-Matrix
388
        Renum=.TRUE.
389
        ! 
390
        call Calc_mixed_frag(Nat,atome,IndZmat0,val_zmat,x,y,z, &
391
             r_cov,fact,frozen,cart,ncart)
392

    
393
        if (debug) WRITE(*,*) 'Back from Calc_Mixed'
394

    
395
        ! How to test this partial zmat ? it must be good as it was automatically generated...
396
        !        IF (.NOT.test_zmat(nat,IndZmat0)) THEN
397
        !           WRITE(IOOUT,*) "I cannot generate a valid Zmat... help"
398
        !           STOP
399
        !        END IF
400

    
401
        ! The IndZmat0 generated by Calc_zmat refers to the order of the input x,y,z
402
        ! However, we need it with the new order so that we can use the zmat 
403
        ! coordinates to generate cartesian coordinates.
404
        ! So, we have to reorder it...
405
        DO I=1,Nat
406
           Order(IndZmat0(I,1))=I
407
           OrderInv(I)=IndZmat0(I,1)
408
           ! I let X0 being the Cartesian coordinates of the reference geometry. Before it was
409
           ! those of the First geometry (which was also the reference one). This might change!
410
           X0(i)=X(IndZmat0(i,1))
411
           Y0(i)=Y(IndZmat0(i,1))
412
           Z0(i)=Z(IndZmat0(i,1))
413
        END DO
414
        IndZmat=0
415
        IndZmat(1:NCart,:)=-1
416

    
417
        DO I=1,Nat
418
           IndZmat(I,1)=Order(IndZmat0(I,1))
419
        END DO
420

    
421
        SELECT CASE (NCart)
422
        CASE (1)
423
           IndZmat(2,2)=Order(IndZmat0(2,2))
424
           IndZmat(3,2)=Order(IndZmat0(3,2))
425
           IndZmat(3,3)=Order(IndZmat0(3,3))
426
           IdX=4
427
        CASE (2)
428
           IndZmat(3,2)=Order(IndZmat0(3,2))
429
           IndZmat(3,3)=Order(IndZmat0(3,3))
430
           Idx=4
431
        CASE DEFAULT 
432
           Idx=NCart+1
433
        END SELECT
434

    
435
        DO I=Idx,Nat
436
           DO J=1,4
437
              IndZmat(I,J)=Order(IndZmat0(I,J))
438
           END DO
439
        end do
440
        IF (DEBUG) THEN
441
           WRITE(IOOUT,*) "DBG PathCreate : IndZmat - Mixed"
442
           DO I=1,Nat
443
              WRITE(IOOUT,*) (IndZmat(I,J),J=1,5)
444
           END DO
445
        END IF
446

    
447
        ! It is the first call, we have to create all internal coordinates
448
        !            Idx=1
449
        !            DO I=1,NCart
450
        !               IntCoordI(1,Idx)=val_zmat(I,1)
451
        !               IntCoordI(1,Idx+1)=val_zmat(I,2)
452
        !               IntCoordI(1,Idx+2)=val_zmat(I,3)
453
        !               Idx=Idx+3
454
        !            END DO
455
        !            SELECT CASE (NCART)
456
        !            CASE (1)
457
        !               IntCoordI(1,Idx)=val_zmat(2,1)
458
        !               IntCoordI(1,Idx+1)=val_zmat(3,1)           
459
        !               IntCoordI(1,Idx+2)=val_zmat(3,2)*Pi/180.
460
        !               Idx=Idx+3
461
        !               IBeg=4
462
        !            CASE (2)
463
        !               IntCoordI(1,Idx)=val_zmat(3,1)           
464
        !               IntCoordI(1,Idx+1)=val_zmat(3,2)*Pi/180.
465
        !               Idx=Idx+2
466
        !               IBeg=4
467
        !            CASE DEFAULT
468
        !               IBeg=NCart+1
469
        !            END SELECT
470
        !            DO i=iBeg,Nat
471
        !               IntCoordI(1,Idx)=val_zmat(i,1)
472
        !               Idx=Idx+1
473
        !               DO j=2,3
474
        !                  IntCoordI(1,Idx)=val_zmat(i,j)*Pi/180.
475
        !                  Idx=Idx+1
476
        !               END DO
477
        !            END DO
478

    
479
        !            If (debug) WRITE(*,*) "Idx, NCoord",Idx-1,NCoord
480

    
481
        DO Igeom=1,NgeomI
482
           x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
483
           y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
484
           z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
485

    
486
           Call ConvXyzMixed(Nat,Ncart,x,y,z,IndZmat0,val_zmatTmp(1,1))
487
           Idx=1
488
           DO I=1,NCart
489
              IntCoordI(IGeom,Idx)=val_zmatTmp(I,1)
490
              IntCoordI(IGeom,Idx+1)=val_zmatTmp(I,2)
491
              IntCoordI(IGeom,Idx+2)=val_zmatTmp(I,3)
492
              Idx=Idx+3
493
           END DO
494

    
495
           SELECT CASE (NCART)
496
           CASE (1)
497
              IntCoordI(IGeom,Idx)=val_zmatTmp(2,1)
498
              Idx=Idx+1
499
              IntCoordI(IGeom,Idx)=val_zmatTmp(3,1)           
500
              IntCoordI(IGeom,Idx+1)=val_zmatTmp(3,2)*Pi/180.
501
              Idx=Idx+2
502
              IBeg=4
503
           CASE (2)
504
              IntCoordI(IGeom,Idx)=val_zmatTmp(3,1)           
505
              IntCoordI(IGeom,Idx+1)=val_zmatTmp(3,2)*Pi/180.
506
              Idx=Idx+2
507
              IBeg=4
508
           CASE DEFAULT
509
              Ibeg=Ncart+1
510
           END SELECT
511
           DO i=IBeg,Nat
512
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,1)
513
              Idx=Idx+1
514
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180.
515
              Idx=Idx+1
516

    
517
              ! Some tests to check that the dihedral angles are similar... this is
518
              ! important if they are close to Pi.
519
              Delta=0.
520

    
521
              if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN
522
                 if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN
523
                    Delta=-1.0d0*TwoPi
524
                 ELSE
525
                    Delta=TwoPi
526
                 END IF
527
                 if (debug) THEN
528
                    WRITE(*,*) Nom(atome(i)), &
529
                         abs(val_zmatTMP(i,3)-val_zmat(i,3)), &
530
                         val_zmatTMP(i,3),val_zmat(i,3),Delta, &
531
                         val_zmatTMP(i,3)-Delta
532
                 END IF
533
              END IF
534
              IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi
535
              Idx=Idx+1
536
           END DO
537

    
538
        END DO
539
     CASE ('CART') 
540
        DO I=1,Nat
541
           Order(I)=I
542
           OrderInv(I)=I
543
        END DO
544
        X0=X
545
        Y0=y
546
        Z0=z
547
     END SELECT
548

    
549
  ELSE ! First iteration of the optimization. Match at L616
550

    
551
     IGeom=1
552
     x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
553
     y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
554
     z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
555

    
556
     SELECT CASE (COORD)
557
     CASE ('HYBRID')
558
        IndZmat0=IndZmat
559
        Call ConvXyzZmat(Nat,x,y,z,IndZmat,val_zmat(1,1))
560
     CASE ('ZMAT')
561
        IndZmat0=IndZmat
562
        val_zmat=0.d0
563
        val_zmat(2,1)=IntCoordI(1,1)
564
        val_zmat(3,1)=IntCoordI(1,2)
565
        val_zmat(3,2)=IntCoordI(1,3)*180./Pi
566

    
567
        Idx=4
568
        DO iat=4,Nat
569
           val_zmat(iat,1)=IntCoordI(1,idx)
570
           Idx=Idx+1
571
           do j=2,3
572
              val_zmat(iat,J)=IntCoordI(1,idx)*180./Pi
573
              Idx=Idx+1
574
           END DO
575
        END DO
576

    
577
     CASE ('MIXED')
578
        IndZmat0=IndZmat
579
        val_zmat=0.d0
580
        Idx=1
581
        DO I=1,NCart
582
           DO J=1,3
583
              val_zmat(i,j)= IntCoordI(1,Idx)
584
              Idx=Idx+1
585
           END DO
586
        END DO
587
        SELECT CASE (NCart)
588
        CASE(1)
589
           val_zmat(2,1)=IntCoordI(1,Idx)
590
           val_zmat(3,1)=IntCoordI(1,Idx+1)
591
           val_zmat(3,2)=IntCoordI(1,Idx+2)*180./Pi
592
           Idx=Idx+3
593
           IBeg=4
594
        CASE(2)
595
           val_zmat(3,1)=IntCoordI(1,Idx)
596
           val_zmat(3,2)=IntCoordI(1,Idx)*180./Pi
597
           Idx=Idx+2
598
           IBeg=4
599
        CASE DEFAULT
600
           IBeg=NCart+1
601
        END SELECT
602

    
603
        DO Iat=IBeg,Nat
604
           val_zmat(iat,1)=IntCoordI(1,idx)
605
           Idx=Idx+1
606
           do j=2,3
607
              val_zmat(iat,J)=IntCoordI(1,idx)*180./Pi
608
              Idx=Idx+1
609
           END DO
610
        END DO
611
     CASE ('BAKER')
612
        ! Nothing to be done here.
613
     CASE ('CART')
614
        ! Nothing to be done here.
615
     CASE DEFAULT
616
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate.L554.STOP"
617
        STOP
618
     END SELECT
619
     X0=x
620
     y0=y
621
     z0=z
622
  END IF  ! First
623

    
624
     if (Print) THEN
625
        WRITE(*,*) "PathCreate - L631 - geometries "
626
        DO I=1,NGeomI
627
           DO J=1,Nat
628
              If (renum) THEN
629
                 Iat=Order(J)
630
                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,Iat)
631
              ELSE
632
                 Iat=OrderInv(J)
633
                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),XyzGeomI(I,1:3,J)
634
              END IF
635
           END DO
636
        END DO
637
     END IF
638

    
639

    
640
  ! Now that we have a zmat, we will generate all the IntCoodI corresponding...
641
  ! First one
642
  IF (COORD.EQ.'HYBRID') THEN     ! Matches at L680
643
     DO IGeom=1,NGeomI
644
        x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
645
        y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
646
        z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
647

    
648
        Call ConvXyzZmat(Nat,x,y,z,IndZmat0,val_zmatTmp(1,1))
649
        IntCoordI(IGeom,1)=val_zmatTmp(2,1)
650
        IntCoordI(IGeom,2)=val_zmatTmp(3,1)
651
        IntCoordI(IGeom,3)=val_zmatTmp(3,2)*Pi/180.
652
        Idx=4
653
        DO i=4,nat
654
           IntCoordI(IGeom,Idx)=val_zmatTmp(i,1)
655
           Idx=Idx+1
656
           IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180.
657
           Idx=Idx+1
658

    
659
           ! Some tests to check that the dihedral angles are similar... this is
660
           ! important if they are close to Pi.
661
           Delta=0.
662

    
663
           if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN
664
              if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN
665
                 Delta=-1.0d0*TwoPi
666
              ELSE
667
                 Delta=TwoPi
668
              END IF
669
              if (debug) THEN
670
                 WRITE(*,*) Nom(atome(i)), &
671
                      abs(val_zmatTMP(i,3)-val_zmat(i,3)), &
672
                      val_zmatTMP(i,3),val_zmat(i,3),Delta, &
673
                      val_zmatTMP(i,3)-Delta
674
              END IF
675
           END IF
676
           IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi
677
           Idx=Idx+1
678
        END DO
679
     END DO
680
  END IF ! Matches Coord=Hybrid
681

    
682
  ! PFL 24 Nov 2008 ->
683
  ! This should not be necessary...
684
  !   ! if we have frozen atoms, we make sure that they do not move between geometries
685
  !   IF (IntFroz.NE.0) THEN
686
  !      if (debug) WRITE(*,*) 'DBG PathCreate L641: Number of frozen coordinates=', IntFroz
687
  !      DO IGeom=2,NGeomI ! For IOpt .GT. 0, NGeomI is equal to NGeomF
688
  !         IntCoordI(IGeom,1:IntFroz)=IntCoordI(1,1:IntFroz)
689
  !      END DO
690
  !   END IF
691
  ! -> PFL 24 Nov 2008
692

    
693
  Idx=min(9,NCoord)
694

    
695
  IF (DEBUG.AND.(COORD.NE.'CART')) THEN
696
     WRITE(*,*) "DBG PathCreate. L685 IntCoordI(IGeom,:)"
697
     DO I=1,NGeomI
698
        WRITE(*,'("Geom:",I5,9(1X,F10.4))') I,IntCoordI(I,1:Idx)
699
     END DO
700
     ! We write it also in terms of Zmat
701
     IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'HYBRID')) THEN
702
        DO I=1,NGeomI
703
           WRITE(*,*) 'Zmat for geom ',I,'/',NGeomI,'Xgeom=',Xgeom(I)
704
           WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(1,1)))
705
           WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(2,1))), &
706
                IndZmat(2,2),IntCoordI(I,1)
707
           WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(3,1))), &
708
                IndZmat(3,2),IntCoordI(I,2),IndZmat(3,3),IntCoordI(I,3)*180./Pi
709
           Idx=4
710
           DO J=4,Nat
711
              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(J,1))), &
712
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3),          &
713
                   IntCoordI(I,Idx+1)*180./Pi, &
714
                   IndZmat(J,4),IntCoordI(I,Idx+2)*180./Pi
715
              Idx=Idx+3
716
           END DO
717

    
718
           XyzTmp2=0.
719
           ! We convert it into Cartesian coordinates
720
           XyzTmp2(2,1)=IntCoordI(I,1)
721
           d=IntCoordI(I,2)
722
           a_val=IntCoordI(I,3)
723
           if (Nat.GE.3) THEN
724
              if (IndZmat(3,2).EQ.1)  THEN
725
                 XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)
726
              ELSE
727
                 XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)
728
              ENDIF
729
              XyzTmp2(3,2)=d*sin(a_val)
730
           ENDIF
731

    
732
           Idx=4
733
           DO iat=4,Nat
734
              val_zmatTMP(iat,1)=IntCoordI(I,idx)
735
              Idx=Idx+1
736
              DO j=2,3
737
                 val_zmatTMP(iat,J)=IntCoordI(I,idx)*180./Pi
738
                 Idx=Idx+1
739
              END DO
740
           END DO
741
           DO iat=4,Nat
742
              call ConvertZmat_cart(iat,IndZmat,val_zmatTMP, &
743
                   XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
744
           END DO
745
           DO Iat=1,nat
746
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(Iat,1))),&
747
                   XyzTmp2(iat,1:3)
748
           END DO
749

    
750
        END DO
751
     END IF
752
     IF (COORD.EQ.'MIXED') THEN
753
        DO I=1,NGeomI
754
           WRITE(*,*) 'Cart+Zmat for geom ',I,'/',NGeomI,'Xgeom=',Xgeom(I)
755
           Idx=1
756
           XyzTmp2=0.
757
           DO J=1,NCART
758
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(1,1))), &
759
                   IntCoordI(I,Idx:Idx+2)
760
              XyzTmp2(J,1:3)=IntCoordI(I,Idx:Idx+2)
761
              Idx=Idx+3
762
           END DO
763
           SELECT CASE (NCart)
764
           CASE (1)
765
              J=2
766
              WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
767
                   IndZmat(J,2),IntCoordI(I,Idx)
768
              Idx=Idx+1
769
              J=3
770
              WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
771
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), &
772
                   IntCoordI(I,Idx+1)*180./Pi
773
              Idx=Idx+2
774
              Ibeg=4
775
           CASE (2)
776
              J=3
777
              WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
778
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), &
779
                   IntCoordI(I,Idx+1)*180./Pi
780
              Idx=Idx+2
781
              Ibeg=4
782
           CASE DEFAULT
783
              IBeg=NCart+1
784
           END SELECT
785
           DO J=IBeg,Nat
786
              IF (J.EQ.2) &
787
                   WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
788
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), &
789
                   IntCoordI(I,Idx+1)*180./Pi, &
790
                   IndZmat(J,4),IntCoordI(I,Idx+2)*180./Pi
791
              Idx=Idx+3
792
           END DO
793

    
794
           ! We convert it into Cartesian coordinates
795
           IntCoordTmp=IntCoordI(I,1:NCoord)
796
           Call Mixed2cart(Nat,IndZmat,IntCoordTmp,XyzTmp2)
797

    
798
           DO Iat=1,nat
799
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(Iat,1))),&
800
                   XyzTmp2(iat,1:3)
801
           END DO
802

    
803
        END DO
804
     END IF
805
  END IF ! matches IF (DEBUG.AND.(COORD.NE.'CART')) THEN
806

    
807
  if (debug) WRITE(*,*) "Before interpolation, PathCreate.f90, L740, IOpt=",IOpt, &
808
       "ISpline=", ISpline 
809

    
810
     if (Print) THEN
811
        WRITE(*,*) "PathCreate - L811 - geometries "
812
        DO I=1,NGeomI
813
           DO J=1,Nat
814
              If (renum) THEN
815
                 Iat=Order(J)
816
                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,Iat)
817
              ELSE
818
                 Iat=OrderInv(J)
819
                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),XyzGeomI(I,1:3,J)
820
              END IF
821
           END DO
822
        END DO
823
     END IF
824

    
825

    
826
  ! Now comes the Interpolation:
827
  IF ((NGeomI>2).AND.(IOpt.GE.ISpline)) THEN
828
     Linear=.FALSE.
829
     ! We have at least three geometries so we can apply cubic spline
830
     SELECT CASE (COORD)
831
     CASE ('ZMAT','HYBRID','MIXED')
832
        DO I=1,NCoord
833
           ! We compute the spline coefficients
834
           call spline(xGeom,IntCoordI(1,I),NGeomI,Inf,Inf, Coef(1,I))
835
           if (debug) THEN
836
              WRITE(*,*) 'Coef Spline for coord',I
837
              WRITE(*,*) Coef(1:NGeomI,I)
838
           END IF
839
        END DO
840
     CASE ('BAKER')
841
        ! We compute the spline coefficients
842
        ! xGeom(NGeomI), IntCoordI(NGeomI,3*Nat-6) declared in Path_module.f90
843
        ! can also be used for the Baker's internal coordinates.
844
        ! we need to get Baker's internal coordinates from cartesian coordinates.
845
        ! as we have subroutine ConvXyzmat(...) in Z matrix case.
846
        ! Coef(NGeomI,NCoord) or (NGeomI,N3at) has INTENT(OUT)
847
        DO I=1,NCoord ! NCoord=3*Nat-6-symmetry_eliminated.
848
           call spline(xGeom,IntCoordI(1,I),NGeomI,Inf,Inf,Coef(1,I))
849
           ! IF (debug) THEN
850
           !   WRITE(*,*) 'Coef Spline for coord',I
851
           !  WRITE(*,*) Coef(1:NGeomI,I)
852
           !    END IF
853
        END DO
854
     CASE ('CART')
855

    
856
!! PFL 2011 June 22 ->
857
!! The alignment procedure from Extrapol_cart has been moved
858
!! here, before computation of the spline coefficients
859

    
860
!  XyzTmp2=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/))
861

    
862
  if (debug) THEN
863
     WRITE(*,*) "DBG PathCreate - Calc Spline coeff CART- Initial geometries"
864
     DO IGeom=1,NGeomI
865
        WRITE(*,*) 'XyzGeomI, IGeom=',IGeom
866
        DO I=1,Nat
867
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),      &
868
                (XyzGeomI(IGeom,J,I),J=1,3)
869
        END DO
870
     END DO
871
  END IF
872

    
873

    
874
! In order to mesure only the relevant distance between two points
875
! we align all geometries on the original one
876

    
877
   DO IGeom=1,NGeomI
878

    
879
      XyzTmp2=Reshape(XyzGeomI(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/))
880
  ! We align this geometry with the original one
881
  ! PFL 17/July/2006: only if we have more than 4 atoms.
882
!  IF (Nat.GT.4) THEN
883
! PFL 24 Nov 2008 ->
884
! If we have frozen atoms we align only those ones.
885
! PFL 8 Feb 2011 ->
886
! I add a flag to see if we should align or not.
887
! For small systems, it might be better to let the user align himself
888
      IF (Align) THEN
889
         if (NFroz.GT.0) THEN
890
            Call AlignPartial(Nat,x0,y0,z0,                     &
891
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
892
                 FrozAtoms,MRot)
893
         ELSE
894
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
895
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
896
                 MRot,rmsd,.TRUE.,.TRUE.)
897
         END IF
898
! <- PFL 24 Nov 2008
899
      END IF
900
! -> PFL 8 Feb 2011
901
!  END IF
902

    
903
  XyzGeomI(IGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
904
 END DO
905

    
906
   if (debug) THEN
907
     WRITE(*,*) "DBG PathCreate - Calc Spline coeff CART- Aligned geometries"
908
      DO J=1, NGeomI
909
         WRITE(IOOUT,*) Nat
910
         WRITE(IOOUT,*) " Aligned geometry ",J,"/",NGeomI
911
           DO i=1,Nat
912
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
913
                   XyzGeomI(J,1:3,I)
914
           END DO
915
        END DO
916
     END IF
917

    
918
!! <- PFL 2011 June 22
919

    
920
        DO I=1,Nat
921
           ! We compute the spline coefficients
922
           call spline(xGeom,XyzGeomI(1,1,I),NGeomI,Inf,Inf, Coef(1,3*I-2))
923
           call spline(xGeom,XyzGeomI(1,2,I),NGeomI,Inf,Inf, Coef(1,3*I-1))
924
           call spline(xGeom,XyzGeomI(1,3,I),NGeomI,Inf,Inf, Coef(1,3*I))
925
           if (debug) THEN
926
              WRITE(*,*) 'Coef Spline for coord',i
927
              WRITE(*,*) Coef(1:NGeomI,I)
928
           END IF
929
        END DO
930
     CASE DEFAULT
931
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate.L783.STOP"
932
        STOP
933
     END SELECT
934
  ELSE ! matches IF ((NGeomI>2).AND.(Iopt.Ge.ISpline)) THEN
935
     ! On n'a que deux images. On fera une interpolation lineaire,
936
     ! qui est un cas particulier avec Coef=0
937
     if (debug.AND.(NGeomI.EQ.2)) WRITE(*,*) "DBG : Only 2 starting geometries"
938
     if (debug.AND.(Iopt.LE.ISpline)) WRITE(*,*) "DBG : Not enough path -> linear int"
939
     Coef=0.d0
940
     Linear=.TRUE.
941
  END IF ! matches IF ((NGeomI>2).AND.(Iopt.Ge.ISpline)) THEN
942
  ! Shall we redistribute the points along the path ?
943
  if (debug) WRITE(*,*) "PathCreate.f90, L795, IOpt=", IOpt, "IReparam=", IReparam
944

    
945
  IF (MOD(IOpt,IReparam).EQ.0) THEN
946
     ! We generate the Path, in two steps:
947
     ! i) we calculate the length of the mass weighted (MW) path
948
     ! ii) we sample the path to get the geometries
949

    
950
     NGeomS=NGeomI
951
     SELECT CASE (COORD)
952
     CASE ('ZMAT','HYBRID')
953
        ! We call the extrapolation once to get the path length
954
        dist=Inf
955
        Call ExtraPol_int(IOpt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
956
        ! we have the length, we calculate the distance between two points
957

    
958
        if (s.LE.1e-6) THEN
959
           IF (OptGeom.LE.0) THEN
960
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
961
              WRITE(*,*) "ERROR !!! Stop !"
962
              STOP
963
           END IF
964
        ELSE
965
           ! we have the length, we calculate the distance between two points
966
           dist=s/Real(NGeomf-1)
967
           if (debug) WRITE(*,*) 'Dbg PathCreate ZMAT s, dist:',s,dist
968
           ! we call it again to obtain the geometries !
969
           Call ExtraPol_int(iopt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
970
        END IF
971
     CASE ('BAKER')
972
        ! We generate the Path, in two steps:
973
        ! i) we calculate the length of the mass weighted (MW) path
974
        ! ii) we sample the path to get the geometries
975

    
976
        ! We call the extrapolation first time to get the path length
977
        ! iopt: Number of the cycles for the optimization.
978
        ! s: a real number and has INTENT(OUT), probably returns the 
979
        ! length of the path.
980
        ! dist=Inf=> ExtraPol_baker has to calculate the length of the path.
981
        ! X0(Nat),Y0(Nat),Z0(Nat): reference geometry.
982
        ! Xgeom: Xgeom(NGeomI) has INTENT(IN)
983
        ! Coef(NGeomI,NCoord) is obtained from spline interpolation subroutine.
984
        ! XgeomF: XGeomF(NGeomF) has INTENT(OUT)
985
        ! Which reference coordinate geometry X0, Y0, Z0 to use???
986
        dist=Inf
987
        !if (debug) WRITE(*,*) "Before the call of ExtraPol_baker in redistribution &
988
        !of the points, in PathCreate.f90, L836"
989
        !Print *, 'Before ExtraPol_baker, PathCreate.90, L848, IntCoordI='
990
        !Print *, IntCoordI
991
        Call ExtraPol_baker(s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
992
        !Print *, 'After ExtraPol_baker, PathCreate.90, L848, IntCoordI='
993
        !Print *, IntCoordI
994
        !if (debug) WRITE(*,*) "After the call of ExtraPol_baker in the redistribution &
995
        ! of the points, in PathCreate.f90, L843"
996

    
997
        ! we have the lenght, we calculate the distance between two points
998
        !Print *, 'After first call of ExtraPol_baker in PathCreate.f90, s=', s
999

    
1000
        if (s.LE.1e-6) THEN
1001
           IF (OptGeom.LE.0) THEN
1002
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
1003
              WRITE(*,*) "ERROR !!! Stop !"
1004
              STOP
1005
           END IF
1006
        ELSE
1007

    
1008
           dist=s/Real(NGeomf-1)
1009
           !if (debug) WRITE(*,*) 'Dbg PathCreate ZMAT s, dist:',s,dist
1010

    
1011
           ! we call it again to obtain the geometries !
1012
           Call ExtraPol_baker(s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
1013
           WRITE(*,*) "At the end of Baker case in the redistribution of the points, in PathCreate.f90, L966"
1014
        END IF
1015

    
1016
        DO IGeom=1,NGeomF
1017
           WRITE(*,*) "UMatF for IGeom=",IGeom
1018
           DO I=1,3*Nat-6
1019
              WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1020
           END DO
1021
        END DO
1022

    
1023
     CASE ('MIXED')
1024
        !           WRITE(*,*) "IOIOIOIOIOIOIOI"
1025
        ! We call the extrapolation once to get the path length
1026
        Dist=Inf
1027
        Call ExtraPol_Mixed(s,dist,x0,y0,z0,xgeom,Coef)
1028

    
1029
        if (s.LE.1e-6) THEN
1030
           IF (OptGeom.LE.0) THEN
1031
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
1032
              WRITE(*,*) "ERROR !!! Stop !"
1033
              STOP
1034
           END IF
1035
        ELSE
1036

    
1037
           ! we have the length, we calculate the distance between two points
1038
           dist=s/Real(NGeomf-1)
1039
           if (debug) WRITE(*,*) 'Dbg PathCreate Mixed s, dist:',s,dist
1040

    
1041
           ! we call it again to obtain the geometries !
1042
           Call ExtraPol_Mixed(s,dist,x0,y0,z0,xgeom,Coef)
1043

    
1044
           IF (debug) THEN
1045
              DO I=1,NGeomF
1046
                 WRITE(*,'(1X,A,12(1X,F10.5))') "IntTangent:",IntTangent(I,1:NCoord)
1047
              END DO
1048
           END IF
1049
        END IF
1050
     CASE ('CART')
1051
        ! We call the extrapolation once to get the path length
1052
        Dist=Inf
1053
        Call ExtraPol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef)
1054

    
1055
        if (s.LE.1e-6) THEN
1056
           IF (OptGeom.LE.0) THEN
1057
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
1058
              WRITE(*,*) "ERROR !!! Stop !"
1059
              STOP
1060
           END IF
1061
        ELSE
1062
           ! we have the lenght, we calculate the distance between two points
1063
           dist=s/Real(NGeomf-1)
1064
           if (debug) WRITE(*,*) 'Debug s, dist:',s,dist
1065

    
1066
           ! we call it again to obtain the geometries !
1067
           Call ExtraPol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef)
1068
        END IF
1069
     END SELECT
1070
  ELSE
1071
     ! PFL 3 Jan 2008 -> 
1072
     ! Tangents are not recalculated if the points are not reparameterized
1073
     ! PFL 6 Apr 2008 -> 
1074
     ! Unless user has asked to reparameterized them !
1075
     IF (mod(Iopt,iReparamT).EQ.0) THEN
1076
        XGeomF=XGeom
1077
        NGeomS=NGeomI
1078
        Call CalcTangent(XgeomF,NGeomS,xgeom,Coef)
1079
        !
1080
     END IF
1081
     ! <- PFL 3 Jan 2008 / 6 Apr 2008
1082
  END IF ! if (mod(Iop,iReparam).EQ.0)
1083
  if (debug) WRITE(*,*) "PathCreate.f90, L885"
1084

    
1085
  ! PFL 22.Nov.2007 ->
1086
  ! We do not write tangents for baker coordinates because they are hard to 
1087
  ! interpret for a chemist.
1088

    
1089
  IF (debug.AND.(COORD.NE."BAKER")) THEN
1090
     ! we write the tangents :-)
1091
     DO I=1,NGeomF
1092
        SELECT CASE(Coord)
1093
        CASE ('ZMAT','MIXED')
1094
           IntCoordTmp=IntTangent(I,:)
1095
           WRITE(*,'(1X,A,12(1X,F10.5))') "IntCoordTmp:",IntCoordTmp(1:NCoord)
1096
           WRITE(*,'(1X,A,12(1X,F10.5))') "IntTangent:",IntTangent(I,1:NCoord)
1097
        CASE ('CART','HYBRID')
1098
           IntCoordTmp=XyzTangent(I,:)
1099
        CASE DEFAULT
1100
           WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate. STOP"
1101
           STOP
1102
        END SELECT
1103
        WRITE(Title,"('Tangents for geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
1104
        Call PrintGeom(Title,Nat,Ncoord,IntCoordTmp,Coord,6,Atome,Order,OrderInv,IndZmat) 
1105
     END DO
1106
  END IF
1107

    
1108

    
1109
  IF (DEBUG.AND.(COORD.NE.'CART')) THEN
1110
     Idx=min(12,NCoord)
1111
     WRITE(*,*) "DBG PathCreate. IntCoordF(IGeom,:)"
1112
     DO I=1,NGeomF
1113
        WRITE(*,'("Geom:",I5,12(1X,F10.4))') I,IntCoordF(I,1:Idx)
1114
     END DO
1115
  END IF
1116

    
1117
  DEALLOCATE(Coef,val_zmat,val_zmatTMP)
1118
  DEALLOCATE(X0,Y0,Z0,Indzmat0)
1119
  DEALLOCATE(X,Y,Z,Xgeom,XgeomF)
1120
  DEALLOCATE(XyzTmp,XyzTmp2)
1121

    
1122
  if (debug) Call header("PathCreate over")
1123

    
1124
END SUBROUTINE PathCreate
1125