root / src / Calc_baker.f90 @ 8
Historique | Voir | Annoter | Télécharger (40,75 ko)
1 |
SUBROUTINE Calc_baker(atome,r_cov,fact,frozen,IGeomRef) |
---|---|
2 |
! |
3 |
! This subroutine analyses a geometry to construct the baker |
4 |
! delocalized internal coordinates |
5 |
! v1.0 |
6 |
! We use only one geometry |
7 |
! |
8 |
Use Path_module, only : max_Z,NMaxL,Pi,Massat,BMat_BakerT, & |
9 |
Nat,NCoord,XyzGeomI,NGeomI,NGeomF,UMat,NPrim,BTransInv, & |
10 |
IntCoordI,Coordinate,CurrentCoord,ScanCoord,BprimT,BBT, & |
11 |
BBT_inv,Xprimitive,XPrimitiveF,UMat_local, & |
12 |
BTransInv_local, UMatF,BTransInvF, & |
13 |
indzmat, dzdc,renum,order,OrderInv |
14 |
! BMat_BakerT(3*Nat,NCoord), NCoord=3*Nat or NFree=3*Nat-6-Symmetry_elimination |
15 |
! depending upon the coordinate choice. IntCoordI(NGeomI,NCoord) where |
16 |
! UMat(NPrim,NCoord), NCoord number of vectors in UMat matrix, i.e. NCoord |
17 |
! Baker coordinates. NPrim is the number of primitive internal coordinates. |
18 |
|
19 |
Use Io_module |
20 |
IMPLICIT NONE |
21 |
|
22 |
! IGeom: Geometry index, IGeomRef: Reference Geometry. |
23 |
INTEGER(KINT), INTENT(IN) :: atome(Nat),IGeomRef |
24 |
REAL(KREAL), INTENT(IN) :: fact |
25 |
REAL(KREAL), INTENT(IN) :: r_cov(0:Max_Z) |
26 |
! Frozen contains the indices of frozen atoms |
27 |
INTEGER(KINT), INTENT(IN) :: Frozen(*) |
28 |
|
29 |
INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:) ! (Nat,0:NMaxL) |
30 |
REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat) |
31 |
! NPrim is the number of primitive coordinates and NCoord is the number |
32 |
! of internal coordinates. BMat is actually (NPrim,3*Nat). |
33 |
REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim) |
34 |
! EigVec(..) contains ALL eigevectors of BMat times BprimT, NOT only Baker Coordinate vectors. |
35 |
REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim) |
36 |
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) |
37 |
REAL(KREAL), ALLOCATABLE :: x_refGeom(:), y_refGeom(:), z_refGeom(:) |
38 |
!REAL(KREAL), ALLOCATABLE :: V(:,:) ! (3*Nat,3*Nat) |
39 |
!REAL(KREAL), ALLOCATABLE :: P(:) ! Used in Baker Coordinates: Matrix Inversion |
40 |
REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:) |
41 |
INTEGER(KINT) :: NbBonds, NbAngles, NbDihedrals, IGeom |
42 |
|
43 |
real(KREAL) :: vx, vy, vz, Norm |
44 |
real(KREAL) :: vx1,vy1,vz1,norm1 |
45 |
real(KREAL) :: vx2,vy2,vz2,norm2 |
46 |
real(KREAL) :: vx3,vy3,vz3,norm3 |
47 |
real(KREAL) :: vx4,vy4,vz4,norm4 |
48 |
real(KREAL) :: vx5,vy5,vz5,norm5 |
49 |
|
50 |
INTEGER(KINT) :: I, J, IAt, K |
51 |
INTEGER(KINT) :: Kat, Lat, L |
52 |
|
53 |
REAL(KREAL) :: sAngleIatIKat, sAngleIIatLat |
54 |
|
55 |
LOGICAL :: debug, bond, AddPrimitiveCoord |
56 |
|
57 |
INTERFACE |
58 |
function valid(string) result (isValid) |
59 |
CHARACTER(*), intent(in) :: string |
60 |
logical :: isValid |
61 |
END function VALID |
62 |
|
63 |
|
64 |
FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2) |
65 |
use Path_module, only : Pi,KINT, KREAL |
66 |
real(KREAL) :: v1x,v1y,v1z,norm1 |
67 |
real(KREAL) :: v2x,v2y,v2z,norm2 |
68 |
real(KREAL) :: angle |
69 |
END FUNCTION ANGLE |
70 |
|
71 |
|
72 |
|
73 |
FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3) |
74 |
use Path_module, only : Pi,KINT, KREAL |
75 |
real(KREAL) :: v1x,v1y,v1z,norm1 |
76 |
real(KREAL) :: v2x,v2y,v2z,norm2 |
77 |
real(KREAL) :: v3x,v3y,v3z,norm3 |
78 |
real(KREAL) :: angle_d,ca,sa |
79 |
END FUNCTION ANGLE_D |
80 |
|
81 |
|
82 |
|
83 |
|
84 |
SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef) |
85 |
! |
86 |
! This subroutine uses the description of a list of Coordinates |
87 |
! to compute the values of the coordinates for a given geometry. |
88 |
! |
89 |
!!!!!!!!!! |
90 |
! Input: |
91 |
! Na: INTEGER, Number of atoms |
92 |
! x,y,z(Na): REAL, cartesian coordinates of the considered geometry |
93 |
! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes |
94 |
! NPrim, INTEGER: Number of coordinates to compute |
95 |
! |
96 |
! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for |
97 |
! a former geometry. Useful for Dihedral angles evolution... |
98 |
|
99 |
!!!!!!!!!!! |
100 |
! Output: |
101 |
! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates |
102 |
! |
103 |
!!!!!!!!! |
104 |
|
105 |
Use VarTypes |
106 |
Use Io_module |
107 |
Use Path_module, only : pi |
108 |
|
109 |
IMPLICIT NONE |
110 |
|
111 |
Type (ListCoord), POINTER :: Coordinate |
112 |
INTEGER(KINT), INTENT(IN) :: Nat,NPrim |
113 |
REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat) |
114 |
REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) |
115 |
REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim) |
116 |
|
117 |
END SUBROUTINE CALC_XPRIM |
118 |
END INTERFACE |
119 |
|
120 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
121 |
! |
122 |
!!!!!!!! PFL Debug |
123 |
! |
124 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
125 |
INTEGER(KINT) :: I1,I2,I3,I4 |
126 |
LOGICAL :: debugPFL |
127 |
REAL(KREAL), ALLOCATABLE :: val_zmat(:,:) |
128 |
INTEGER(KINT), ALLOCATABLE :: indzmat0(:,:) |
129 |
|
130 |
|
131 |
debug=valid("Calc_baker") |
132 |
debugPFL=valid("BakerPFL") |
133 |
|
134 |
if (debug) Call Header(" Entering Cal_baker ") |
135 |
|
136 |
IF (Nat.le.2) THEN |
137 |
WRITE(*,*) "I do not work for less than 2 atoms :-p" |
138 |
RETURN |
139 |
END IF |
140 |
|
141 |
IF (debug) THEN |
142 |
WRITE(*,*) "DBG Calc_baker, geom" |
143 |
DO I=1,Nat |
144 |
WRITE(*,'(1X,I5,":",I3,3(1X,F15.3))') I,atome(I),XyzGeomI(IGeomRef,1:3,I) |
145 |
END DO |
146 |
END IF |
147 |
|
148 |
|
149 |
|
150 |
|
151 |
! we construct the list of bonds, valence angles and dihedral angles using this connectivity |
152 |
ALLOCATE(Coordinate) |
153 |
NULLIFY(Coordinate%next) |
154 |
NbBonds=0 |
155 |
NbAngles=0 |
156 |
NbDihedrals=0 |
157 |
|
158 |
CurrentCoord => Coordinate |
159 |
|
160 |
ALLOCATE(Liaisons(Nat,0:NMaxL),Geom(3,Nat),x(Nat),y(Nat),z(Nat)) |
161 |
ALLOCATE(x_refGeom(Nat),y_refGeom(Nat),z_refGeom(Nat)) |
162 |
|
163 |
x_refGeom(1:Nat) = XyzGeomI(IGeomRef,1,1:Nat) |
164 |
y_refGeom(1:Nat) = XyzGeomI(IGeomRef,2,1:Nat) |
165 |
z_refGeom(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat) |
166 |
|
167 |
!!!!!!!!!!!!!!!!!!!! |
168 |
! |
169 |
! PFL Debug |
170 |
! |
171 |
!!!!!!!!!!!!!!!!!!!! |
172 |
if (debugPFL) THEN |
173 |
WRITE(*,*) "DBG PFL Calc_BAKER - Calculating Zmat" |
174 |
if (.not.ALLOCATED(indzmat)) THEN |
175 |
ALLOCATE(indzmat(Nat,5)) |
176 |
END IF |
177 |
IF (.NOT.ALLOCATED(DzDc)) THEN |
178 |
ALLOCATE(DzDc(3,nat,3,Nat)) |
179 |
END IF |
180 |
IF (.NOT.ALLOCATED(Val_zmat)) THEN |
181 |
ALLOCATE(val_zmat(nat,3)) |
182 |
END IF |
183 |
IF (.NOT.ALLOCATED(indzmat0)) THEN |
184 |
ALLOCATE(indzmat0(nat,5)) |
185 |
END IF |
186 |
! We construct a Zmat |
187 |
IF (Frozen(1).NE.0) THEN |
188 |
Renum=.TRUE. |
189 |
!!! remplcaer indzmat par indzmat0 puis renumerotation |
190 |
call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat, & |
191 |
x_refGeom,y_refGeom,z_refGeom, & |
192 |
r_cov,fact,frozen) |
193 |
ELSE |
194 |
!no zmat... we create our own |
195 |
call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat, & |
196 |
x_refGeom,y_refGeom,z_refGeom, & |
197 |
r_cov,fact) |
198 |
END IF |
199 |
|
200 |
DO I=1,Nat |
201 |
Order(IndZmat0(I,1))=I |
202 |
OrderInv(I)=IndZmat0(I,1) |
203 |
END DO |
204 |
IndZmat(1,1)=Order(IndZmat0(1,1)) |
205 |
IndZmat(2,1)=Order(IndZmat0(2,1)) |
206 |
IndZmat(2,2)=Order(IndZmat0(2,2)) |
207 |
IndZmat(3,1)=Order(IndZmat0(3,1)) |
208 |
IndZmat(3,2)=Order(IndZmat0(3,2)) |
209 |
IndZmat(3,3)=Order(IndZmat0(3,3)) |
210 |
DO I=4,Nat |
211 |
DO J=1,4 |
212 |
IndZmat(I,J)=Order(IndZmat0(I,J)) |
213 |
END DO |
214 |
end do |
215 |
|
216 |
|
217 |
WRITE(*,*) "DBG PFL CalcBaker, IndZmat0" |
218 |
DO I=1,Nat |
219 |
WRITE(*,*) I, indZmat0(I,:) |
220 |
end do |
221 |
|
222 |
WRITE(*,*) "DBG PFL CalcBaker, IndZmat" |
223 |
DO I=1,Nat |
224 |
WRITE(*,*) I, indZmat(I,:) |
225 |
end do |
226 |
|
227 |
|
228 |
DO I=1,Nat |
229 |
I1=indzmat0(I,1) |
230 |
I2=indzmat0(I,2) |
231 |
I3=indzmat0(I,3) |
232 |
I4=indzmat0(I,4) |
233 |
! Adding BOND between at1 and at2 |
234 |
WRITE(*,*) "DBG Indzmat",I,I1,I2,I3,I4 |
235 |
if (I2.NE.0) THEN |
236 |
NbBonds=NbBonds+1 |
237 |
! values of the coordinate (bond) is calculated for the reference geometry |
238 |
! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) |
239 |
! are determined using all geometries. vecteur is defined in bib_oper.f |
240 |
Call vecteur(I1,I2,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
241 |
! CurrentCoord%value=Norm2 |
242 |
WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Bond",NbBonds," Norm2,val_zmat",Norm2,val_zmat(I,1) |
243 |
CurrentCoord%value=val_zmat(I,1) |
244 |
CurrentCoord%SignDihedral = 0 |
245 |
CurrentCoord%At1=I1 |
246 |
CurrentCoord%At2=I2 |
247 |
CurrentCoord%Type="BOND" |
248 |
IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
249 |
ALLOCATE(CurrentCoord%next) |
250 |
NULLIFY(CurrentCoord%next%next) |
251 |
END IF |
252 |
CurrentCoord => CurrentCoord%next |
253 |
END IF !IE.NE.0 |
254 |
if (I3.NE.0) THEN |
255 |
|
256 |
NbAngles=NbAngles+1 |
257 |
! values of the coordinate (bond) is calculated for the reference geometry |
258 |
! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
259 |
! determined using all geometries. |
260 |
Call vecteur(i2,I3,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
261 |
Call vecteur(i2,I1,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
262 |
! CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
263 |
WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Angle",NbAngles," angle,val_zmat", & |
264 |
angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),val_zmat(I,2) |
265 |
CurrentCoord%value=val_zmat(I,2)*Pi/180. |
266 |
CurrentCoord%SignDihedral = 0 |
267 |
CurrentCoord%At1=I1 |
268 |
CurrentCoord%At2=I2 |
269 |
CurrentCoord%At3=I3 |
270 |
CurrentCoord%Type="ANGLE" |
271 |
IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
272 |
ALLOCATE(CurrentCoord%next) |
273 |
NULLIFY(CurrentCoord%next%next) |
274 |
END IF |
275 |
CurrentCoord => CurrentCoord%next |
276 |
END IF ! matches IF (I3.NE.0) |
277 |
if (I4.NE.0) THEN |
278 |
NbDihedrals=NbDihedrals+1 |
279 |
Call vecteur(I2,I1,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
280 |
Call vecteur(I2,I3,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
281 |
Call vecteur(I3,I4,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3) |
282 |
CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, & |
283 |
vx4,vy4,vz4,norm4) |
284 |
CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, & |
285 |
vx5,vy5,vz5,norm5) |
286 |
! CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
287 |
! vx2,vy2,vz2,norm2)*Pi/180. |
288 |
WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Dihedral",NbDihedrals, & |
289 |
" angle_d,val_zmat", & |
290 |
angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
291 |
vx2,vy2,vz2,norm2) & |
292 |
,val_zmat(I,3) |
293 |
|
294 |
CurrentCoord%value = val_zmat(I,3)*Pi/180. |
295 |
CurrentCoord%SignDihedral = int(sign(1.D0,val_zmat(I,3))) |
296 |
CurrentCoord%At1=I1 |
297 |
CurrentCoord%At2=I2 |
298 |
CurrentCoord%At3=I3 |
299 |
CurrentCoord%At4=I4 |
300 |
CurrentCoord%Type="DIHEDRAL" |
301 |
IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
302 |
ALLOCATE(CurrentCoord%next) |
303 |
NULLIFY(CurrentCoord%next%next) |
304 |
END IF |
305 |
CurrentCoord => CurrentCoord%next |
306 |
END IF ! matches IF (I4.NE.0) |
307 |
END DO ! matches DO I=1,Nat |
308 |
deallocate(indzmat0) |
309 |
|
310 |
ELSE !! if (debugPFL) |
311 |
|
312 |
DO IGeom=1, NGeomI |
313 |
!IGeom=1 ! need when using only first geometry. |
314 |
x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |
315 |
y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |
316 |
z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |
317 |
|
318 |
! First step: we calculate the connectivity |
319 |
Liaisons=0 |
320 |
|
321 |
Call CalcCnct(Nat,atome,x,y,z,LIAISONS,r_cov,fact) |
322 |
|
323 |
IF (debug) THEN |
324 |
WRITE(*,*) 'Connectivity done' |
325 |
DO I=1,Nat |
326 |
WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,J),J=0,NMaxL) |
327 |
WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,:)) |
328 |
END DO |
329 |
DO I=1,Nat |
330 |
WRITE(*,'(1X,I5,":",I3," r_cov=",F15.6)') I,atome(I),r_cov(atome(I)) |
331 |
END DO |
332 |
END IF |
333 |
|
334 |
DO I=1,Nat |
335 |
IF (Liaisons(I,0).NE.0) THEN |
336 |
DO J=1,Liaisons(I,0) |
337 |
! We add the bond |
338 |
Iat=Liaisons(I,J) |
339 |
IF (Iat.GT.I) THEN |
340 |
! Checking the uniqueness of bond. |
341 |
! Duplicate bonds should not be added. |
342 |
AddPrimitiveCoord=.True. |
343 |
ScanCoord=>Coordinate |
344 |
DO WHILE (Associated(ScanCoord%next)) |
345 |
IF (ScanCoord%Type .EQ. 'BOND') THEN |
346 |
IF (ScanCoord%At1 .EQ. Iat) THEN |
347 |
IF (ScanCoord%At2 .EQ. I) THEN |
348 |
AddPrimitiveCoord=.False. |
349 |
END IF |
350 |
END IF |
351 |
END IF |
352 |
ScanCoord => ScanCoord%next |
353 |
END DO |
354 |
|
355 |
IF (AddPrimitiveCoord) THEN |
356 |
NbBonds=NbBonds+1 |
357 |
! values of the coordinate (bond) is calculated for the reference geometry |
358 |
! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) |
359 |
! are determined using all geometries. vecteur is defined in bib_oper.f |
360 |
Call vecteur(I,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
361 |
CurrentCoord%value=Norm2 |
362 |
CurrentCoord%SignDihedral = 0 |
363 |
CurrentCoord%At1=Iat |
364 |
CurrentCoord%At2=I |
365 |
CurrentCoord%Type="BOND" |
366 |
IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
367 |
ALLOCATE(CurrentCoord%next) |
368 |
NULLIFY(CurrentCoord%next%next) |
369 |
END IF |
370 |
CurrentCoord => CurrentCoord%next |
371 |
END IF ! matches IF (AddPrimitiveCoord) THEN |
372 |
END IF ! matches IF (Iat.GT.I) THEN |
373 |
|
374 |
! Call Annul(Liaisons,IAt,I,Nat,NMaxL) |
375 |
! Liaisons(Iat,0)=Liaisons(Iat,0)-1 |
376 |
|
377 |
! We add the valence angles |
378 |
DO K=J+1,Liaisons(I,0) |
379 |
Kat=Liaisons(I,K) |
380 |
IF (Kat.GT.I) THEN |
381 |
! Checking the uniqueness of valence angle. |
382 |
! Duplicate bonds should not be added. |
383 |
AddPrimitiveCoord=.True. |
384 |
ScanCoord=>Coordinate |
385 |
DO WHILE (Associated(ScanCoord%next)) |
386 |
IF (ScanCoord%Type .EQ. 'ANGLE') THEN |
387 |
IF (ScanCoord%At1 .EQ. Iat) THEN |
388 |
IF (ScanCoord%At2 .EQ. I) THEN |
389 |
IF (ScanCoord%At3 .EQ. Kat) THEN |
390 |
AddPrimitiveCoord=.False. |
391 |
END IF |
392 |
END IF |
393 |
END IF |
394 |
END IF |
395 |
ScanCoord => ScanCoord%next |
396 |
END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
397 |
|
398 |
IF (AddPrimitiveCoord) THEN |
399 |
IF (debug) WRITE(*,*) "Adding angle:",Iat,I,Kat |
400 |
NbAngles=NbAngles+1 |
401 |
! values of the coordinate (bond) is calculated for the reference geometry |
402 |
! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
403 |
! determined using all geometries. |
404 |
Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
405 |
Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
406 |
CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
407 |
CurrentCoord%SignDihedral = 0 |
408 |
sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.)) |
409 |
CurrentCoord%At1=Iat |
410 |
CurrentCoord%At2=I |
411 |
CurrentCoord%At3=KAt |
412 |
CurrentCoord%Type="ANGLE" |
413 |
IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
414 |
ALLOCATE(CurrentCoord%next) |
415 |
NULLIFY(CurrentCoord%next%next) |
416 |
END IF |
417 |
CurrentCoord => CurrentCoord%next |
418 |
END IF ! matches IF (AddPrimitiveCoord) THEN |
419 |
ELSE ! matches IF (Kat.GT.I) THEN |
420 |
! Kat is less than I. If there is a bond between I and Kat, then this angle |
421 |
! has already been taken into account. |
422 |
if (debug) WRITE(*,*) "Testing angle:",Iat,I,Kat |
423 |
Bond=.FALSE. |
424 |
DO L=1,Liaisons(Iat,0) |
425 |
IF (Liaisons(Iat,L).EQ.Kat) Bond=.TRUE. |
426 |
END DO |
427 |
IF (.NOT.Bond) THEN |
428 |
IF (debug) WRITE(*,*) "Not Bond Adding angle:",Iat,I,Kat |
429 |
|
430 |
! Checking the uniqueness of valence angle. |
431 |
! Duplicate bonds should not be added. |
432 |
AddPrimitiveCoord=.True. |
433 |
ScanCoord=>Coordinate |
434 |
DO WHILE (Associated(ScanCoord%next)) |
435 |
IF (ScanCoord%Type .EQ. 'ANGLE') THEN |
436 |
IF (ScanCoord%At1 .EQ. Iat) THEN |
437 |
IF (ScanCoord%At2 .EQ. I) THEN |
438 |
IF (ScanCoord%At3 .EQ. Kat) THEN |
439 |
AddPrimitiveCoord=.False. |
440 |
END IF |
441 |
END IF |
442 |
END IF |
443 |
END IF |
444 |
ScanCoord => ScanCoord%next |
445 |
END DO |
446 |
|
447 |
IF (AddPrimitiveCoord) THEN |
448 |
NbAngles=NbAngles+1 |
449 |
! values of the coordinate (bond) is calculated for the reference geometry |
450 |
! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) |
451 |
! are determined using all geometries. |
452 |
Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
453 |
Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
454 |
CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
455 |
CurrentCoord%SignDihedral = 0 |
456 |
CurrentCoord%At1=Iat |
457 |
CurrentCoord%At2=I |
458 |
CurrentCoord%At3=KAt |
459 |
CurrentCoord%Type="ANGLE" |
460 |
IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
461 |
ALLOCATE(CurrentCoord%next) |
462 |
NULLIFY(CurrentCoord%next%next) |
463 |
END IF |
464 |
CurrentCoord => CurrentCoord%next |
465 |
END IF ! matches IF (AddPrimitiveCoord) THEN |
466 |
END IF ! matches IF (.NOT.Bond) THEN |
467 |
END IF ! matches IF (Kat.GT.I) THEN, after else |
468 |
END DO ! matches DO K=J+1,Liaisons(I,0) |
469 |
END DO ! matches DO J=1,Liaisons(I,0) |
470 |
END IF ! matches IF (Liaisons(I,0).NE.0) THEN |
471 |
|
472 |
! We add dihedral angles. We do not concentrate on necessary angles, ie those that are |
473 |
! needed to build the molecule. Instead we collect all of them, and the choice of the best |
474 |
! ones will be done when diagonalizing the Baker matrix. |
475 |
IF (Liaisons(I,0).GE.3) Then |
476 |
DO J=1,Liaisons(I,0) |
477 |
Iat=Liaisons(I,J) |
478 |
! values of the coordinate (bond) is calculated for the reference geometry |
479 |
! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
480 |
! determined using all geometries. |
481 |
Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
482 |
|
483 |
DO K=J+1,Liaisons(I,0) |
484 |
Kat=Liaisons(I,K) |
485 |
Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
486 |
sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.)) |
487 |
|
488 |
DO L=K+1,Liaisons(I,0) |
489 |
Lat=Liaisons(I,L) |
490 |
if (debug) WRITE(*,*) "0-Dihe Kat,I,Iat:",Kat,I,Iat,angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2) |
491 |
Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3) |
492 |
sAngleIIatLat=abs(sin(angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2)*Pi/180.)) |
493 |
if (debug) WRITE(*,*) "0-Dihe I,Iat,Lat:",angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2) |
494 |
IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN |
495 |
if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat |
496 |
|
497 |
! Checking for the uniqueness of valence angle. |
498 |
! Duplicate bonds should not be added. |
499 |
AddPrimitiveCoord=.True. |
500 |
ScanCoord=>Coordinate |
501 |
DO WHILE (Associated(ScanCoord%next)) |
502 |
IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN |
503 |
IF (ScanCoord%At1 .EQ. Kat) THEN |
504 |
IF (ScanCoord%At2 .EQ. I) THEN |
505 |
IF (ScanCoord%At3 .EQ. Iat) THEN |
506 |
IF (ScanCoord%At4 .EQ. Lat) THEN |
507 |
AddPrimitiveCoord=.False. |
508 |
END IF |
509 |
END IF |
510 |
END IF |
511 |
END IF |
512 |
END IF |
513 |
ScanCoord => ScanCoord%next |
514 |
END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
515 |
|
516 |
! IF (AddPrimitiveCoord) THEN |
517 |
! NbDihedrals=NbDihedrals+1 |
518 |
! CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, & |
519 |
! vx5,vy5,vz5,norm5) |
520 |
!CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, & |
521 |
! vx4,vy4,vz4,norm4) |
522 |
! CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
523 |
! vx2,vy2,vz2,norm2)*Pi/180. |
524 |
! CurrentCoord%SignDihedral = 0 |
525 |
! CurrentCoord%At1=Kat |
526 |
! CurrentCoord%At2=I |
527 |
! CurrentCoord%At3=IAt |
528 |
! CurrentCoord%At4=LAt |
529 |
!WRITE(*,'(1X,":(dihed 1)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') & |
530 |
! ScanCoord%At1,ScanCoord%At2,ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
531 |
! CurrentCoord%Type="DIHEDRAL" |
532 |
! IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
533 |
! ALLOCATE(CurrentCoord%next) |
534 |
! NULLIFY(CurrentCoord%next%next) |
535 |
! END IF |
536 |
! CurrentCoord => CurrentCoord%next |
537 |
!END IF ! matches IF (AddPrimitiveCoord) THEN |
538 |
|
539 |
END IF ! matches IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN |
540 |
END DO ! matches DO L=K+1,Liaisons(I,0) |
541 |
END DO ! matches DO K=J+1,Liaisons(I,0) |
542 |
END DO ! matches DO J=1,Liaisons(I,0) |
543 |
END IF !matches IF (Liaisons(I,0).GE.3) Then |
544 |
|
545 |
! we now add other dihedrals |
546 |
DO J=1,Liaisons(I,0) |
547 |
Iat=Liaisons(I,J) |
548 |
If (Iat.LE.I) Cycle |
549 |
! values of the coordinate (bond) is calculated for the reference geometry |
550 |
! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
551 |
! determined using all geometries. |
552 |
Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
553 |
|
554 |
DO K=1,Liaisons(I,0) |
555 |
IF (Liaisons(I,K).EQ.Iat) Cycle |
556 |
Kat=Liaisons(I,K) |
557 |
Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
558 |
sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.)) |
559 |
|
560 |
DO L=1,Liaisons(Iat,0) |
561 |
Lat=Liaisons(Iat,L) |
562 |
IF (Lat.eq.I) Cycle |
563 |
|
564 |
if (debug) WRITE(*,*) "Dihe Kat,I,Iat:",angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),sAngleIatIKat |
565 |
Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx,vy,vz,Norm) |
566 |
|
567 |
sAngleIIatLat=abs(sin(angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2)*Pi/180.)) |
568 |
if (debug) WRITE(*,*) "Dihe I,Iat,Lat:",angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2) |
569 |
! we add the dihedral angle Kat-I-Iat-Lat |
570 |
if (debug) WRITE(*,*) "Trying dihedral:",Kat,I,Iat,Lat,sAngleIatIKat,sAngleIIatLat |
571 |
IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) & |
572 |
.AND.(sAngleIIatLat.ge.0.01d0)) THEN |
573 |
if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat |
574 |
|
575 |
! Checking for the uniqueness of valence angle. |
576 |
! Duplicate bonds should not be added. |
577 |
AddPrimitiveCoord=.True. |
578 |
ScanCoord=>Coordinate |
579 |
DO WHILE (Associated(ScanCoord%next)) |
580 |
IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN |
581 |
IF (ScanCoord%At1 .EQ. Kat) THEN |
582 |
IF (ScanCoord%At2 .EQ. I) THEN |
583 |
IF (ScanCoord%At3 .EQ. Iat) THEN |
584 |
IF (ScanCoord%At4 .EQ. Lat) THEN |
585 |
AddPrimitiveCoord=.False. |
586 |
END IF |
587 |
END IF |
588 |
END IF |
589 |
END IF |
590 |
END IF |
591 |
ScanCoord => ScanCoord%next |
592 |
END DO |
593 |
|
594 |
IF (AddPrimitiveCoord) THEN |
595 |
NbDihedrals=NbDihedrals+1 |
596 |
Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3) |
597 |
CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, & |
598 |
vx5,vy5,vz5,norm5) |
599 |
CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, & |
600 |
vx4,vy4,vz4,norm4) |
601 |
CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
602 |
vx2,vy2,vz2,norm2)*Pi/180. |
603 |
CurrentCoord%At1=Kat |
604 |
CurrentCoord%At2=I |
605 |
CurrentCoord%At3=IAt |
606 |
CurrentCoord%At4=LAt |
607 |
CurrentCoord%Type="DIHEDRAL" |
608 |
CurrentCoord%SignDihedral=int(sign(1.d0,angle_d(vx4,vy4,vz4,norm4, & |
609 |
vx5,vy5,vz5,norm5, vx2,vy2,vz2,norm2))) |
610 |
WRITE(*,'(1X,":(dihed 2)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') ScanCoord%At1,& |
611 |
ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
612 |
IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
613 |
ALLOCATE(CurrentCoord%next) |
614 |
NULLIFY(CurrentCoord%next%next) |
615 |
END IF |
616 |
CurrentCoord => CurrentCoord%next |
617 |
END IF ! matches IF (AddPrimitiveCoord) THEN |
618 |
END IF ! matches IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) ..... |
619 |
END DO ! matches DO L=1,Liaisons(Iat,0) |
620 |
END DO ! matches DO K=1,Liaisons(I,0) |
621 |
END DO ! matches DO J=1,Liaisons(I,0) |
622 |
END DO ! end of loop over all atoms, DO I=1, Nat |
623 |
END DO ! matches DO IGeom=1, NGeomI |
624 |
|
625 |
|
626 |
|
627 |
END IF ! if (debugPFL) |
628 |
DEALLOCATE(Liaisons) |
629 |
|
630 |
WRITE(*,*) "I have found" |
631 |
WRITE(*,*) NbBonds, " bonds" |
632 |
WRITE(*,*) NbAngles, " angles" |
633 |
WRITE(*,*) NbDihedrals, " dihedrals" |
634 |
|
635 |
WRITE(*,*) 'All in all...' |
636 |
ScanCoord=>Coordinate |
637 |
I=0 |
638 |
DO WHILE (Associated(ScanCoord%next)) |
639 |
I=I+1 |
640 |
SELECT CASE (ScanCoord%Type) |
641 |
CASE ('BOND') |
642 |
!WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value |
643 |
WRITE(*,'(1X,I3,":(bond )",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value |
644 |
CASE ('ANGLE') |
645 |
!WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, & |
646 |
! ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi |
647 |
WRITE(*,'(1X,I3,":(angle)",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, & |
648 |
ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi |
649 |
CASE ('DIHEDRAL') |
650 |
!WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,& |
651 |
! ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
652 |
WRITE(*,'(1X,I3,":(dihed)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,& |
653 |
ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
654 |
END SELECT |
655 |
ScanCoord => ScanCoord%next |
656 |
END DO |
657 |
|
658 |
! NPrim is the number of primitive coordinates. |
659 |
NPrim=NbBonds+NbAngles+NbDihedrals |
660 |
|
661 |
! Now calculating values of all primitive bonds for all geometries: |
662 |
ALLOCATE(Xprimitive(NgeomI,NPrim),XPrimitiveF(NGeomF,NPrim)) |
663 |
ALLOCATE(UMatF(NGeomF,NPrim,NCoord)) |
664 |
ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat)) |
665 |
|
666 |
! We calculate the Primitive values for the first geometry, this plays a special role |
667 |
! as it will serve as a reference for other geometries ! |
668 |
|
669 |
x(1:Nat) = XyzGeomI(1,1,1:Nat) |
670 |
y(1:Nat) = XyzGeomI(1,2,1:Nat) |
671 |
z(1:Nat) = XyzGeomI(1,3,1:Nat) |
672 |
|
673 |
Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(1,:)) |
674 |
|
675 |
DO IGeom=2, NGeomI |
676 |
|
677 |
x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |
678 |
y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |
679 |
z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |
680 |
|
681 |
Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(IGeom,:),Xprimitive(1,:)) |
682 |
|
683 |
! ScanCoord=>Coordinate |
684 |
! I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |
685 |
! DO WHILE (Associated(ScanCoord%next)) |
686 |
! I=I+1 |
687 |
! SELECT CASE (ScanCoord%Type) |
688 |
! CASE ('BOND') |
689 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) |
690 |
! Xprimitive(IGeom,I) = Norm2 |
691 |
! CASE ('ANGLE') |
692 |
! Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1) |
693 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) |
694 |
! Xprimitive(IGeom,I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
695 |
! CASE ('DIHEDRAL') |
696 |
|
697 |
|
698 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1) |
699 |
! Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2) |
700 |
! Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3) |
701 |
! Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, & |
702 |
! vx4,vy4,vz4,norm4) |
703 |
! Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, & |
704 |
! vx5,vy5,vz5,norm5) |
705 |
|
706 |
! DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
707 |
! vx2,vy2,vz2,norm2) |
708 |
! Xprimitive(IGeom,I) = DiheTmp*Pi/180. |
709 |
! ! We treat large dihedral angles differently as +180=-180 mathematically and physically |
710 |
! ! but this causes lots of troubles when converting baker to cart. |
711 |
! ! So we ensure that large dihedral angles always have the same sign |
712 |
! if (abs(ScanCoord%SignDihedral).NE.1) THEN |
713 |
! ScanCoord%SignDihedral=Int(Sign(1.D0,DiheTmp)) |
714 |
! ELSE |
715 |
! If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN |
716 |
! Xprimitive(IGeom,I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi |
717 |
! END IF |
718 |
! END IF |
719 |
! ! Another test... will all this be consistent ??? |
720 |
! ! We want the shortest path, so we check that the change in dihedral angles is less than Pi: |
721 |
! IF (IGeom.GT.1) THEN |
722 |
! IF (abs(Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN |
723 |
! if ((Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN |
724 |
! Xprimitive(IGeom,I)=Xprimitive(IGeom,I)-2.d0*Pi |
725 |
! ELSE |
726 |
! Xprimitive(IGeom,I)=Xprimitive(IGeom,I)+2.d0*Pi |
727 |
! END IF |
728 |
! END IF |
729 |
! END IF |
730 |
! END SELECT |
731 |
! ScanCoord => ScanCoord%next |
732 |
! END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
733 |
END DO ! matches DO IGeom=1, NGeomI |
734 |
|
735 |
IF (debug) THEN |
736 |
WRITE(*,*) "DBG Calc_Baker Xprimitives for all geometries" |
737 |
DO IGeom=1, NGeomI |
738 |
WRITE(*,*) "Geometry ",IGeom |
739 |
ScanCoord=>Coordinate |
740 |
I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |
741 |
DO WHILE (Associated(ScanCoord%next)) |
742 |
I=I+1 |
743 |
SELECT CASE (ScanCoord%Type) |
744 |
CASE ('BOND') |
745 |
WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive(IGeom,I) |
746 |
CASE ('ANGLE') |
747 |
WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, & |
748 |
ScanCoord%At2, ScanCoord%At3,Xprimitive(IGeom,I)*180./Pi |
749 |
CASE ('DIHEDRAL') |
750 |
WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,& |
751 |
ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive(IGeom,I)*180./Pi |
752 |
END SELECT |
753 |
ScanCoord => ScanCoord%next |
754 |
END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
755 |
END DO ! matches DO IGeom=1, NGeomI |
756 |
END IF |
757 |
! Here reference geometry is assigned to Geom. Earlier x,y,z(Nat) were assigned from XyzGeomI(...) |
758 |
! before the call of this (Calc_baker) subroutine and passed to this subroutine as arguments of the |
759 |
! subroutine. |
760 |
! PFL 19 Feb 2008 -> |
761 |
! In order to reproduce the B matrix in the Baker article, one has |
762 |
! to divide Geom by a0. |
763 |
! However, this is inconsitent with our way of storing geometries |
764 |
! so we do not do it ! |
765 |
|
766 |
Geom(1,:)=XyzGeomI(IGeomRef,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat) |
767 |
Geom(2,:)=XyzGeomI(IGeomRef,2,1:Nat) |
768 |
Geom(3,:)=XyzGeomI(IGeomRef,3,1:Nat) |
769 |
|
770 |
! <- PFL 19 Feb 2008 |
771 |
|
772 |
! We now have to compute dq/dx... |
773 |
ALLOCATE(BprimT(3*Nat,NPrim)) |
774 |
BprimT=0.d0 |
775 |
|
776 |
ScanCoord=>Coordinate |
777 |
I=0 |
778 |
DO WHILE (Associated(ScanCoord%next)) |
779 |
I=I+1 |
780 |
SELECT CASE (ScanCoord%Type) |
781 |
CASE ('BOND') |
782 |
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2, & |
783 |
Geom,BprimT(1,I)) |
784 |
CASE ('ANGLE') |
785 |
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2, & |
786 |
ScanCoord%At3,Geom,BprimT(1,I)) |
787 |
CASE ('DIHEDRAL') |
788 |
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2, & |
789 |
ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I)) |
790 |
! CALL CONSTRAINTS_TORSION_DER(Nat,ScanCoord%At1,ScanCoord%AT2, & |
791 |
! ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I)) |
792 |
|
793 |
END SELECT |
794 |
ScanCoord => ScanCoord%next |
795 |
END DO |
796 |
|
797 |
|
798 |
!!!!!!!!!!!!!!!!!!!! |
799 |
! |
800 |
! PFL Debug |
801 |
! |
802 |
!!!!!!!!!!!!!!!!!!!! |
803 |
|
804 |
if (debugPFL) THEN |
805 |
WRITe(*,*) "DBG PFL Calc_Baker =====" |
806 |
DzDc=0.d0 |
807 |
WRITE(*,*) "PFL DBG, DzdC" |
808 |
do I=1,Nat |
809 |
Geom(1,i)=XyzGeomI(IGeomRef,1,orderInv(i)) ! XyzGeomI(NGeomI,3,Nat) |
810 |
Geom(2,i)=XyzGeomI(IGeomRef,2,OrderInv(i)) |
811 |
Geom(3,i)=XyzGeomI(IGeomRef,3,OrderInv(i)) |
812 |
WRITE(*,*) I,OrderInv(I),Geom(:,I) |
813 |
end do |
814 |
CALL CalcBMat_int (nat, Geom, indzmat, dzdc, massat,atome) |
815 |
|
816 |
WRITE(*,*) "PFL DBG, BPrimT" |
817 |
DO I=1,NPrim |
818 |
WRITE(*,'(50(1X,F12.6))') BPrimT(:,I) |
819 |
END DO |
820 |
WRITe(*,*) "DBG PFL Calc_Baker =====" |
821 |
END IF |
822 |
|
823 |
!!!!!!!!!!!!!!!!!!!! |
824 |
! |
825 |
! PFL Debug |
826 |
! |
827 |
!!!!!!!!!!!!!!!!!!!! |
828 |
|
829 |
DEALLOCATE(Geom) |
830 |
|
831 |
! BprimT(3*Nat,NPrim) |
832 |
! We now compute G=B(BT) matrix |
833 |
ALLOCATE(Gmat(NPrim,NPrim)) |
834 |
GMat=0.d0 |
835 |
DO I=1,NPrim |
836 |
DO J=1,3*Nat |
837 |
GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0))) |
838 |
END DO |
839 |
END DO |
840 |
|
841 |
!!!!!!!!!!!!!!!!!!!! |
842 |
! |
843 |
! PFL Debug ---> |
844 |
! |
845 |
!!!!!!!!!!!!!!!!!!!! |
846 |
if (debugPFL) THEN |
847 |
GMat=0.d0 |
848 |
DO I=1,NPrim |
849 |
GMat(I,I)=1. |
850 |
END DO |
851 |
END IF |
852 |
|
853 |
!!!!!!!!!!!!!!!!!!!! |
854 |
! |
855 |
! <--- PFL Debug |
856 |
! |
857 |
!!!!!!!!!!!!!!!!!!!! |
858 |
|
859 |
!DO I=1,NPrim |
860 |
! WRITE(*,'(20(1X,F6.3))') GMat(1:min(20,NPrim),I) |
861 |
!END DO |
862 |
|
863 |
! We diagonalize G |
864 |
ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim)) |
865 |
Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim) |
866 |
Call Trie(NPrim,EigVal,EigVec,NPrim) |
867 |
DO I=1,NPrim |
868 |
WRITE(*,'(1X,"Vector ",I3,": e=",F8.3)') I,EigVal(i) |
869 |
WRITE(*,'(20(1X,F8.4))') EigVec(1:min(20,NPrim),I) |
870 |
END DO |
871 |
|
872 |
!DO I=1,NPrim |
873 |
! WRITE(*,'(1X,"Vector ",I3,(20(F8.3)))') I,EigVec(1:min(20,NPrim),I) |
874 |
! END DO |
875 |
|
876 |
!DO I=1,NPrim |
877 |
! WRITE(*,'(1X,"Vector ",I3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,& |
878 |
! F8.3,F8.3,F8.3,F8.3,F8.3)') I,EigVec(1,I),EigVec(4,I),EigVec(6,I),EigVec(14,I),EigVec(16,I),& |
879 |
!EigVec(2,I),EigVec(3,I),EigVec(5,I),EigVec(12,I),EigVec(13,I),EigVec(15,I),EigVec(7,I),& |
880 |
!EigVec(8,I),EigVec(9,I),EigVec(10,I),EigVec(11,I),EigVec(17,I) |
881 |
!END DO |
882 |
|
883 |
! We construct the new DzDc(b=baker) matrix |
884 |
! UMat is nonredundant vector set, i.e. set of eigenvectors of BB^T corresponding |
885 |
! to eigenvalues > zero. |
886 |
ALLOCATE(UMat(NPrim,NCoord)) |
887 |
ALLOCATE(UMat_local(NPrim,NCoord)) |
888 |
ALLOCATE(UMat_tmp(NPrim,3*Nat-6)) |
889 |
! BMat_BakerT(3*Nat,NCoord), allocated in Path.f90, |
890 |
! NCoord=3*Nat-6 |
891 |
BMat_BakerT = 0.d0 |
892 |
J=0 |
893 |
UMat=0.d0 |
894 |
UMat_tmp=0.d0 |
895 |
DO I=1,NPrim |
896 |
IF (abs(eigval(I))>=1e-6) THEN |
897 |
J=J+1 |
898 |
DO K=1,NPrim |
899 |
!DzDcb(:,J)=DzDcb(:,J)+Eigvec(K,I)*BprimT(:,K) ! original |
900 |
! BprimT is transpose of B^prim. |
901 |
! B = UMat^T B^prim, B^T = (B^prim)^T UMat |
902 |
BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I) |
903 |
!Dzdcb(:,J)=DzDcb(:,J)+Eigvec(K,:)*BprimT(I,K) |
904 |
END DO |
905 |
IF(J .GT. 3*Nat-6) THEN |
906 |
WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. 1e-6(=UMat) (=' & |
907 |
,J,') exceeded 3*Nat-6=',3*Nat-6, & |
908 |
'Stopping the calculation.' |
909 |
STOP |
910 |
END IF |
911 |
UMat_tmp(:,J) = Eigvec(:,I) |
912 |
END IF |
913 |
END DO |
914 |
|
915 |
if (debug) THEN |
916 |
WRITE(*,*) "DBG Calc_Baker: UMat" |
917 |
DO I=1,3*Nat-6 |
918 |
WRITE(*,'(50(1X,F12.8))') UMat_tmp(:,I) |
919 |
END DO |
920 |
END IF |
921 |
|
922 |
! THIS IS TO BE CORRECTED: |
923 |
UMat = UMat_tmp |
924 |
!J=0 |
925 |
!DO I=1, 3*Nat-6 |
926 |
! IF ((I .NE. 6) .AND. (I .NE. 7) .AND. (I .NE. 9)) Then |
927 |
! J=J+1 |
928 |
! Print *, 'J=', J, ' I=', I |
929 |
! UMat(:,J) = UMat_tmp(:,I) |
930 |
!END IF |
931 |
!END DO |
932 |
! BMat_BakerT=0.d0 |
933 |
!DO I=1,NCoord |
934 |
! DO J=1,NPrim |
935 |
! B = UMat^T B^prim, B^T = (B^prim)^T UMat |
936 |
! BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat(J,I) |
937 |
! END DO |
938 |
!END DO |
939 |
! TO BE CORRECTED, ENDS HERE. |
940 |
|
941 |
IntCoordI=0.d0 |
942 |
DO IGeom=1, NGeomI |
943 |
DO I=1, NPrim |
944 |
! Transpose of UMat is needed below, that is why UMat(I,:). |
945 |
IntCoordI(IGeom,:) = IntCoordI(IGeom,:) + UMat(I,:)*Xprimitive(IGeom,I) |
946 |
END DO |
947 |
END DO |
948 |
|
949 |
if (debug) THEN |
950 |
WRITE(*,*) "DBG Calc_Baker IntCoordI" |
951 |
DO IGeom=1, NGeomI |
952 |
WRITE(*,'(1X,I5,":",30(1X,F10.6))') IGeom,IntCoordI(IGeom,:) |
953 |
END DO |
954 |
END IF |
955 |
|
956 |
|
957 |
! the following might be used to get rid of some coordinates that |
958 |
! do not respect the symmetry of the problem. |
959 |
|
960 |
! ! We look at the group symmetry of the molecule. Nothing complicate here |
961 |
! ! we just look for planar symmetry. |
962 |
! ! We first place it into the standard orientation |
963 |
! Call Std_ori(Nat,x,y,z,a,b,c,massat) |
964 |
|
965 |
! ! Is the molecule planar |
966 |
! DO I=1,Nat |
967 |
! WRITE(*,*) nom(atome(i)),x(i),y(i),z(i) |
968 |
! END DO |
969 |
DEALLOCATE(BprimT,GMat,EigVal,EigVec) ! BprimT is B^prim^T |
970 |
! Calculation of BTransInv starts here: |
971 |
! Calculation of BBT(3*Nat-6,3*Nat-6)=BB^T: |
972 |
! BMat_BakerT(3*Nat,NCoord) is Transpose of B = UMat^TB^prim |
973 |
ALLOCATE(BBT(NCoord,NCoord)) |
974 |
BBT = 0.d0 |
975 |
DO I=1, NCoord |
976 |
DO J=1, 3*Nat |
977 |
! BBT(:,I) forms BB^T |
978 |
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) |
979 |
END DO |
980 |
END DO |
981 |
|
982 |
!ALLOCATE(V(3*Nat,3*Nat)) |
983 |
!Call GINVSE(BBT,3*Nat,3*Nat,BBT_inv,3*Nat,3*Nat,3*Nat,1e-6,V,3*Nat,FAIL,NOUT) |
984 |
!DEALLOCATE(V) |
985 |
ALLOCATE(BBT_inv(NCoord,NCoord)) |
986 |
! GenInv is in Mat_util.f90 |
987 |
Call GenInv(NCoord,BBT,BBT_inv,NCoord) |
988 |
|
989 |
! Calculation of (B^T)^-1 = (BB^T)^-1B: |
990 |
!ALLOCATE(BTransInv(3*Nat-6,3*Nat)) ! allocated in Path.f90 |
991 |
BTransInv = 0.d0 |
992 |
DO I=1, 3*Nat |
993 |
DO J=1, NCoord |
994 |
BTransInv(:,I) = BTransInv(:,I) + BBT_inv(:,J)*BMat_BakerT(I,J) |
995 |
END DO |
996 |
END DO |
997 |
BTransInv_local = BTransInv |
998 |
UMat_local = UMat |
999 |
DEALLOCATE(BBT,BBT_inv,UMat_tmp) |
1000 |
DEALLOCATE(x_refGeom,y_refGeom,z_refGeom,x,y,z) |
1001 |
IF (debug) WRITE(*,*) "DBG Calc_baker over." |
1002 |
END SUBROUTINE Calc_baker |