Révision 7 src/Calc_baker.f90

Calc_baker.f90 (revision 7)
78 78
       real(KREAL) ::  angle_d,ca,sa
79 79
     END FUNCTION ANGLE_D
80 80

  
81
 
82 81

  
83 82

  
83

  
84 84
     SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
85 85
       !     
86 86
       ! This subroutine uses the description of a list of Coordinates
......
118 118
  END INTERFACE
119 119

  
120 120
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121
!
121
  !
122 122
!!!!!!!! PFL Debug
123
!
123
  !
124 124
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
125 125
  INTEGER(KINT) :: I1,I2,I3,I4
126 126
  LOGICAL :: debugPFL
......
130 130

  
131 131
  debug=valid("Calc_baker")
132 132
  debugPFL=valid("BakerPFL")
133
  if (debug) WRITE(*,*) "============= Entering Cal_baker =============="
134 133

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

  
135 136
  IF (Nat.le.2) THEN
136 137
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
137 138
     RETURN
......
164 165
  z_refGeom(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat) 
165 166

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

  
199
         DO I=1,Nat
200
           Order(IndZmat0(I,1))=I
201
           OrderInv(I)=IndZmat0(I,1)
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))
202 213
        END DO
203
        IndZmat(1,1)=Order(IndZmat0(1,1))
204
        IndZmat(2,1)=Order(IndZmat0(2,1))
205
        IndZmat(2,2)=Order(IndZmat0(2,2))
206
        IndZmat(3,1)=Order(IndZmat0(3,1))
207
        IndZmat(3,2)=Order(IndZmat0(3,2))
208
        IndZmat(3,3)=Order(IndZmat0(3,3))
209
        DO I=4,Nat
210
           DO J=1,4
211
              IndZmat(I,J)=Order(IndZmat0(I,J))
212
           END DO
213
        end do
214
     end do
214 215

  
215 216

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

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

  
226 227

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

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

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

  
309 310
  ELSE !! if (debugPFL)
310 311

  
311
  DO IGeom=1, NGeomI
312
     !IGeom=1 ! need when using only first geometry.
313
     x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
314
     y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
315
     z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) 
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) 
316 317

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

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

  
322
     IF (debug) THEN 
323
        WRITE(*,*) 'Connectivity done'
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

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

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

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

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

  
376
              ! We add the valence angles
377
              DO K=J+1,Liaisons(I,0)
378
                 Kat=Liaisons(I,K)
379
                 IF (Kat.GT.I) THEN
380
                    ! Checking the uniqueness of valence angle.
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.
381 341
                    ! Duplicate bonds should not be added.
382 342
                    AddPrimitiveCoord=.True.
383 343
                    ScanCoord=>Coordinate
384 344
                    DO WHILE (Associated(ScanCoord%next))
385
                       IF (ScanCoord%Type .EQ. 'ANGLE') THEN
345
                       IF (ScanCoord%Type .EQ. 'BOND') THEN
386 346
                          IF (ScanCoord%At1 .EQ. Iat) THEN
387 347
                             IF (ScanCoord%At2 .EQ. I) THEN
388
                                IF (ScanCoord%At3 .EQ. Kat) THEN
389
                                   AddPrimitiveCoord=.False.
390
                                END IF
348
                                AddPrimitiveCoord=.False.
391 349
                             END IF
392 350
                          END IF
393 351
                       END IF
394 352
                       ScanCoord => ScanCoord%next
395
                    END DO ! matches DO WHILE (Associated(ScanCoord%next))
353
                    END DO
396 354

  
397
                    IF (AddPrimitiveCoord) THEN      
398
                       IF (debug) WRITE(*,*) "Adding angle:",Iat,I,Kat
399
                       NbAngles=NbAngles+1
355
                    IF (AddPrimitiveCoord) THEN
356
                       NbBonds=NbBonds+1
400 357
                       ! values of the coordinate (bond) is calculated for the reference geometry
401
                       ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
402
                       ! determined using all geometries.
403
                       Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
404
                       Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
405
                       CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
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
406 362
                       CurrentCoord%SignDihedral = 0
407
                       sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.))
408 363
                       CurrentCoord%At1=Iat
409 364
                       CurrentCoord%At2=I
410
                       CurrentCoord%At3=KAt
411
                       CurrentCoord%Type="ANGLE"
365
                       CurrentCoord%Type="BOND"
412 366
                       IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN
413 367
                          ALLOCATE(CurrentCoord%next)
414 368
                          NULLIFY(CurrentCoord%next%next)
415 369
                       END IF
416 370
                       CurrentCoord => CurrentCoord%next
417 371
                    END IF ! matches IF (AddPrimitiveCoord) THEN
418
                 ELSE ! matches IF (Kat.GT.I) THEN
419
                    ! Kat is less than I. If there is a bond between I and Kat, then this angle
420
                    ! has already been taken into account.
421
                    if (debug) WRITE(*,*) "Testing angle:",Iat,I,Kat
422
                    Bond=.FALSE.
423
                    DO L=1,Liaisons(Iat,0)
424
                       IF (Liaisons(Iat,L).EQ.Kat) Bond=.TRUE.
425
                    END DO
426
                    IF (.NOT.Bond) THEN
427
                       IF (debug) WRITE(*,*) "Not Bond Adding angle:",Iat,I,Kat
372
                 END IF ! matches IF (Iat.GT.I) THEN
428 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
429 381
                       ! Checking the uniqueness of valence angle.
430
                       ! Duplicate bonds should not be added.        
382
                       ! Duplicate bonds should not be added.
431 383
                       AddPrimitiveCoord=.True.
432 384
                       ScanCoord=>Coordinate
433 385
                       DO WHILE (Associated(ScanCoord%next))
......
441 393
                             END IF
442 394
                          END IF
443 395
                          ScanCoord => ScanCoord%next
444
                       END DO
396
                       END DO ! matches DO WHILE (Associated(ScanCoord%next))
445 397

  
446
                       IF (AddPrimitiveCoord) THEN  
398
                       IF (AddPrimitiveCoord) THEN      
399
                          IF (debug) WRITE(*,*) "Adding angle:",Iat,I,Kat
447 400
                          NbAngles=NbAngles+1
448 401
                          ! values of the coordinate (bond) is calculated for the reference geometry
449
                          ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) 
450
                          ! are determined using all geometries.
402
                          ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are 
403
                          ! determined using all geometries.
451 404
                          Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
452 405
                          Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
453 406
                          CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
454
                           CurrentCoord%SignDihedral = 0
407
                          CurrentCoord%SignDihedral = 0
408
                          sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.))
455 409
                          CurrentCoord%At1=Iat
456 410
                          CurrentCoord%At2=I
457 411
                          CurrentCoord%At3=KAt
......
461 415
                             NULLIFY(CurrentCoord%next%next)
462 416
                          END IF
463 417
                          CurrentCoord => CurrentCoord%next
464
                       END IF  ! matches IF (AddPrimitiveCoord) THEN
465
                    END IF ! matches IF (.NOT.Bond) THEN
466
                 END IF ! matches IF (Kat.GT.I) THEN, after else
467
              END DO ! matches DO K=J+1,Liaisons(I,0)
468
           END DO ! matches DO J=1,Liaisons(I,0)
469
        END IF ! matches IF (Liaisons(I,0).NE.0) THEN
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
470 429

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

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

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

  
496 575
                       ! Checking for the uniqueness of valence angle.
......
510 589
                             END IF
511 590
                          END IF
512 591
                          ScanCoord => ScanCoord%next
513
                       END DO ! matches DO WHILE (Associated(ScanCoord%next))
592
                       END DO
514 593

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

  
538
                    END IF ! matches IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN
539
                 END DO ! matches DO L=K+1,Liaisons(I,0)
540
              END DO ! matches DO K=J+1,Liaisons(I,0)
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)
541 621
           END DO ! matches DO J=1,Liaisons(I,0)
542
        END IF !matches  IF (Liaisons(I,0).GE.3) Then
622
        END DO ! end of loop over all atoms, DO I=1, Nat
623
     END DO ! matches DO IGeom=1, NGeomI
543 624

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

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

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

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

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

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

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

  
624

  
625

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

  
629 630
  WRITE(*,*) "I have found"
......
662 663
  ALLOCATE(UMatF(NGeomF,NPrim,NCoord))
663 664
  ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
664 665

  
665
! We calculate the Primitive values for the first geometry, this plays a special role
666
! as it will serve as a reference for other geometries !
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 !
667 668

  
668 669
  x(1:Nat) = XyzGeomI(1,1,1:Nat)
669 670
  y(1:Nat) = XyzGeomI(1,2,1:Nat)
......
673 674

  
674 675
  DO IGeom=2, NGeomI
675 676

  
676
        x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
677
        y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
678
        z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
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)
679 680

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

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

  
696 697

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

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

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

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

  
769
! <- PFL 19 Feb 2008
770
  ! <- PFL 19 Feb 2008
770 771

  
771 772
  ! We now have to compute dq/dx...
772 773
  ALLOCATE(BprimT(3*Nat,NPrim))
......
786 787
     CASE ('DIHEDRAL')
787 788
        CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2, &
788 789
             ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))
789
!        CALL CONSTRAINTS_TORSION_DER(Nat,ScanCoord%At1,ScanCoord%AT2, &
790
!             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))
791 792

  
792 793
     END SELECT
793 794
     ScanCoord => ScanCoord%next
......
795 796

  
796 797

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

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

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

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

  
828 829
  DEALLOCATE(Geom)
......
838 839
  END DO
839 840

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

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

  
858 859
  !DO I=1,NPrim
......
876 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,&
877 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),&
878 879
  !EigVec(2,I),EigVec(3,I),EigVec(5,I),EigVec(12,I),EigVec(13,I),EigVec(15,I),EigVec(7,I),&
879
       !EigVec(8,I),EigVec(9,I),EigVec(10,I),EigVec(11,I),EigVec(17,I)
880
  !EigVec(8,I),EigVec(9,I),EigVec(10,I),EigVec(11,I),EigVec(17,I)
880 881
  !END DO
881 882

  
882 883
  ! We construct the new DzDc(b=baker) matrix

Formats disponibles : Unified diff