root / src / Extrapol_mixed.f90 @ 7
Historique  Voir  Annoter  Télécharger (6,42 ko)
1 
SUBROUTINE Extrapol_mixed(s,dist,x0,y0,z0,xgeom,Coef) 

2  
3 
! This subroutine constructs the path, and if dist<>Infinity, it samples 
4 
! the path to obtain geometries. 
5 
! Basically, you call it twice: i) dist=infinity, it will calculate the 
6 
! length of the path 
7 
! ii) dist finite, it will give you the images you want along the path. 
8 
! 
9 
! For now, it gives equidistant geometries 
10 
!!!!!!!!!!!!!!!! 
11 
! 
12 
! v2.0 
13 
! Uses Mix2cart for conversions. 
14  
15  
16 
use Path_module 
17 
use Io_module 
18  
19  
20 
IMPLICIT NONE 
21  
22  
23 
REAL(KREAL), INTENT(OUT) :: s 
24 
REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat) 
25 
REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord) 
26  
27 
INTEGER(KINT) :: IdxGeom, I, J, K, Idx 
28 
REAL(KREAL) :: Rmsd, MRot(3, 3), ds, u 
29  
30 
REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3) 
31 
REAL(KREAL), ALLOCATABLE :: IntCoordTmp(:),DerInt(:) ! (NCoord) 
32  
33  
34 
LOGICAL :: debug, print 
35 
LOGICAL, EXTERNAL :: valid 
36  
37 
!We will calculate the length of the path, in MW coordinates... 
38 
! this is done is a stupid way: we interpolate the zmatrix values, 
39 
! convert them into cartesian, weight the cartesian 
40 
! and calculate the evolution of the distance ! 
41 
! We have to follow the same procedure for every geometry 
42 
! so even for the first one, we have to convert it from zmat 
43 
! to cartesian ! 
44  
45 
debug=valid("pathcreate") 
46 
print=valid("printgeom") 
47  
48 
if (debug) Call Header("Entering Extrapol_mixed") 
49 
if (debug) WRITE(*,*) "NFroz,NCart",NFroz,NCart 
50 
ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),IntCoordTmp(NCoord),DerInt(NCoord)) 
51 
IdxGeom=1 
52 

53 
IntCoordTmp=IntCoordI(1,1:NCoord) 
54 
call Mixed2Cart(Nat,IndZmat,IntCoordTmp,XyzTmp2(1,1)) 
55  
56 
XyzTmp=XyzTmp2 
57  
58 
! We align this geometry with the original one 
59 
! PFL 01/feb/2007: useless when dealing with more than three cartesian atoms 
60 
! PFL 17/July/2006: only if we have more than 4 atoms. 
61 
IF ((NCart.LT.3).AND.(Nat.GT.4)) THEN 
62 
! IF (Nat.GT.4) THEN 
63 
Call CalcRmsd(Nat, x0,y0,z0, & 
64 
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & 
65 
MRot,rmsd,FRot,FAlign) 
66 
END IF 
67  
68 
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) 
69 
IntCoordF(IdxGeom,:)=IntCoordI(1,:) 
70  
71 
if (print.AND.(Dist.LE.1e20)) THEN 
72 
WRITE(IOOUT,'(1X,I5)') Nat 
73 
WRITE(IOOUT,*) "# Cartesian Coordinates for geom",IdxGeom 
74 
DO i=1,Nat 
75 
If (Renum) THEN 
76 
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & 
77 
(XyzTmp2(Order(I),J),J=1,3) 
78 
ELSE 
79 
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & 
80 
(XyzTmp2(I,J),J=1,3) 
81 
END IF 
82 
END DO 
83 
END IF 
84  
85 
! Calculate tangents for first geometry 
86 
u=0.d0 
87 
DO Idx=1,NCoord 
88 
call splintder(u,IntCoordTmp(Idx),DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx)) 
89 
END DO 
90 
IntTangent(IdxGeom,:)=DerInt 
91  
92  
93 
! First geometry already initialized 
94  
95 
s=0. 
96 
! valzmat=0.d0 
97  
98 
DO K=1,NMaxPtPath 
99 
u=real(K)/NMaxPtPath*(NGeomI1.) 
100  
101 
XYZTmp2=0. 
102  
103 
! We generate the interpolated coordinates 
104 
DO Idx=1,NCoord 
105 
call splintder(u,IntCoordTmp(Idx),DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx)) 
106 
END DO 
107  
108 
! We convert it into Cartesian coordinates 
109 
call Mixed2Cart(Nat,IndZmat,IntCoordTmp,XyzTmp2(1,1)) 
110  
111 
! We calculate ds 
112 
ds=0. 
113 
DO I=1,Nat 
114 
DO J=1,3 
115 
ds=ds+MassAt(I)*(XYZTMp2(I,J)XYZTmp(I,J))**2 
116 
XYZTmp(I,J)=XyzTMP2(I,J) 
117 
ENDDO 
118 
ENDDO 
119 
s=s+sqrt(ds) 
120  
121 
! if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist 
122 
if (s>=dist) THEN 
123 
s=sdist 
124 
IdxGeom=IdxGeom+1 
125 
SGeom(IdxGeom)=s+IdxGeom*dist 
126 
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) 
127  
128 
IntCoordF(IdxGeom,:)=IntCoordTmp 
129 
IntTangent(IdxGeom,:)=DerInt 
130  
131 
if (print) THEN 
132 
WRITE(IOOUT,'(1X,A,12(1X,F10.5))') "#Internal Coord",IntCoordTmp(1:NCoord) 
133 
WRITE(IOOUT,'(1X,A,12(1X,F10.5))') "#Internal tan",IntTangent(IdxGeom,1:NCoord) 
134 
WRITE(IOOUT,'(1X,I5)') Nat 
135 
WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K,u,Xgeom(NGeomI) 
136  
137 
! PFL 01/feb/2007: useless when dealing with more than three cartesian atoms 
138 
! PFL 17/July/2006: only if we have more than 4 atoms. 
139 
IF ((NCart.LT.3).AND.(Nat.GT.4)) THEN 
140 
! IF (Nat.GT.4) THEN 
141 
Call CalcRmsd(Nat, x0,y0,z0, & 
142 
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & 
143 
MRot,rmsd,FRot,FAlign) 
144 
END IF 
145  
146 
DO i=1,Nat 
147 
If (Renum) THEN 
148 
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & 
149 
(XyzTmp2(Order(I),J),J=1,3) 
150 
ELSE 
151 
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & 
152 
(XyzTmp2(I,J),J=1,3) 
153 
END IF 
154 
END DO 
155  
156 
END IF 
157 
END IF 
158 
END DO 
159  
160  
161 
if ((s>=0.1*dist).AND.(IdxGeom.EQ.NGeomF)) THEN 
162 
WRITE(*,*) "** PathCreate ***" 
163 
WRITE(*,*) "Distribution of points along the path is wrong." 
164 
WRITE(*,*) "Increase value of NMaxPtPath in the input file" 
165 
WRITE(*,*) "Present value is:", NMaxPtPath 
166 
STOP 
167 
END IF 
168  
169 
IdxGeom=NGeomF 
170  
171 
! We have to add the last geometry. We copy the last geom of Initial geometries. 
172 
IntCoordF(IdxGeom,:)=IntCoordI(NGeomI,:) 
173 
IntTangent(IdxGeom,:)=DerInt 
174  
175 
! we convert it to cartesian geom 
176 
call Mixed2Cart(Nat,IndZmat,IntCoordF(IdxGeom,1),XyzTmp2(1,1)) 
177 
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) 
178  
179 
if (print) THEN 
180 
WRITE(IOOUT,'(1X,I5)') Nat 
181 
WRITE(IOOUT,*) "# Cartesian coord for last Geometry ",IdxGeom,K,u,Xgeom(NGeomI) 
182 
! PFL 01/feb/2007: useless when dealing with more than three cartesian atoms 
183 
! PFL 17/July/2006: only if we have more than 4 atoms. 
184 
IF ((NCART.LT.3).AND.(Nat.GT.4)) THEN 
185 
! IF (Nat.GT.4) THEN 
186 
Call CalcRmsd(Nat, x0,y0,z0, & 
187 
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & 
188 
MRot,rmsd,FRot,FAlign) 
189 
END IF 
190  
191 
DO i=1,Nat 
192 
If (Renum) THEN 
193 
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & 
194 
(XyzTmp2(Order(I),J),J=1,3) 
195 
ELSE 
196 
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))), & 
197 
(XyzTmp2(I,J),J=1,3) 
198 
END IF 
199 
END DO 
200 
END IF 
201  
202  
203 
if (debug) WRITE(*,*) 's final =',s 
204  
205 
DEALLOCATE(XyzTmp,XyzTmp2) 
206  
207 
if (debug) Call Header("Extrapol_Mixed Over") 
208  
209 

210 
END SUBROUTINE EXTRAPOL_MIXED 