## root / src / CalcDist_cart.f90 @ 7

History | View | Annotate | Download (12.3 kB)

1 |
SUBROUTINE CalcDist_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef) |
---|---|

2 | |

3 |
! This subroutine computes the curvilinear distance of the images on the path |

4 |
! It is adapted from Extrapol_cart |

5 | |

6 |
use Path_module |

7 |
use Io_module |

8 | |

9 | |

10 |
IMPLICIT NONE |

11 | |

12 | |

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

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

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

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

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

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

19 |
INTEGER(KINT) :: IdxGeom, I, J, K, Idx,Iat, IGeom |

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

21 |
INTEGER(KINT) :: NSpline |

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

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

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

25 |
REAL(KREAL) :: a_val, d |

26 | |

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

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

29 | |

30 | |

31 |
LOGICAL :: debug, print, printspline |

32 |
LOGICAL, EXTERNAL :: valid |

33 | |

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

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

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

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

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

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

40 |
! to cartesian ! |

41 | |

42 |
debug=valid("pathcreate") |

43 |
print=valid("printgeom") |

44 | |

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

46 | |

47 |
if (debug) Call Header("Entering CalcDist_cart") |

48 | |

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

50 |
NSpline=int(NMaxPtPath/100) |

51 | |

52 |
if (printspline) THEN |

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

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

55 |
OPEN(IOTMP,FILE=FileSpline) |

56 | |

57 |
END IF |

58 | |

59 | |

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

61 | |

62 |
! XyzTmp2=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/)) |

63 | |

64 |
if (debug) THEN |

65 |
WRITE(*,*) "DBG Extrapol_cart Initial geometries" |

66 |
DO IGeom=1,NGeomI |

67 |
WRITE(*,*) 'XyzGeomI, IGeom=',IGeom |

68 |
DO I=1,Nat |

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

70 |
(XyzGeomI(IGeom,J,I),J=1,3) |

71 |
END DO |

72 |
END DO |

73 |
END IF |

74 | |

75 | |

76 |
! In order to mesure only the relevant distance between two points |

77 |
! we align all geometries on the original one |

78 | |

79 |
DO IGeom=1,NGeomI |

80 | |

81 |
XyzTmp2=Reshape(XyzGeomI(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)) |

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

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

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

85 |
! PFL 24 Nov 2008 -> |

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

87 |
! PFL 8 Feb 2011 -> |

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

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

90 |
IF (Align) THEN |

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

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

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

94 |
FrozAtoms,MRot) |

95 |
ELSE |

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

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

98 |
MRot,rmsd,.TRUE.,.TRUE.) |

99 |
END IF |

100 |
! <- PFL 24 Nov 2008 |

101 |
END IF |

102 |
! -> PFL 8 Feb 2011 |

103 |
! END IF |

104 | |

105 |
XyzGeomI(IGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/)) |

106 |
END DO |

107 | |

108 |
if (print) THEN |

109 |
WRITE(*,*) "Aligned geometries" |

110 |
DO J=1, NGeomI |

111 |
WRITE(IOOUT,*) Nat |

112 |
WRITE(IOOUT,*) " Aligned geometry ",J,"/",NGeomI |

113 |
DO i=1,Nat |

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

115 |
XyzGeomI(J,1:3,I) |

116 |
END DO |

117 |
END DO |

118 |
END IF |

119 | |

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

121 | |

122 | |

123 |
! We initialize the first geometry |

124 | |

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

126 | |

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

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

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

130 |
! PFL 24 Nov 2008 -> |

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

132 |
! PFL 8 Feb 2011 -> |

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

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

135 |
IF (Align) THEN |

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

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

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

139 |
FrozAtoms,MRot) |

140 |
ELSE |

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

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

143 |
MRot,rmsd,.TRUE.,.TRUE.) |

144 |
END IF |

145 |
! <- PFL 24 Nov 2008 |

146 |
END IF |

147 |
! -> PFL 8 Feb 2011 |

148 |
! END IF |

149 | |

150 | |

151 | |

152 |
s=0.d0 |

153 |
SGeom(1)=0.d0 |

154 | |

155 |
if (printspline) THEN |

156 |
u=0.d0 |

157 |
DO Iat=1,Nat |

158 |
! We generate the interpolated coordinates |

159 |
if (Linear) THEN |

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

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

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

163 | |

164 |
ELSE |

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

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

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

168 |
END IF |

169 |
END DO |

170 | |

171 |
WRITE(IOTMP,*) Nat |

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

173 |
DO Iat=1,Nat |

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

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

176 |
END DO |

177 |
END IF |

178 | |

179 |
IdxGeom=1 |

180 | |

181 |
DO K=1,NMaxPtPath |

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

183 | |

184 |
DO Iat=1,Nat |

185 |
! We generate the interpolated coordinates |

186 |
if (Linear) THEN |

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

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

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

190 | |

191 |
ELSE |

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

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

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

195 |
END IF |

196 |
END DO |

197 | |

198 | |

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

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

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

202 |
! PFL 24 Nov 2008 -> |

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

204 |
! PFL 8 Feb 2011 -> |

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

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

207 |
IF (Align) THEN |

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

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

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

211 |
FrozAtoms,MRot) |

212 |
ELSE |

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

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

215 |
MRot,rmsd,.TRUE.,.TRUE.) |

216 |
END IF |

217 |
! <- PFL 24 Nov 2008 |

218 |
END IF |

219 |
! -> PFL 8 Feb 2011 |

220 |
! END IF |

221 | |

222 | |

223 | |

224 |
ds=0. |

225 |
DO I=1,Nat |

226 |
DO J=1,3 |

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

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

229 |
ENDDO |

230 |
ENDDO |

231 |
s=s+sqrt(ds) |

232 | |

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

234 | |

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

236 |
WRITE(IOTMP,*) Nat |

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

238 |
DO Iat=1,Nat |

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

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

241 |
END DO |

242 |
END IF |

243 | |

244 | |

245 |
if (s>=dist) THEN |

246 | |

247 |
if (debug) THEN |

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

249 |
DO i=1,Nat |

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

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

252 |
END DO |

253 |
END IF |

254 | |

255 |
s=s-dist |

256 |
IdxGeom=IdxGeom+1 |

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

258 | |

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

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

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

262 |
! PFL 24 Nov 2008 -> |

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

264 |
! PFL 8 Feb 2011 -> |

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

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

267 |
IF (Align) THEN |

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

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

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

271 |
FrozAtoms,MRot) |

272 |
ELSE |

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

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

275 |
MRot,rmsd,.TRUE.,.TRUE.) |

276 |
END IF |

277 |
! <- PFL 24 Nov 2008 |

278 |
END IF |

279 |
! -> PFL 8 Feb 2011 |

280 |
! END IF |

281 | |

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

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

284 | |

285 |
if (print) THEN |

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

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

288 | |

289 |
DO i=1,Nat |

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

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

292 |
END DO |

293 | |

294 |
END IF |

295 |
END IF ! s>= dist |

296 |
ENDDO ! K |

297 | |

298 | |

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

300 |
s=s-dist |

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

302 | |

303 |
if (printspline) THEN |

304 |
WRITE(IOTMP,*) Nat |

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

306 |
DO Iat=1,Nat |

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

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

309 |
END DO |

310 |
END IF |

311 | |

312 | |

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

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

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

316 |
! PFL 24 Nov 2008 -> |

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

318 |
! PFL 8 Feb 2011 -> |

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

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

321 |
IF (Align) THEN |

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

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

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

325 |
FrozAtoms,MRot) |

326 |
ELSE |

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

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

329 |
MRot,rmsd,.TRUE.,.TRUE.) |

330 |
END IF |

331 |
! <- PFL 24 Nov 2008 |

332 |
END IF |

333 |
! -> PFL 8 Feb 2011 |

334 |
! END IF |

335 | |

336 |
IdxGeom=IdxGeom+1 |

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

338 | |

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

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

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

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

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

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

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

346 |
STOP |

347 |
END IF |

348 | |

349 | |

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

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

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

353 |
! PFL 24 Nov 2008 -> |

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

355 |
! PFL 8 Feb 2011 -> |

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

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

358 |
IF (Align) THEN |

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

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

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

362 |
FrozAtoms,MRot) |

363 |
ELSE |

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

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

366 |
MRot,rmsd,.TRUE.,.TRUE.) |

367 |
END IF |

368 |
! <- PFL 24 Nov 2008 |

369 |
END IF |

370 |
! -> PFL 8 Feb 2011 |

371 |
! END IF |

372 | |

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

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

375 | |

376 |
if (print) THEN |

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

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

379 | |

380 |
DO i=1,Nat |

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

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

383 |
END DO |

384 |
END IF |

385 |
END IF |

386 | |

387 |
DEALLOCATE(XyzTmp,XyzTmp2) |

388 | |

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

390 | |

391 |
if (printspline) CLOSE(IOTMP) |

392 | |

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

394 | |

395 |
END SUBROUTINE EXTRAPOL_CART |