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(*,*) "0Dihe 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(*,*) "0Dihe 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(*,*) "0Dihe 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(*,*) "0Dihe 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 KatIIatLat 

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 KatIIatLat 

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(IGeom1,I)).GE.Pi) THEN 

722 
! if ((Xprimitive(IGeom,I)XPrimitive(IGeom1,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(IGeom1,I)).GE.Pi) THEN


723 
! if ((Xprimitive(IGeom,I)XPrimitive(IGeom1,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