root / src / Calc_baker_allGeomF.f90 @ 8
History  View  Annotate  Download (7.3 kB)
1 
SUBROUTINE Calc_baker_allGeomF() 

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 : BMat_BakerT,Nat,NCoord,UMatF, & 
9 
NPrim,BTransInvF,Coordinate, & 
10 
ScanCoord,BprimT,BBT,BBT_inv,XprimitiveF, & 
11 
NgeomF,XyzGeomF 
12 
! BMat_BakerT(3*Nat,NCoord), NCoord=3*Nat or NFree=3*Nat6Symmetry_elimination 
13 
! depending upon the coordinate choice. IntCoordI(NGeomI,NCoord) where 
14 
! UMatF(NGeomI,NPrim,NCoord), NCoord number of vectors in UMat matrix, i.e. NCoord 
15 
! Baker coordinates. NPrim is the number of primitive internal coordinates. 
16 

17 
Use Io_module 
18 
IMPLICIT NONE 
19  
20 
REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat) 
21 
! NPrim is the number of primitive coordinates and NCoord is the number 
22 
! of internal coordinates. BMat is actually (NPrim,3*Nat). 
23 
REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim) 
24 
! EigVec(..) contains ALL eigevectors of BMat times BprimT, NOT only Baker Coordinate vectors. 
25 
REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim) 
26 
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) 
27 
REAL(KREAL), ALLOCATABLE :: XPrimRef(:) ! NPrim 
28 
INTEGER(KINT) :: IGeom 
29 

30  
31 
INTEGER(KINT) :: I, J, K 
32  
33  
34 
LOGICAL :: debug 
35 
LOGICAL :: DebugPFL 
36  
37 
INTERFACE 
38 
function valid(string) result (isValid) 
39 
CHARACTER(*), intent(in) :: string 
40 
logical :: isValid 
41 
END function VALID 
42  
43 
FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2) 
44 
use Path_module, only : Pi,KINT, KREAL 
45 
real(KREAL) :: v1x,v1y,v1z,norm1 
46 
real(KREAL) :: v2x,v2y,v2z,norm2 
47 
real(KREAL) :: angle 
48 
END FUNCTION ANGLE 
49  
50 
FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3) 
51 
use Path_module, only : Pi,KINT, KREAL 
52 
real(KREAL) :: v1x,v1y,v1z,norm1 
53 
real(KREAL) :: v2x,v2y,v2z,norm2 
54 
real(KREAL) :: v3x,v3y,v3z,norm3 
55 
real(KREAL) :: angle_d,ca,sa 
56 
END FUNCTION ANGLE_D 
57  
58  
59  
60 
SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef) 
61 
! 
62 
! This subroutine uses the description of a list of Coordinates 
63 
! to compute the values of the coordinates for a given geometry. 
64 
! 
65 
!!!!!!!!!! 
66 
! Input: 
67 
! Na: INTEGER, Number of atoms 
68 
! x,y,z(Na): REAL, cartesian coordinates of the considered geometry 
69 
! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes 
70 
! NPrim, INTEGER: Number of coordinates to compute 
71 
! 
72 
! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for 
73 
! a former geometry. Useful for Dihedral angles evolution... 
74  
75 
!!!!!!!!!!! 
76 
! Output: 
77 
! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates 
78 
! 
79 
!!!!!!!!! 
80  
81 
Use VarTypes 
82 
Use Io_module 
83 
Use Path_module, only : pi 
84  
85 
IMPLICIT NONE 
86  
87 
Type (ListCoord), POINTER :: Coordinate 
88 
INTEGER(KINT), INTENT(IN) :: Nat,NPrim 
89 
REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat) 
90 
REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
91 
REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim) 
92  
93 
END SUBROUTINE CALC_XPRIM 
94 
END INTERFACE 
95  
96  
97  
98  
99 
debug=valid("Calc_baker_allGeomF") 
100 
debugPFL=valid("bakerPFL") 
101 
if (debug) WRITE(*,*) '============ Entering Calc_baker_allGeomF =============' 
102  
103 
ALLOCATE(Geom(3,Nat),x(Nat),y(Nat),z(Nat)) 
104 
ALLOCATE(XPrimRef(NPrim)) 
105 

106 
! Now calculating values of all primitive bonds for all final geometries: 
107 
DO IGeom=1, NGeomF 
108 
x(1:Nat) = XyzGeomF(IGeom,1,1:Nat) 
109 
y(1:Nat) = XyzGeomF(IGeom,2,1:Nat) 
110 
z(1:Nat) = XyzGeomF(IGeom,3,1:Nat) 
111 
XPrimREf=XPrimitiveF(IGeom,:) 
112 
Call Calc_XPrim(nat,x,y,z,Coordinate,NPrim,XPrimitiveF(IGeom,:),XPrimRef) 
113 
END DO ! matches DO IGeom=1, NGeomF 
114 

115 
ALLOCATE(BprimT(3*Nat,NPrim)) 
116 
ALLOCATE(Gmat(NPrim,NPrim)) 
117 
ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim)) 
118 
ALLOCATE(BBT(NCoord,NCoord)) 
119 
ALLOCATE(BBT_inv(NCoord,NCoord)) 
120 
BTransInvF = 0.d0 
121  
122 
DO IGeom=1, NGeomF 
123 
Geom(1,:)=XyzGeomF(IGeom,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat) 
124 
Geom(2,:)=XyzGeomF(IGeom,2,1:Nat) 
125 
Geom(3,:)=XyzGeomF(IGeom,3,1:Nat) 
126  
127 
BprimT=0.d0 
128 
ScanCoord=>Coordinate 
129 
I=0 
130 
DO WHILE (Associated(ScanCoord%next)) 
131 
I=I+1 
132 
SELECT CASE (ScanCoord%Type) 
133 
CASE ('BOND') 
134 
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2, & 
135 
Geom,BprimT(1,I)) 
136 
CASE ('ANGLE') 
137 
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2, & 
138 
ScanCoord%At3,Geom,BprimT(1,I)) 
139 
CASE ('DIHEDRAL') 
140 
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2, & 
141 
ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I)) 
142 
END SELECT 
143 
ScanCoord => ScanCoord%next 
144 
END DO 
145  
146 
! BprimT(3*Nat,NPrim) 
147 
! We now compute G=B(BT) matrix 
148 
GMat=0.d0 
149 
DO I=1,NPrim 
150 
DO J=1,3*Nat 
151 
GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0))) 
152 
END DO 
153 
END DO 
154 

155 
! Diagonalize G 
156 
EigVal=0.d0 
157 
EigVec=0.d0 
158 
Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim) 
159 
Call Trie(NPrim,EigVal,EigVec,NPrim) 
160 
DO I=1,NPrim 
161 
!WRITE(*,'(1X,"Vector ",I3,": e=",F8.3)') I,EigVal(i) 
162 
!WRITE(*,'(20(1X,F8.4))') EigVec(1:min(20,NPrim),I) 
163 
END DO 
164  
165 
! UMatF is nonredundant vector set, i.e. set of eigenvectors of BB^T 
166 
! corresponding to eigenvalues > zero. 
167 
! BMat_BakerT(3*Nat,NCoord), allocated in Path.f90, 
168 
! NCoord=3*Nat6 
169 
BMat_BakerT = 0.d0 
170 
J=0 
171 
DO I=1,NPrim 
172 
IF (abs(eigval(I))>=1e6) THEN 
173 
J=J+1 
174 
DO K=1,NPrim 
175 
! BprimT is transpose of B^prim. 
176 
! B = UMatF^T B^prim, B^T = (B^prim)^T UMatF 
177 
BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I) 
178 
END DO 
179 
IF(J .GT. 3*Nat6) THEN 
180 
WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. 1e6(=UMatF) (=' & 
181 
,J,') exceeded 3*Nat6=',3*Nat6, & 
182 
'Stopping the calculation.' 
183 
STOP 
184 
END IF 
185 
UMatF(IGeom,:,J) = Eigvec(:,I) 
186 
END IF 
187 
END DO 
188  
189 
!!!!!!!!!!!!!!!!!!!! 
190 
! 
191 
! Debug purposes 
192 
! 
193 
if (debugPFL) THEN 
194 
UMatF(IGeom,:,:)=0. 
195 
DO J=1,3*Nat6 
196 
UMatF(IGeom,J,J)=1. 
197 
END DO 
198 
END IF 
199  
200 

201 
!DO I=1, NPrim ! This loop is not needed because we already have IntCoordF 
202 
! from interpolation. 
203 
! Transpose of UMatF is needed below, that is why UMatF(IGeom,I,:). 
204 
! IntCoordF(IGeom,:) = IntCoordF(IGeom,:) + UMat(IGeom,I,:)*XprimitiveF(IGeom,I) 
205 
!END DO 
206  
207 
! Calculation of BTransInvF starts here: 
208 
! Calculation of BBT(3*Nat6,3*Nat6)=BB^T: 
209 
! BMat_BakerT(3*Nat,NCoord) is Transpose of B = UMatF^TB^prim 
210 

211 
BBT = 0.d0 
212 
DO I=1, NCoord 
213 
DO J=1, 3*Nat 
214 
! BBT(:,I) forms BB^T 
215 
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) 
216 
END DO 
217 
END DO 
218 

219 
Call GenInv(NCoord,BBT,BBT_inv,NCoord) ! GenInv is in Mat_util.f90 
220  
221 
! Calculation of (B^T)^1 = (BB^T)^1B: 
222 
DO I=1, 3*Nat 
223 
DO J=1, NCoord 
224 
BTransInvF(IGeom,:,I) = BTransInvF(IGeom,:,I) + BBT_inv(:,J)*BMat_BakerT(I,J) 
225 
END DO 
226 
END DO 
227 

228 
END DO !matches DO IGeom=1, NGeomF 
229 

230 
DEALLOCATE(BBT,BBT_inv,BprimT,GMat,EigVal,EigVec) 
231 
DEALLOCATE(Geom,x,y,z,XprimRef) 
232 

233 
IF (debug) WRITE(*,*) "DBG Calc_baker_allGeomF over." 
234 
END SUBROUTINE Calc_baker_allGeomF 