## root / src / Extrapol_cart.f90 @ 5

History | View | Annotate | Download (10.7 kB)

1 |
SUBROUTINE Extrapol_cart(iopt,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 length of the path |

6 |
! ii) dist finite, it will give you the images you want along the path. |

7 |
! |

8 |
! For now, it gives equidistant geometries |

9 |
! |

10 | |

11 |
use Path_module |

12 |
use Io_module |

13 | |

14 | |

15 |
IMPLICIT NONE |

16 | |

17 | |

18 |
REAL(KREAL), INTENT(OUT) :: s |

19 |
! X0(Nat),Y0(Nat),Z0(Nat): reference geometry. |

20 |
REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat) |

21 |
REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord) |

22 |
! Iopt: Number of the cycles for the optimization |

23 |
INTEGER(KINT), INTENT(IN) :: Iopt |

24 |
INTEGER(KINT) :: IdxGeom, I, J, K, Iat |

25 |
! NSpline is the number of points along the interpolating path |

26 |
INTEGER(KINT) :: NSpline |

27 |
! FileSpline: Filename to save the interpolating path coordinates |

28 |
CHARACTER(LCHARS) :: FileSpline,TmpChar |

29 |
REAL(KREAL) :: Rmsd, MRot(3, 3), ds, u |

30 | |

31 |
REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3) |

32 |
REAL(KREAL), ALLOCATABLE :: DerXyz(:,:) ! Nat,3 |

33 | |

34 | |

35 |
LOGICAL :: debug, print, printspline |

36 |
LOGICAL, EXTERNAL :: valid |

37 | |

38 |
!We will calculate the length of the path, in MW coordinates... |

39 |
! this is done is a stupid way: we interpolate the zmatrix values, |

40 |
! convert them into cartesian, weight the cartesian |

41 |
! and calculate the evolution of the distance ! |

42 |
! We have to follow the same procedure for every geometry |

43 |
! so even for the first one, we have to convert it from zmat |

44 |
! to cartesian ! |

45 | |

46 |
debug=valid("pathcreate") |

47 |
print=valid("printgeom") |

48 | |

49 |
printspline=(valid("printspline").AND.(dist<=1e30)) |

50 | |

51 |
if (debug) Call Header("Entering Extrapol_cart") |

52 | |

53 |
! We want 100 points along the interpolating path |

54 |
NSpline=int(NMaxPtPath/100) |

55 | |

56 |
if (printspline) THEN |

57 |
WRITE(TmpChar,'(I5)') Iopt |

58 |
FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar)) |

59 |
OPEN(IOTMP,FILE=FileSpline) |

60 | |

61 |
END IF |

62 | |

63 | |

64 |
ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),DerXyz(Nat,3)) |

65 | |

66 |
! PFL 2011 June 22 |

67 |
! There was an alignment procedure there. |

68 |
! It has been moved to PathCreate before the |

69 |
! computation of the Spline coefficient. |

70 | |

71 |
if (debug) THEN |

72 |
WRITE(*,*) "Extrapol_cart- L72 - geometries " |

73 |
DO I=1,NGeomI |

74 |
WRITE(*,*) "Geom ",I |

75 |
DO J=1,Nat |

76 |
If (renum) THEN |

77 |
Iat=Order(J) |

78 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,Iat) |

79 |
ELSE |

80 |
Iat=OrderInv(J) |

81 |
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),XyzGeomI(I,1:3,J) |

82 |
END IF |

83 |
END DO |

84 |
END DO |

85 |
END IF |

86 | |

87 | |

88 |
XyzGeomF(1,:,:)=XyzGeomI(1,:,:) |

89 | |

90 | |

91 |
! We initialize the first geometry |

92 | |

93 |
XyzTmp=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/)) |

94 | |

95 | |

96 | |

97 |
s=0.d0 |

98 |
SGeom(1)=0.d0 |

99 | |

100 |
if (printspline) THEN |

101 |
u=0.d0 |

102 |
DO Iat=1,Nat |

103 |
! We generate the interpolated coordinates |

104 |
if (Linear) THEN |

105 |
call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat)) |

106 |
call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat)) |

107 |
call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat)) |

108 | |

109 |
ELSE |

110 |
call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2)) |

111 |
call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1)) |

112 |
call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat)) |

113 |
END IF |

114 |
END DO |

115 | |

116 |
WRITE(IOTMP,*) Nat |

117 |
WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s |

118 |
DO Iat=1,Nat |

119 |
WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)), & |

120 |
(XyzTmp2(Iat,1:3)) |

121 |
END DO |

122 |
END IF |

123 | |

124 |
IdxGeom=1 |

125 | |

126 |
DO K=1,NMaxPtPath |

127 |
u=real(K)/NMaxPtPath*(NGeomI-1.) |

128 | |

129 |
DO Iat=1,Nat |

130 |
! We generate the interpolated coordinates |

131 |
if (Linear) THEN |

132 |
call LinearInt(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat)) |

133 |
call LinearInt(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat)) |

134 |
call LinearInt(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat)) |

135 | |

136 |
ELSE |

137 |
call splintder(u,XyzTmp2(Iat,1),DerXyz(Iat,1),NGeomI,xgeom(1),XyzGeomI(1,1,Iat),Coef(1,3*Iat-2)) |

138 |
call splintder(u,XyzTmp2(Iat,2),DerXyz(Iat,2),NGeomI,xgeom(1),XyzGeomI(1,2,Iat),Coef(1,3*Iat-1)) |

139 |
call splintder(u,XyzTmp2(Iat,3),DerXyz(Iat,3),NGeomI,xgeom(1),XyzGeomI(1,3,Iat),Coef(1,3*Iat)) |

140 |
END IF |

141 |
END DO |

142 | |

143 | |

144 |
! PFL 2011 June 22: We now align on the previous one |

145 |
! as we do for internal coord interpolation |

146 | |

147 |
! We align this geometry with the previous one |

148 | |

149 |
! PFL 17/July/2006: only if we have more than 4 atoms. |

150 |
! IF (Nat.GT.4) THEN |

151 |
! PFL 24 Nov 2008 -> |

152 |
! If we have frozen atoms we align only those ones. |

153 |
! PFL 8 Feb 2011 -> |

154 |
! I add a flag to see if we should align or not. |

155 |
! For small systems, it might be better to let the user align himself |

156 |
IF (Align) THEN |

157 |
if (NFroz.GT.0) THEN |

158 |
Call AlignPartial(Nat, & |

159 |
xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3), & |

160 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

161 |
FrozAtoms,MRot) |

162 |
ELSE |

163 |
Call CalcRmsd(Nat, & |

164 |
xyzTmp(1,1),xyzTmp(1,2),xyzTMP(1,3), & |

165 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

166 |
MRot,rmsd,.TRUE.,.TRUE.) |

167 |
END IF |

168 |
! <- PFL 24 Nov 2008 |

169 |
END IF |

170 |
! -> PFL 8 Feb 2011 |

171 |
! END IF |

172 | |

173 | |

174 | |

175 |
ds=0. |

176 |
DO I=1,Nat |

177 |
DO J=1,3 |

178 |
ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2 |

179 |
XYZTmp(I,J)=XyzTMP2(I,J) |

180 |
ENDDO |

181 |
ENDDO |

182 |
s=s+sqrt(ds) |

183 | |

184 |
! if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist |

185 | |

186 |
if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN |

187 |
WRITE(IOTMP,*) Nat |

188 |
WRITE(IOTMP,'(1X,A,1X,F15.6)') "Debug Spline - Coord=Cart,s=",s |

189 |
DO Iat=1,Nat |

190 |
WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)), & |

191 |
(XyzTmp2(Iat,1:3)) |

192 |
END DO |

193 |
END IF |

194 | |

195 | |

196 |
if (s>=dist) THEN |

197 | |

198 |
if (debug) THEN |

199 |
WRITE(*,*) "DBG Interpol_cart",s |

200 |
DO i=1,Nat |

201 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |

202 |
(XyzTmp2(I,J),J=1,3) |

203 |
END DO |

204 |
END IF |

205 | |

206 |
s=s-dist |

207 |
IdxGeom=IdxGeom+1 |

208 |
SGeom(IdxGeom)=s+dist*(IdxGeom-1) |

209 | |

210 |
! We align this geometry with the original one |

211 |
! PFL 17/July/2006: only if we have more than 4 atoms. |

212 |
! IF (Nat.GT.4) THEN |

213 |
! PFL 24 Nov 2008 -> |

214 |
! If we have frozen atoms we align only those ones. |

215 |
! PFL 8 Feb 2011 -> |

216 |
! I add a flag to see if we should align or not. |

217 |
! For small systems, it might be better to let the user align himself |

218 |
IF (Align) THEN |

219 |
if (NFroz.GT.0) THEN |

220 |
Call AlignPartial(Nat,x0,y0,z0, & |

221 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

222 |
FrozAtoms,MRot) |

223 |
ELSE |

224 |
Call CalcRmsd(Nat, x0,y0,z0, & |

225 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

226 |
MRot,rmsd,.TRUE.,.TRUE.) |

227 |
END IF |

228 |
! <- PFL 24 Nov 2008 |

229 |
END IF |

230 |
! -> PFL 8 Feb 2011 |

231 |
! END IF |

232 | |

233 |
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) |

234 |
XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/)) |

235 | |

236 |
if (print) THEN |

237 |
WRITE(IOOUT,'(1X,I5)') Nat |

238 |
WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K |

239 | |

240 |
DO i=1,Nat |

241 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |

242 |
(XyzTmp2(I,J),J=1,3) |

243 |
END DO |

244 | |

245 |
END IF |

246 |
END IF ! s>= dist |

247 |
ENDDO ! K |

248 | |

249 | |

250 |
if (s>=0.9*dist) THEN |

251 |
s=s-dist |

252 |
XyzTmp2=Reshape(XyzGeomI(NGeomI,:,:),(/Nat,3/),ORDER=(/2,1/)) |

253 | |

254 |
if (printspline) THEN |

255 |
WRITE(IOTMP,*) Nat |

256 |
WRITE(IOTMP,*) "Debug Spline - Coord=Cart,s=",s |

257 |
DO Iat=1,Nat |

258 |
WRITE(IOTMP,'(1X,A2,3(1X,F15.6))') Nom(Atome(Iat)), & |

259 |
(XyzTmp2(Iat,1:3)) |

260 |
END DO |

261 |
END IF |

262 | |

263 | |

264 |
! We align this geometry with the original one |

265 |
! PFL 17/July/2006: only if we have more than 4 atoms. |

266 |
! IF (Nat.GT.4) THEN |

267 |
! PFL 24 Nov 2008 -> |

268 |
! If we have frozen atoms we align only those ones. |

269 |
! PFL 8 Feb 2011 -> |

270 |
! I add a flag to see if we should align or not. |

271 |
! For small systems, it might be better to let the user align himself |

272 |
IF (Align) THEN |

273 |
if (NFroz.GT.0) THEN |

274 |
Call AlignPartial(Nat,x0,y0,z0, & |

275 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

276 |
FrozAtoms,MRot) |

277 |
ELSE |

278 |
Call CalcRmsd(Nat, x0,y0,z0, & |

279 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

280 |
MRot,rmsd,.TRUE.,.TRUE.) |

281 |
END IF |

282 |
! <- PFL 24 Nov 2008 |

283 |
END IF |

284 |
! -> PFL 8 Feb 2011 |

285 |
! END IF |

286 | |

287 |
IdxGeom=IdxGeom+1 |

288 |
SGeom(IdxGeom)=s+dist*(IdxGeom-1) |

289 | |

290 |
IF (IdxGeom.GT.NGeomF) THEN |

291 |
WRITE(IOOUT,*) "!!! ERROR in Extrapol_cart !!!!" |

292 |
WRITE(IOOUT,*) "Too many structures. Increase NMaxPath" |

293 |
WRITE(*,*) "** PathCreate ***" |

294 |
WRITE(*,*) "Distribution of points along the path is wrong." |

295 |
WRITE(*,*) "Increase value of NMaxPtPath in the input file" |

296 |
WRITE(*,*) "Present value is:", NMaxPtPath |

297 |
STOP |

298 |
END IF |

299 | |

300 | |

301 |
! We align this geometry with the original one |

302 |
! PFL 17/July/2006: only if we have more than 4 atoms. |

303 |
! IF (Nat.GT.4) THEN |

304 |
! PFL 24 Nov 2008 -> |

305 |
! If we have frozen atoms we align only those ones. |

306 |
! PFL 8 Feb 2011 -> |

307 |
! I add a flag to see if we should align or not. |

308 |
! For small systems, it might be better to let the user align himself |

309 |
IF (Align) THEN |

310 |
if (NFroz.GT.0) THEN |

311 |
Call AlignPartial(Nat,x0,y0,z0, & |

312 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

313 |
FrozAtoms,MRot) |

314 |
ELSE |

315 |
Call CalcRmsd(Nat, x0,y0,z0, & |

316 |
xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3), & |

317 |
MRot,rmsd,.TRUE.,.TRUE.) |

318 |
END IF |

319 |
! <- PFL 24 Nov 2008 |

320 |
END IF |

321 |
! -> PFL 8 Feb 2011 |

322 |
! END IF |

323 | |

324 |
XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) |

325 |
XyzTangent(IdxGeom,:)=Reshape(DerXyz(:,:),(/3*Nat/)) |

326 | |

327 |
if (print) THEN |

328 |
WRITE(IOOUT,'(1X,I5)') Nat |

329 |
WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K |

330 | |

331 |
DO i=1,Nat |

332 |
WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)), & |

333 |
(XyzTmp2(I,J),J=1,3) |

334 |
END DO |

335 |
END IF |

336 |
END IF |

337 | |

338 |
DEALLOCATE(XyzTmp,XyzTmp2) |

339 | |

340 |
if (debug) WRITE(*,*) 's final =',s |

341 | |

342 |
if (printspline) CLOSE(IOTMP) |

343 | |

344 |
if (debug) Call Header("Extrapol_cart over") |

345 | |

346 |
END SUBROUTINE EXTRAPOL_CART |