root / src / Calc_Xprim.f90 @ 7
History  View  Annotate  Download (6.7 kB)
1 
SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef) 

2 
! 
3 
! This subroutine uses the description of a list of Coordinates 
4 
! to compute the values of the coordinates for a given geometry. 
5 
! 
6 
!!!!!!!!!! 
7 
! Input: 
8 
! Na: INTEGER, Number of atoms 
9 
! x,y,z(Na): REAL, cartesian coordinates of the considered geometry 
10 
! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes 
11 
! NPrim, INTEGER: Number of coordinates to compute 
12 
! 
13 
! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for 
14 
! a former geometry. Useful for Dihedral angles evolution... 
15  
16 
!!!!!!!!!!! 
17 
! Output: 
18 
! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates 
19 
! 
20 
!!!!!!!!! 
21  
22 
Use VarTypes 
23 
Use Io_module 
24 
Use Path_module, only : pi 
25  
26 
IMPLICIT NONE 
27  
28 
Type (ListCoord), POINTER :: Coordinate 
29 
INTEGER(KINT), INTENT(IN) :: Nat,NPrim 
30 
REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat) 
31 
REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
32 
REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim) 
33  
34  
35 
Type (ListCoord), POINTER :: ScanCoord 
36  
37 
real(KREAL) :: vx1,vy1,vz1,norm1 
38 
real(KREAL) :: vx2,vy2,vz2,norm2 
39 
real(KREAL) :: vx3,vy3,vz3,norm3 
40 
real(KREAL) :: vx4,vy4,vz4,norm4 
41 
real(KREAL) :: vx5,vy5,vz5,norm5 
42  
43 
INTEGER(KINT) :: I 
44  
45 
REAL(KREAL) :: DiheTmp 
46  
47 
LOGICAL :: debug, debugPFL 
48  
49 
INTERFACE 
50 
function valid(string) result (isValid) 
51 
CHARACTER(*), intent(in) :: string 
52 
logical :: isValid 
53 
END function VALID 
54  
55 
FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2) 
56 
use Path_module, only : Pi,KINT, KREAL 
57 
real(KREAL) :: v1x,v1y,v1z,norm1 
58 
real(KREAL) :: v2x,v2y,v2z,norm2 
59 
real(KREAL) :: angle 
60 
END FUNCTION ANGLE 
61  
62 
FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2) 
63 
use Path_module, only : Pi,KINT, KREAL 
64 
real(KREAL) :: v1x,v1y,v1z,norm1 
65 
real(KREAL) :: v2x,v2y,v2z,norm2 
66 
real(KREAL) :: SinAngle 
67 
END FUNCTION SINANGLE 
68  
69  
70 
FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3) 
71 
use Path_module, only : Pi,KINT, KREAL 
72 
real(KREAL) :: v1x,v1y,v1z,norm1 
73 
real(KREAL) :: v2x,v2y,v2z,norm2 
74 
real(KREAL) :: v3x,v3y,v3z,norm3 
75 
real(KREAL) :: angle_d,ca,sa 
76 
END FUNCTION ANGLE_D 
77  
78 
END INTERFACE 
79  
80  
81 
debug=valid("Calc_Xprim") 
82 
debugPFL=valid("BakerPFL") 
83  
84 
if (debug) Call Header("Entering Cal_XPrim") 
85  
86  
87 
IF (debug) THEN 
88 
WRITE(*,*) "DBG Calc_Xprim, geom" 
89 
DO I=1,Nat 
90 
WRITE(*,'(1X,I5,3(1X,F15.3))') I,X(I),y(i),z(i) 
91 
END DO 
92 
IF (Present(XPrimRef)) THEN 
93 
WRITE(*,*) "XPrimRef" 
94 
WRITE(*,'(15(1X,F10.6))') XPrimRef 
95 
END IF 
96 
WRITE(*,*) "NPrim:",NPrim 
97 
END IF 
98  
99  
100 
WRITE(*,*) "Coordinate:",associated(Coordinate),coordinate%Type 
101  
102 
ScanCoord => Coordinate 
103  
104 
WRITE(*,*) "ScanCoord:",associated(ScanCoord),ScanCoord%Type 
105  
106 
! WRITE(*,*) "coucou" 
107 
I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). 
108 
DO WHILE (Associated(ScanCoord%next)) 
109 
I=I+1 
110 
! WRITE(*,*) i 
111 
SELECT CASE (ScanCoord%Type) 
112 
CASE ('BOND') 
113 
Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) 
114 
Xprimitive(I) = Norm2 
115 
CASE ('ANGLE') 
116 
Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1) 
117 
Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) 
118 
Xprimitive(I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. 
119 
CASE ('DIHEDRAL') 
120 
Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1) 
121 
Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2) 
122 
Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3) 
123 
Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, & 
124 
vx4,vy4,vz4,norm4) 
125 
Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, & 
126 
vx5,vy5,vz5,norm5) 
127 

128 
DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & 
129 
vx2,vy2,vz2,norm2) 
130 
Xprimitive(I) = DiheTmp*Pi/180. 
131 
! PFL 15th March 2008 > 
132 
! I think that the test on changes less than Pi should be enough 
133 
!! We treat large dihedral angles differently as +180=180 mathematically and physically 
134 
!! but this causes lots of troubles when converting baker to cart. 
135 
!! So we ensure that large dihedral angles always have the same sign 
136 
! if (abs(ScanCoord%SignDihedral).NE.1) THEN 
137 
! ScanCoord%SignDihedral=Int(Sign(1.D0,DiheTmp)) 
138 
! ELSE 
139 
! If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN 
140 
! Xprimitive(I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi 
141 
! END IF 
142 
! END IF 
143 
!!!! < PFL 15th March 2008 
144 
! Another test... will all this be consistent ??? 
145 
! We want the shortest path, so we check that the change in dihedral angles is less than Pi: 
146 
IF (Present(XPrimRef)) THEN 
147 
IF (abs(Xprimitive(I)XPrimRef(I)).GE.Pi) THEN 
148 
if (debug) THEN 
149 
WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I) 
150 
WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi 
151 
END IF 
152 
if ((Xprimitive(I)XPrimRef(I)).GE.Pi) THEN 
153 
Xprimitive(I)=Xprimitive(I)2.d0*Pi 
154 
ELSE 
155 
Xprimitive(I)=Xprimitive(I)+2.d0*Pi 
156 
END IF 
157 
END IF 
158 
if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi 
159 
END IF 
160 
CASE DEFAULT 
161 
WRITE(*,*) "Scancoord%type unknown for I",I,scancoord%type 
162 
END SELECT 
163 
ScanCoord => ScanCoord%next 
164 
END DO ! matches DO WHILE (Associated(ScanCoord%next)) 
165  
166  
167  
168 
IF (debug) THEN 
169 
WRITE(*,*) "DBG Calc_Xprim Values" 
170 
WRITE(*,*) "Found ",I," primitives" 
171 

172 
ScanCoord=>Coordinate 
173 
I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). 
174 
DO WHILE (Associated(ScanCoord%next)) 
175 
I=I+1 
176 
SELECT CASE (ScanCoord%Type) 
177 
CASE ('BOND') 
178 
WRITE(*,'(1X,I3,":",I5,"  ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive(I) 
179 
CASE ('ANGLE') 
180 
WRITE(*,'(1X,I3,":",I5,"  ",I5,"  ",I5," = ",F15.6)') I,ScanCoord%At1, & 
181 
ScanCoord%At2, ScanCoord%At3,Xprimitive(I)*180./Pi 
182 
CASE ('DIHEDRAL') 
183 
WRITE(*,'(1X,I3,":",I5,"  ",I5,"  ",I5,"  ",I5," = ",F15.6)') I,ScanCoord%At1,& 
184 
ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive(I)*180./Pi 
185 
END SELECT 
186 
ScanCoord => ScanCoord%next 
187 
END DO ! matches DO WHILE (Associated(ScanCoord%next)) 
188 
END IF 
189  
190 
if (debug) Call Header(" Cal_XPrim OVER ") 
191  
192 
END SUBROUTINE Calc_Xprim 