Statistiques
| Révision :

root / src / Calc_baker.f90 @ 7

Historique | Voir | Annoter | Télécharger (40,79 ko)

1
SUBROUTINE Calc_baker(atome,r_cov,fact,frozen,IGeomRef)
2
  !     
3
  ! This subroutine analyses a geometry to construct the baker
4
  ! delocalized internal coordinates
5
  ! v1.0
6
  ! We use only one geometry
7
  !     
8
  Use Path_module, only : max_Z,NMaxL,Nom,MaxFroz,NFroz,Pi,Massat,a0,BMat_BakerT, &
9
       Nat,NCoord,XyzGeomI,NGeomI,NGeomF,UMat,NPrim,BTransInv, &
10
       IntCoordI,Coordinate,CurrentCoord,ScanCoord,BprimT,BBT, &
11
       BBT_inv,Xprimitive,XPrimitiveF,Symmetry_elimination,UMat_local, &
12
       BTransInv_local, UMatF,BTransInvF,      &
13
       indzmat, dzdc,renum,order,OrderInv
14
  ! BMat_BakerT(3*Nat,NCoord), NCoord=3*Nat or NFree=3*Nat-6-Symmetry_elimination
15
  ! depending upon the coordinate choice. IntCoordI(NGeomI,NCoord) where 
16
  ! UMat(NPrim,NCoord), NCoord number of vectors in UMat matrix, i.e. NCoord
17
  ! Baker coordinates. NPrim is the number of primitive internal coordinates.
18

    
19
  Use Io_module
20
  IMPLICIT NONE
21

    
22
  ! IGeom: Geometry index, IGeomRef: Reference Geometry.
23
  INTEGER(KINT), INTENT(IN) :: atome(Nat),IGeomRef
24
  REAL(KREAL), INTENT(IN) :: fact
25
  REAL(KREAL), INTENT(IN)  :: r_cov(0:Max_Z)
26
  ! Frozen contains the indices of frozen atoms
27
  INTEGER(KINT), INTENT(IN) :: Frozen(*)
28

    
29
  INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:) ! (Nat,0:NMaxL)
30
  REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat)
31
  ! NPrim is the number of primitive coordinates and NCoord is the number
32
  ! of internal coordinates. BMat is actually (NPrim,3*Nat).
33
  REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim)
34
  ! EigVec(..) contains ALL eigevectors of BMat times BprimT, NOT only Baker Coordinate vectors.
35
  REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim)
36
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
37
  REAL(KREAL), ALLOCATABLE :: x_refGeom(:), y_refGeom(:), z_refGeom(:)
38
  !REAL(KREAL), ALLOCATABLE :: V(:,:) ! (3*Nat,3*Nat)
39
  !REAL(KREAL), ALLOCATABLE :: P(:) ! Used in Baker Coordinates: Matrix Inversion
40
  REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:)
41
  INTEGER(KINT) :: NbBonds, NbAngles, NbDihedrals, IGeom
42

    
43
  real(KREAL) :: vx, vy, vz, Norm
44
  real(KREAL) :: vx1,vy1,vz1,norm1
45
  real(KREAL) :: vx2,vy2,vz2,norm2
46
  real(KREAL) :: vx3,vy3,vz3,norm3
47
  real(KREAL) :: vx4,vy4,vz4,norm4
48
  real(KREAL) :: vx5,vy5,vz5,norm5
49

    
50
  INTEGER(KINT) :: I, J, IAt, K
51
  INTEGER(KINT) :: Kat, Lat, L
52

    
53
  REAL(KREAL) :: sAngleIatIKat, sAngleIIatLat
54

    
55
  LOGICAL :: debug, bond, AddPrimitiveCoord
56

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

    
63

    
64
     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
65
       use Path_module, only :  Pi,KINT, KREAL
66
       real(KREAL) ::  v1x,v1y,v1z,norm1
67
       real(KREAL) ::  v2x,v2y,v2z,norm2
68
       real(KREAL) ::  angle
69
     END FUNCTION ANGLE
70

    
71

    
72

    
73
     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
74
       use Path_module, only :  Pi,KINT, KREAL
75
       real(KREAL) ::  v1x,v1y,v1z,norm1
76
       real(KREAL) ::  v2x,v2y,v2z,norm2
77
       real(KREAL) ::  v3x,v3y,v3z,norm3
78
       real(KREAL) ::  angle_d,ca,sa
79
     END FUNCTION ANGLE_D
80

    
81

    
82

    
83

    
84
     SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
85
       !     
86
       ! This subroutine uses the description of a list of Coordinates
87
       ! to compute the values of the coordinates for a given geometry.
88
       !
89
!!!!!!!!!!
90
       ! Input:
91
       !   Na: INTEGER, Number of atoms
92
       !  x,y,z(Na): REAL, cartesian coordinates of the considered geometry
93
       ! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes
94
       ! NPrim, INTEGER: Number of coordinates to compute
95
       !
96
       ! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for
97
       ! a former geometry. Useful for Dihedral angles evolution...
98

    
99
!!!!!!!!!!!
100
       ! Output:
101
       ! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates
102
       !
103
!!!!!!!!!
104

    
105
       Use VarTypes
106
       Use Io_module
107
       Use Path_module, only : pi
108

    
109
       IMPLICIT NONE
110

    
111
       Type (ListCoord), POINTER :: Coordinate
112
       INTEGER(KINT), INTENT(IN) :: Nat,NPrim
113
       REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat)
114
       REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
115
       REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim)
116

    
117
     END SUBROUTINE CALC_XPRIM
118
  END INTERFACE
119

    
120
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121
  !
122
!!!!!!!! PFL Debug
123
  !
124
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
125
  INTEGER(KINT) :: I1,I2,I3,I4
126
  LOGICAL :: debugPFL
127
  REAL(KREAL), ALLOCATABLE :: val_zmat(:,:)
128
  INTEGER(KINT), ALLOCATABLE :: indzmat0(:,:)
129

    
130

    
131
  debug=valid("Calc_baker")
132
  debugPFL=valid("BakerPFL")
133

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

    
136
  IF (Nat.le.2) THEN
137
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
138
     RETURN
139
  END IF
140

    
141
  IF (debug) THEN
142
     WRITE(*,*) "DBG Calc_baker, geom"
143
     DO I=1,Nat
144
        WRITE(*,'(1X,I5,":",I3,3(1X,F15.3))') I,atome(I),XyzGeomI(IGeomRef,1:3,I)
145
     END DO
146
  END IF
147

    
148

    
149

    
150

    
151
  ! we construct the list of bonds, valence angles and dihedral angles using this connectivity
152
  ALLOCATE(Coordinate)
153
  NULLIFY(Coordinate%next)
154
  NbBonds=0
155
  NbAngles=0
156
  NbDihedrals=0
157

    
158
  CurrentCoord => Coordinate
159

    
160
  ALLOCATE(Liaisons(Nat,0:NMaxL),Geom(3,Nat),x(Nat),y(Nat),z(Nat))
161
  ALLOCATE(x_refGeom(Nat),y_refGeom(Nat),z_refGeom(Nat))
162

    
163
  x_refGeom(1:Nat) = XyzGeomI(IGeomRef,1,1:Nat)
164
  y_refGeom(1:Nat) = XyzGeomI(IGeomRef,2,1:Nat)
165
  z_refGeom(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat) 
166

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

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

    
216

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

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

    
227

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

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

    
294
           CurrentCoord%value = val_zmat(I,3)*Pi/180.
295
           CurrentCoord%SignDihedral = int(sign(1.D0,val_zmat(I,3)))
296
           CurrentCoord%At1=I1
297
           CurrentCoord%At2=I2
298
           CurrentCoord%At3=I3
299
           CurrentCoord%At4=I4
300
           CurrentCoord%Type="DIHEDRAL"
301
           IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
302
              ALLOCATE(CurrentCoord%next)
303
              NULLIFY(CurrentCoord%next%next)
304
           END IF
305
           CurrentCoord => CurrentCoord%next
306
        END IF ! matches IF (I4.NE.0)
307
     END DO ! matches DO I=1,Nat
308
     deallocate(indzmat0)
309

    
310
  ELSE !! if (debugPFL)
311

    
312
     DO IGeom=1, NGeomI
313
        !IGeom=1 ! need when using only first geometry.
314
        x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
315
        y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
316
        z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) 
317

    
318
        ! First step: we calculate the connectivity
319
        Liaisons=0
320

    
321
        Call CalcCnct(Nat,atome,x,y,z,LIAISONS,r_cov,fact)
322

    
323
        IF (debug) THEN 
324
           WRITE(*,*) 'Connectivity done'
325
           DO I=1,Nat
326
              WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,J),J=0,NMaxL)
327
              WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,:))
328
           END DO
329
           DO I=1,Nat
330
              WRITE(*,'(1X,I5,":",I3," r_cov=",F15.6)') I,atome(I),r_cov(atome(I))
331
           END DO
332
        END IF
333

    
334
        DO I=1,Nat
335
           IF (Liaisons(I,0).NE.0) THEN
336
              DO J=1,Liaisons(I,0)
337
                 ! We add the bond
338
                 Iat=Liaisons(I,J)
339
                 IF (Iat.GT.I) THEN
340
                    ! Checking the uniqueness of bond.
341
                    ! Duplicate bonds should not be added.
342
                    AddPrimitiveCoord=.True.
343
                    ScanCoord=>Coordinate
344
                    DO WHILE (Associated(ScanCoord%next))
345
                       IF (ScanCoord%Type .EQ. 'BOND') THEN
346
                          IF (ScanCoord%At1 .EQ. Iat) THEN
347
                             IF (ScanCoord%At2 .EQ. I) THEN
348
                                AddPrimitiveCoord=.False.
349
                             END IF
350
                          END IF
351
                       END IF
352
                       ScanCoord => ScanCoord%next
353
                    END DO
354

    
355
                    IF (AddPrimitiveCoord) THEN
356
                       NbBonds=NbBonds+1
357
                       ! values of the coordinate (bond) is calculated for the reference geometry
358
                       ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.)
359
                       ! are determined using all geometries. vecteur is defined in bib_oper.f
360
                       Call vecteur(I,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
361
                       CurrentCoord%value=Norm2
362
                       CurrentCoord%SignDihedral = 0
363
                       CurrentCoord%At1=Iat
364
                       CurrentCoord%At2=I
365
                       CurrentCoord%Type="BOND"
366
                       IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
367
                          ALLOCATE(CurrentCoord%next)
368
                          NULLIFY(CurrentCoord%next%next)
369
                       END IF
370
                       CurrentCoord => CurrentCoord%next
371
                    END IF ! matches IF (AddPrimitiveCoord) THEN
372
                 END IF ! matches IF (Iat.GT.I) THEN
373

    
374
                 !           Call Annul(Liaisons,IAt,I,Nat,NMaxL)
375
                 !           Liaisons(Iat,0)=Liaisons(Iat,0)-1
376

    
377
                 ! We add the valence angles
378
                 DO K=J+1,Liaisons(I,0)
379
                    Kat=Liaisons(I,K)
380
                    IF (Kat.GT.I) THEN
381
                       ! Checking the uniqueness of valence angle.
382
                       ! Duplicate bonds should not be added.
383
                       AddPrimitiveCoord=.True.
384
                       ScanCoord=>Coordinate
385
                       DO WHILE (Associated(ScanCoord%next))
386
                          IF (ScanCoord%Type .EQ. 'ANGLE') THEN
387
                             IF (ScanCoord%At1 .EQ. Iat) THEN
388
                                IF (ScanCoord%At2 .EQ. I) THEN
389
                                   IF (ScanCoord%At3 .EQ. Kat) THEN
390
                                      AddPrimitiveCoord=.False.
391
                                   END IF
392
                                END IF
393
                             END IF
394
                          END IF
395
                          ScanCoord => ScanCoord%next
396
                       END DO ! matches DO WHILE (Associated(ScanCoord%next))
397

    
398
                       IF (AddPrimitiveCoord) THEN      
399
                          IF (debug) WRITE(*,*) "Adding angle:",Iat,I,Kat
400
                          NbAngles=NbAngles+1
401
                          ! values of the coordinate (bond) is calculated for the reference geometry
402
                          ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
403
                          ! determined using all geometries.
404
                          Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
405
                          Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
406
                          CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
407
                          CurrentCoord%SignDihedral = 0
408
                          sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.))
409
                          CurrentCoord%At1=Iat
410
                          CurrentCoord%At2=I
411
                          CurrentCoord%At3=KAt
412
                          CurrentCoord%Type="ANGLE"
413
                          IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
414
                             ALLOCATE(CurrentCoord%next)
415
                             NULLIFY(CurrentCoord%next%next)
416
                          END IF
417
                          CurrentCoord => CurrentCoord%next
418
                       END IF ! matches IF (AddPrimitiveCoord) THEN
419
                    ELSE ! matches IF (Kat.GT.I) THEN
420
                       ! Kat is less than I. If there is a bond between I and Kat, then this angle
421
                       ! has already been taken into account.
422
                       if (debug) WRITE(*,*) "Testing angle:",Iat,I,Kat
423
                       Bond=.FALSE.
424
                       DO L=1,Liaisons(Iat,0)
425
                          IF (Liaisons(Iat,L).EQ.Kat) Bond=.TRUE.
426
                       END DO
427
                       IF (.NOT.Bond) THEN
428
                          IF (debug) WRITE(*,*) "Not Bond Adding angle:",Iat,I,Kat
429

    
430
                          ! Checking the uniqueness of valence angle.
431
                          ! Duplicate bonds should not be added.        
432
                          AddPrimitiveCoord=.True.
433
                          ScanCoord=>Coordinate
434
                          DO WHILE (Associated(ScanCoord%next))
435
                             IF (ScanCoord%Type .EQ. 'ANGLE') THEN
436
                                IF (ScanCoord%At1 .EQ. Iat) THEN
437
                                   IF (ScanCoord%At2 .EQ. I) THEN
438
                                      IF (ScanCoord%At3 .EQ. Kat) THEN
439
                                         AddPrimitiveCoord=.False.
440
                                      END IF
441
                                   END IF
442
                                END IF
443
                             END IF
444
                             ScanCoord => ScanCoord%next
445
                          END DO
446

    
447
                          IF (AddPrimitiveCoord) THEN  
448
                             NbAngles=NbAngles+1
449
                             ! values of the coordinate (bond) is calculated for the reference geometry
450
                             ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) 
451
                             ! are determined using all geometries.
452
                             Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
453
                             Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
454
                             CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
455
                             CurrentCoord%SignDihedral = 0
456
                             CurrentCoord%At1=Iat
457
                             CurrentCoord%At2=I
458
                             CurrentCoord%At3=KAt
459
                             CurrentCoord%Type="ANGLE"
460
                             IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
461
                                ALLOCATE(CurrentCoord%next)
462
                                NULLIFY(CurrentCoord%next%next)
463
                             END IF
464
                             CurrentCoord => CurrentCoord%next
465
                          END IF  ! matches IF (AddPrimitiveCoord) THEN
466
                       END IF ! matches IF (.NOT.Bond) THEN
467
                    END IF ! matches IF (Kat.GT.I) THEN, after else
468
                 END DO ! matches DO K=J+1,Liaisons(I,0)
469
              END DO ! matches DO J=1,Liaisons(I,0)
470
           END IF ! matches IF (Liaisons(I,0).NE.0) THEN
471

    
472
           ! We add dihedral angles. We do not concentrate on necessary angles, ie those that are
473
           ! needed to build the molecule. Instead we collect all of them, and the choice of the best
474
           ! ones will be done when diagonalizing the Baker matrix.
475
           IF (Liaisons(I,0).GE.3) Then
476
              DO J=1,Liaisons(I,0)
477
                 Iat=Liaisons(I,J)
478
                 ! values of the coordinate (bond) is calculated for the reference geometry
479
                 ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
480
                 ! determined using all geometries.
481
                 Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
482

    
483
                 DO K=J+1,Liaisons(I,0)
484
                    Kat=Liaisons(I,K)
485
                    Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
486
                    sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.))
487

    
488
                    DO L=K+1,Liaisons(I,0)
489
                       Lat=Liaisons(I,L)
490
                       if (debug) WRITE(*,*) "0-Dihe Kat,I,Iat:",Kat,I,Iat,angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)
491
                       Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
492
                       sAngleIIatLat=abs(sin(angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2)*Pi/180.))
493
                       if (debug) WRITE(*,*) "0-Dihe I,Iat,Lat:",angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2)
494
                       IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN
495
                          if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat
496

    
497
                          ! Checking for the uniqueness of valence angle.
498
                          ! Duplicate bonds should not be added.
499
                          AddPrimitiveCoord=.True.
500
                          ScanCoord=>Coordinate
501
                          DO WHILE (Associated(ScanCoord%next))
502
                             IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN
503
                                IF (ScanCoord%At1 .EQ. Kat) THEN
504
                                   IF (ScanCoord%At2 .EQ. I) THEN
505
                                      IF (ScanCoord%At3 .EQ. Iat) THEN
506
                                         IF (ScanCoord%At4 .EQ. Lat) THEN
507
                                            AddPrimitiveCoord=.False.
508
                                         END IF
509
                                      END IF
510
                                   END IF
511
                                END IF
512
                             END IF
513
                             ScanCoord => ScanCoord%next
514
                          END DO ! matches DO WHILE (Associated(ScanCoord%next))
515

    
516
                          ! IF (AddPrimitiveCoord) THEN
517
                          !   NbDihedrals=NbDihedrals+1
518
                          !  CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
519
                          !      vx5,vy5,vz5,norm5)
520
                          !CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
521
                          !    vx4,vy4,vz4,norm4)
522
                          !         CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
523
                          !             vx2,vy2,vz2,norm2)*Pi/180.
524
                          !          CurrentCoord%SignDihedral = 0
525
                          !       CurrentCoord%At1=Kat
526
                          !      CurrentCoord%At2=I
527
                          !     CurrentCoord%At3=IAt
528
                          !    CurrentCoord%At4=LAt
529
                          !WRITE(*,'(1X,":(dihed 1)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') &
530
                          !        ScanCoord%At1,ScanCoord%At2,ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
531
                          !         CurrentCoord%Type="DIHEDRAL"
532
                          !        IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
533
                          !          ALLOCATE(CurrentCoord%next)
534
                          !         NULLIFY(CurrentCoord%next%next)
535
                          !     END IF
536
                          !    CurrentCoord => CurrentCoord%next
537
                          !END IF ! matches IF (AddPrimitiveCoord) THEN
538

    
539
                       END IF ! matches IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN
540
                    END DO ! matches DO L=K+1,Liaisons(I,0)
541
                 END DO ! matches DO K=J+1,Liaisons(I,0)
542
              END DO ! matches DO J=1,Liaisons(I,0)
543
           END IF !matches  IF (Liaisons(I,0).GE.3) Then
544

    
545
           ! we now add other dihedrals
546
           DO J=1,Liaisons(I,0)
547
              Iat=Liaisons(I,J)
548
              If (Iat.LE.I) Cycle
549
              ! values of the coordinate (bond) is calculated for the reference geometry
550
              ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
551
              ! determined using all geometries.
552
              Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
553

    
554
              DO K=1,Liaisons(I,0)
555
                 IF (Liaisons(I,K).EQ.Iat) Cycle 
556
                 Kat=Liaisons(I,K)
557
                 Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
558
                 sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.))
559

    
560
                 DO L=1,Liaisons(Iat,0)
561
                    Lat=Liaisons(Iat,L)
562
                    IF (Lat.eq.I) Cycle
563

    
564
                    if (debug) WRITE(*,*) "Dihe Kat,I,Iat:",angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),sAngleIatIKat
565
                    Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx,vy,vz,Norm)
566

    
567
                    sAngleIIatLat=abs(sin(angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2)*Pi/180.))
568
                    if (debug) WRITE(*,*) "Dihe I,Iat,Lat:",angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2)
569
                    ! we add the dihedral angle Kat-I-Iat-Lat
570
                    if (debug) WRITE(*,*) "Trying dihedral:",Kat,I,Iat,Lat,sAngleIatIKat,sAngleIIatLat
571
                    IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) &
572
                         .AND.(sAngleIIatLat.ge.0.01d0)) THEN
573
                       if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat
574

    
575
                       ! Checking for the uniqueness of valence angle.
576
                       ! Duplicate bonds should not be added.
577
                       AddPrimitiveCoord=.True.
578
                       ScanCoord=>Coordinate
579
                       DO WHILE (Associated(ScanCoord%next))
580
                          IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN
581
                             IF (ScanCoord%At1 .EQ. Kat) THEN
582
                                IF (ScanCoord%At2 .EQ. I) THEN
583
                                   IF (ScanCoord%At3 .EQ. Iat) THEN
584
                                      IF (ScanCoord%At4 .EQ. Lat) THEN
585
                                         AddPrimitiveCoord=.False.
586
                                      END IF
587
                                   END IF
588
                                END IF
589
                             END IF
590
                          END IF
591
                          ScanCoord => ScanCoord%next
592
                       END DO
593

    
594
                       IF (AddPrimitiveCoord) THEN
595
                          NbDihedrals=NbDihedrals+1
596
                          Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
597
                          CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
598
                               vx5,vy5,vz5,norm5)
599
                          CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
600
                               vx4,vy4,vz4,norm4)
601
                          CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
602
                               vx2,vy2,vz2,norm2)*Pi/180.
603
                          CurrentCoord%At1=Kat
604
                          CurrentCoord%At2=I
605
                          CurrentCoord%At3=IAt
606
                          CurrentCoord%At4=LAt
607
                          CurrentCoord%Type="DIHEDRAL"
608
                          CurrentCoord%SignDihedral=int(sign(1.d0,angle_d(vx4,vy4,vz4,norm4,   &
609
                               vx5,vy5,vz5,norm5, vx2,vy2,vz2,norm2)))
610
                          WRITE(*,'(1X,":(dihed 2)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') ScanCoord%At1,&
611
                               ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
612
                          IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
613
                             ALLOCATE(CurrentCoord%next)
614
                             NULLIFY(CurrentCoord%next%next)
615
                          END IF
616
                          CurrentCoord => CurrentCoord%next
617
                       END IF ! matches IF (AddPrimitiveCoord) THEN
618
                    END IF ! matches IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) .....
619
                 END DO ! matches DO L=1,Liaisons(Iat,0)
620
              END DO ! matches DO K=1,Liaisons(I,0)
621
           END DO ! matches DO J=1,Liaisons(I,0)
622
        END DO ! end of loop over all atoms, DO I=1, Nat
623
     END DO ! matches DO IGeom=1, NGeomI
624

    
625

    
626

    
627
  END IF ! if (debugPFL)
628
  DEALLOCATE(Liaisons)
629

    
630
  WRITE(*,*) "I have found"
631
  WRITE(*,*) NbBonds, " bonds"
632
  WRITE(*,*) NbAngles, " angles"
633
  WRITE(*,*) NbDihedrals, " dihedrals"
634

    
635
  WRITE(*,*) 'All in all...'
636
  ScanCoord=>Coordinate
637
  I=0
638
  DO WHILE (Associated(ScanCoord%next))
639
     I=I+1
640
     SELECT CASE (ScanCoord%Type)
641
     CASE ('BOND')
642
        !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value
643
        WRITE(*,'(1X,I3,":(bond )",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value
644
     CASE ('ANGLE')
645
        !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &
646
        !    ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi
647
        WRITE(*,'(1X,I3,":(angle)",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &
648
             ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi   
649
     CASE ('DIHEDRAL')
650
        !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&
651
        !    ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
652
        WRITE(*,'(1X,I3,":(dihed)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&
653
             ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi
654
     END SELECT
655
     ScanCoord => ScanCoord%next
656
  END DO
657

    
658
  ! NPrim is the number of primitive coordinates.
659
  NPrim=NbBonds+NbAngles+NbDihedrals
660

    
661
  ! Now calculating values of all primitive bonds for all geometries:
662
  ALLOCATE(Xprimitive(NgeomI,NPrim),XPrimitiveF(NGeomF,NPrim))
663
  ALLOCATE(UMatF(NGeomF,NPrim,NCoord))
664
  ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
665

    
666
  ! We calculate the Primitive values for the first geometry, this plays a special role
667
  ! as it will serve as a reference for other geometries !
668

    
669
  x(1:Nat) = XyzGeomI(1,1,1:Nat)
670
  y(1:Nat) = XyzGeomI(1,2,1:Nat)
671
  z(1:Nat) = XyzGeomI(1,3,1:Nat)
672

    
673
  Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(1,:))
674

    
675
  DO IGeom=2, NGeomI
676

    
677
     x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
678
     y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
679
     z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
680

    
681
     Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(IGeom,:),Xprimitive(1,:))
682

    
683
     !      ScanCoord=>Coordinate
684
     !      I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
685
     !      DO WHILE (Associated(ScanCoord%next))
686
     !         I=I+1
687
     !         SELECT CASE (ScanCoord%Type)
688
     !         CASE ('BOND')
689
     !            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
690
     !            Xprimitive(IGeom,I) = Norm2  
691
     !         CASE ('ANGLE')
692
     !            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1)
693
     !            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
694
     !            Xprimitive(IGeom,I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
695
     !         CASE ('DIHEDRAL')
696

    
697

    
698
     !            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
699
     !            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)
700
     !            Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)
701
     !            Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
702
     !                 vx4,vy4,vz4,norm4)
703
     !            Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
704
     !                 vx5,vy5,vz5,norm5)
705

    
706
     !                  DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
707
     !                          vx2,vy2,vz2,norm2)
708
     !                  Xprimitive(IGeom,I) = DiheTmp*Pi/180.
709
     ! ! We treat large dihedral angles differently as +180=-180 mathematically and physically
710
     ! ! but this causes lots of troubles when converting baker to cart.
711
     ! ! So we ensure that large dihedral angles always have the same sign
712
     !                     if (abs(ScanCoord%SignDihedral).NE.1) THEN
713
     !                         ScanCoord%SignDihedral=Int(Sign(1.D0,DiheTmp))
714
     !                     ELSE
715
     !                         If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN
716
     !                           Xprimitive(IGeom,I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi
717
     !                        END IF
718
     !                   END IF
719
     ! ! Another test... will all this be consistent ??? 
720
     ! ! We want the shortest path, so we check that the change in dihedral angles is less than Pi:
721
     !                   IF (IGeom.GT.1) THEN
722
     !                      IF (abs(Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN
723
     !                         if ((Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN
724
     !                            Xprimitive(IGeom,I)=Xprimitive(IGeom,I)-2.d0*Pi
725
     !                         ELSE
726
     !                            Xprimitive(IGeom,I)=Xprimitive(IGeom,I)+2.d0*Pi
727
     !                         END IF
728
     !                      END IF
729
     !                   END IF
730
     !         END SELECT
731
     !         ScanCoord => ScanCoord%next
732
     !      END DO ! matches DO WHILE (Associated(ScanCoord%next))
733
  END DO ! matches DO IGeom=1, NGeomI
734

    
735
  IF (debug) THEN
736
     WRITE(*,*) "DBG Calc_Baker Xprimitives for all geometries"
737
     DO IGeom=1, NGeomI
738
        WRITE(*,*) "Geometry ",IGeom
739
        ScanCoord=>Coordinate
740
        I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
741
        DO WHILE (Associated(ScanCoord%next))
742
           I=I+1
743
           SELECT CASE (ScanCoord%Type)
744
           CASE ('BOND')
745
              WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive(IGeom,I)
746
           CASE ('ANGLE')
747
              WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &
748
                   ScanCoord%At2, ScanCoord%At3,Xprimitive(IGeom,I)*180./Pi
749
           CASE ('DIHEDRAL')
750
              WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&
751
                   ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive(IGeom,I)*180./Pi
752
           END SELECT
753
           ScanCoord => ScanCoord%next
754
        END DO ! matches DO WHILE (Associated(ScanCoord%next))
755
     END DO ! matches DO IGeom=1, NGeomI
756
  END IF
757
  ! Here reference geometry is assigned to Geom. Earlier x,y,z(Nat) were assigned from XyzGeomI(...)
758
  ! before the call of this (Calc_baker) subroutine and passed to this subroutine as arguments of the
759
  ! subroutine.
760
  ! PFL 19 Feb 2008 ->
761
  ! In order to reproduce the B matrix in the Baker article, one has
762
  ! to divide Geom by a0. 
763
  ! However, this is inconsitent with our way of storing geometries
764
  ! so we do not do it !
765

    
766
  Geom(1,:)=XyzGeomI(IGeomRef,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat)
767
  Geom(2,:)=XyzGeomI(IGeomRef,2,1:Nat)
768
  Geom(3,:)=XyzGeomI(IGeomRef,3,1:Nat)
769

    
770
  ! <- PFL 19 Feb 2008
771

    
772
  ! We now have to compute dq/dx...
773
  ALLOCATE(BprimT(3*Nat,NPrim))
774
  BprimT=0.d0
775

    
776
  ScanCoord=>Coordinate
777
  I=0
778
  DO WHILE (Associated(ScanCoord%next))
779
     I=I+1
780
     SELECT CASE (ScanCoord%Type)
781
     CASE ('BOND') 
782
        CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2, &
783
             Geom,BprimT(1,I))
784
     CASE ('ANGLE')
785
        CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2, &
786
             ScanCoord%At3,Geom,BprimT(1,I))
787
     CASE ('DIHEDRAL')
788
        CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2, &
789
             ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))
790
        !        CALL CONSTRAINTS_TORSION_DER(Nat,ScanCoord%At1,ScanCoord%AT2, &
791
        !             ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))
792

    
793
     END SELECT
794
     ScanCoord => ScanCoord%next
795
  END DO
796

    
797

    
798
!!!!!!!!!!!!!!!!!!!!
799
  !
800
  !   PFL Debug
801
  !
802
!!!!!!!!!!!!!!!!!!!!
803

    
804
  if (debugPFL) THEN
805
     WRITe(*,*) "DBG PFL Calc_Baker ====="
806
     DzDc=0.d0
807
     WRITE(*,*) "PFL DBG, DzdC"
808
     do I=1,Nat
809
        Geom(1,i)=XyzGeomI(IGeomRef,1,orderInv(i)) ! XyzGeomI(NGeomI,3,Nat)
810
        Geom(2,i)=XyzGeomI(IGeomRef,2,OrderInv(i))
811
        Geom(3,i)=XyzGeomI(IGeomRef,3,OrderInv(i))
812
        WRITE(*,*) I,OrderInv(I),Geom(:,I)
813
     end do
814
     CALL CalcBMat_int (nat, Geom, indzmat, dzdc, massat,atome)     
815

    
816
     WRITE(*,*) "PFL DBG, BPrimT"
817
     DO I=1,NPrim
818
        WRITE(*,'(50(1X,F12.6))') BPrimT(:,I)
819
     END DO
820
     WRITe(*,*) "DBG PFL Calc_Baker ====="
821
  END IF
822

    
823
!!!!!!!!!!!!!!!!!!!!
824
  !
825
  !   PFL Debug
826
  !
827
!!!!!!!!!!!!!!!!!!!!
828

    
829
  DEALLOCATE(Geom)
830

    
831
  ! BprimT(3*Nat,NPrim)
832
  ! We now compute G=B(BT) matrix
833
  ALLOCATE(Gmat(NPrim,NPrim))
834
  GMat=0.d0
835
  DO I=1,NPrim
836
     DO J=1,3*Nat
837
        GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0)))
838
     END DO
839
  END DO
840

    
841
!!!!!!!!!!!!!!!!!!!!
842
  !
843
  !   PFL Debug --->
844
  !
845
!!!!!!!!!!!!!!!!!!!!
846
  if (debugPFL) THEN
847
     GMat=0.d0
848
     DO I=1,NPrim
849
        GMat(I,I)=1.
850
     END DO
851
  END IF
852

    
853
!!!!!!!!!!!!!!!!!!!!
854
  !
855
  !  <--- PFL Debug
856
  !
857
!!!!!!!!!!!!!!!!!!!!
858

    
859
  !DO I=1,NPrim
860
  !  WRITE(*,'(20(1X,F6.3))') GMat(1:min(20,NPrim),I)
861
  !END DO
862

    
863
  ! We diagonalize G
864
  ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim))
865
  Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim)
866
  Call Trie(NPrim,EigVal,EigVec,NPrim)
867
  DO I=1,NPrim
868
     WRITE(*,'(1X,"Vector ",I3,": e=",F8.3)') I,EigVal(i)
869
     WRITE(*,'(20(1X,F8.4))') EigVec(1:min(20,NPrim),I)
870
  END DO
871

    
872
  !DO I=1,NPrim
873
  ! WRITE(*,'(1X,"Vector ",I3,(20(F8.3)))') I,EigVec(1:min(20,NPrim),I)
874
  ! END DO
875

    
876
  !DO I=1,NPrim
877
  !  WRITE(*,'(1X,"Vector ",I3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,&
878
  ! F8.3,F8.3,F8.3,F8.3,F8.3)') I,EigVec(1,I),EigVec(4,I),EigVec(6,I),EigVec(14,I),EigVec(16,I),&
879
  !EigVec(2,I),EigVec(3,I),EigVec(5,I),EigVec(12,I),EigVec(13,I),EigVec(15,I),EigVec(7,I),&
880
  !EigVec(8,I),EigVec(9,I),EigVec(10,I),EigVec(11,I),EigVec(17,I)
881
  !END DO
882

    
883
  ! We construct the new DzDc(b=baker) matrix
884
  ! UMat is nonredundant vector set, i.e. set of eigenvectors of BB^T corresponding
885
  ! to eigenvalues > zero.
886
  ALLOCATE(UMat(NPrim,NCoord))
887
  ALLOCATE(UMat_local(NPrim,NCoord))
888
  ALLOCATE(UMat_tmp(NPrim,3*Nat-6))
889
  ! BMat_BakerT(3*Nat,NCoord), allocated in Path.f90, 
890
  ! NCoord=3*Nat-6
891
  BMat_BakerT = 0.d0
892
  J=0
893
  UMat=0.d0
894
  UMat_tmp=0.d0
895
  DO I=1,NPrim
896
     IF (abs(eigval(I))>=1e-6) THEN
897
        J=J+1
898
        DO K=1,NPrim
899
           !DzDcb(:,J)=DzDcb(:,J)+Eigvec(K,I)*BprimT(:,K) ! original
900
           ! BprimT is transpose of B^prim.
901
           ! B = UMat^T B^prim, B^T = (B^prim)^T UMat
902
           BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I)
903
           !Dzdcb(:,J)=DzDcb(:,J)+Eigvec(K,:)*BprimT(I,K)       
904
        END DO
905
        IF(J .GT. 3*Nat-6) THEN
906
           WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. 1e-6(=UMat) (=' &
907
                ,J,') exceeded 3*Nat-6=',3*Nat-6, &
908
                'Stopping the calculation.'
909
           STOP
910
        END IF
911
        UMat_tmp(:,J) =  Eigvec(:,I)
912
     END IF
913
  END DO
914

    
915
  if (debug) THEN
916
     WRITE(*,*) "DBG Calc_Baker: UMat"
917
     DO I=1,3*Nat-6
918
        WRITE(*,'(50(1X,F12.8))') UMat_tmp(:,I)
919
     END DO
920
  END IF
921

    
922
  ! THIS IS TO BE CORRECTED:
923
  UMat = UMat_tmp
924
  !J=0
925
  !DO I=1, 3*Nat-6
926
  !  IF ((I .NE. 6) .AND. (I .NE. 7) .AND. (I .NE. 9)) Then
927
  !   J=J+1
928
  !  Print *, 'J=', J, ' I=', I
929
  ! UMat(:,J) = UMat_tmp(:,I)
930
  !END IF
931
  !END DO
932
  ! BMat_BakerT=0.d0
933
  !DO I=1,NCoord
934
  !  DO J=1,NPrim
935
  ! B = UMat^T B^prim, B^T = (B^prim)^T UMat
936
  !   BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat(J,I)
937
  !    END DO
938
  !END DO
939
  ! TO BE CORRECTED, ENDS HERE.
940

    
941
  IntCoordI=0.d0
942
  DO IGeom=1, NGeomI
943
     DO I=1, NPrim
944
        ! Transpose of UMat is needed below, that is why UMat(I,:).
945
        IntCoordI(IGeom,:) = IntCoordI(IGeom,:) + UMat(I,:)*Xprimitive(IGeom,I)
946
     END DO
947
  END DO
948

    
949
  if (debug) THEN
950
     WRITE(*,*) "DBG Calc_Baker IntCoordI"
951
     DO IGeom=1, NGeomI
952
        WRITE(*,'(1X,I5,":",30(1X,F10.6))') IGeom,IntCoordI(IGeom,:)
953
     END DO
954
  END IF
955

    
956

    
957
  ! the following might be used to get rid of some coordinates that 
958
  ! do not respect the symmetry of the problem.
959

    
960
  ! ! We look at the group symmetry of the molecule. Nothing complicate here
961
  ! ! we just look for planar symmetry.
962
  ! ! We first place it into the standard orientation
963
  !   Call Std_ori(Nat,x,y,z,a,b,c,massat)
964

    
965
  ! ! Is the molecule planar
966
  !   DO I=1,Nat
967
  !      WRITE(*,*) nom(atome(i)),x(i),y(i),z(i)
968
  !      END DO
969
  DEALLOCATE(BprimT,GMat,EigVal,EigVec) ! BprimT is B^prim^T
970
  ! Calculation of BTransInv starts here:
971
  ! Calculation of BBT(3*Nat-6,3*Nat-6)=BB^T:
972
  ! BMat_BakerT(3*Nat,NCoord) is Transpose of B = UMat^TB^prim
973
  ALLOCATE(BBT(NCoord,NCoord))
974
  BBT = 0.d0
975
  DO I=1, NCoord
976
     DO J=1, 3*Nat
977
        ! BBT(:,I) forms BB^T
978
        BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)
979
     END DO
980
  END DO
981

    
982
  !ALLOCATE(V(3*Nat,3*Nat))
983
  !Call GINVSE(BBT,3*Nat,3*Nat,BBT_inv,3*Nat,3*Nat,3*Nat,1e-6,V,3*Nat,FAIL,NOUT)
984
  !DEALLOCATE(V)
985
  ALLOCATE(BBT_inv(NCoord,NCoord))
986
  ! GenInv is in Mat_util.f90
987
  Call GenInv(NCoord,BBT,BBT_inv,NCoord)
988

    
989
  ! Calculation of (B^T)^-1 = (BB^T)^-1B:
990
  !ALLOCATE(BTransInv(3*Nat-6,3*Nat))    ! allocated in Path.f90
991
  BTransInv = 0.d0
992
  DO I=1, 3*Nat
993
     DO J=1, NCoord
994
        BTransInv(:,I) = BTransInv(:,I) + BBT_inv(:,J)*BMat_BakerT(I,J)
995
     END DO
996
  END DO
997
  BTransInv_local = BTransInv
998
  UMat_local = UMat
999
  DEALLOCATE(BBT,BBT_inv,UMat_tmp)
1000
  DEALLOCATE(x_refGeom,y_refGeom,z_refGeom,x,y,z)
1001
  IF (debug) WRITE(*,*) "DBG Calc_baker over."
1002
END SUBROUTINE Calc_baker