## root / src / PathCreate.f90 @ 5

History | View | Annotate | Download (39.5 kB)

1 |
!================================================================ |
---|---|

2 |
! |

3 |
! This subroutine create a path using starting geometries |

4 |
! it is 'original' in the sense that the path is generated |

5 |
! by interpolating internal coordinates instead of cartesian. |

6 |
! Comes from the Extrapol.f subroutine of Cart suite. |

7 |
! Rewritten in F90 to be more flexible |

8 |
! |

9 |
!================================================================ |

10 | |

11 |
SUBROUTINE PathCreate(IOpt) |

12 | |

13 |
use Io_module |

14 |
use Path_module, only : Nat, NGeomI, NCoord, NGeomF, IGeomRef, NMaxL, IReparam, IReparamT, & |

15 |
Coord, Frozen, Cart, NCart, NFroz, XyzGeomI, atome, r_cov, fact, & |

16 |
IndZmat, Renum, Order, OrderInv, IntCoordI, IntCoordF,Pi, Nom, & |

17 |
ISpline, IntTangent, NPrim,Xprimitive_t, OptGeom, & |

18 |
UMatF, UMat_local, XyzTangent, Linear, Align, FrozAtoms,AtName |

19 |
! BMat_BakerT (3*Nat,NCoord), XyzGeomI(NGeomI,3,Nat), IGeomRef=-1 (default value) |

20 |
! IntCoordI(NGeomI,NCoord) |

21 | |

22 |
IMPLICIT NONE |

23 | |

24 |
REAL(KREAL), PARAMETER :: Inf=1e32 |

25 |
INTEGER(KINT), INTENT(IN) :: IOpt |

26 | |

27 |
REAL(KREAL), ALLOCATABLE :: val_zmat(:,:), val_zmatTmp(:,:) ! (3,Nat) |

28 |
REAL(KREAL), ALLOCATABLE :: X0(:), Y0(:), Z0(:) ! Nat |

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

30 |
! REAL(KREAL), ALLOCATABLE :: XyzTmp3(:,:) ! (3,Nat) |

31 |
REAL(KREAL), ALLOCATABLE :: IntCoordTmp(:) ! (NCoord) |

32 |
INTEGER(KINT), ALLOCATABLE :: IndZmat0(:,:) !Nat,5 |

33 |
REAL(KREAL), ALLOCATABLE :: Coef(:,:) ! (NGeomI,NCoord) |

34 |
! NPrim=Number of primitive coordinates. |

35 |
REAL(KREAL), ALLOCATABLE :: XGeom(:) ! NGeomI |

36 |
REAL(KREAL), ALLOCATABLE :: XGeomF(:) ! NGeomF |

37 |
INTEGER(KINT) :: NGeomS |

38 | |

39 |
REAL(KREAL) :: Rmsd, MRot(3, 3), Delta |

40 |
REAL(KREAL), PARAMETER :: TwoPi=360.d0 |

41 |
REAL(KREAL) ::s, dist, a_val, d |

42 |
INTEGER(KINT) :: I, J, igeom, iat |

43 |
INTEGER(KINT) :: Itmp, Jat |

44 |
INTEGER(KINT) :: Idx,IBeg |

45 |
CHARACTER(LCHARS) :: Title |

46 |
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) |

47 |
LOGICAL :: debug, print |

48 |
LOGICAL, SAVE :: First=.TRUE. |

49 |
! Frozen contains the indices of frozen atoms |

50 |
LOGICAL, ALLOCATABLE :: FCart(:) ! Nat |

51 |
INTEGER, ALLOCATABLE :: AtCart(:) !Nat |

52 |
INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:),LiaisonsIni(:,:) ! (Nat,0:NMaxL) |

53 |
INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !nat |

54 |
INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !nat,nat |

55 |
INTEGER(KINT) :: NbFrag,NFragRef |

56 | |

57 | |

58 |
INTERFACE |

59 |
function valid(string) result (isValid) |

60 |
CHARACTER(*), intent(in) :: string |

61 |
logical :: isValid |

62 |
END function VALID |

63 | |

64 |
FUNCTION test_zmat(na,ind_zmat) |

65 | |

66 |
USE Path_module |

67 | |

68 |
logical :: test_zmat |

69 |
integer(KINT) :: na |

70 |
integer(KINT) :: ind_zmat(Na,5) |

71 |
END FUNCTION TEST_ZMAT |

72 | |

73 | |

74 |
SUBROUTINE die(routine, msg, file, line, unit) |

75 | |

76 |
Use VarTypes |

77 |
Use io_module |

78 | |

79 |
implicit none |

80 | |

81 |
character(len=*), intent(in) :: routine, msg |

82 |
character(len=*), intent(in), optional :: file |

83 |
integer(KINT), intent(in), optional :: line, unit |

84 | |

85 |
END SUBROUTINE die |

86 | |

87 |
SUBROUTINE Warning(routine, msg, file, line, unit) |

88 | |

89 |
Use VarTypes |

90 |
Use io_module |

91 | |

92 |
implicit none |

93 | |

94 |
character(len=*), intent(in) :: routine, msg |

95 |
character(len=*), intent(in), optional :: file |

96 |
integer(KINT), intent(in), optional :: line, unit |

97 | |

98 |
END SUBROUTINE Warning |

99 | |

100 |
END INTERFACE |

101 | |

102 | |

103 | |

104 |
debug=valid("pathcreate") |

105 |
print=valid("printgeom") |

106 |
if (debug) Call header("Entering PathCreate") |

107 |
if (debug.AND.First) WRITE(*,*) "= First call =" |

108 | |

109 |
if (debug) WRITE(*,*) "Iopt,IReparam=",Iopt,IReparam |

110 | |

111 |
ALLOCATE(Coef(NGeomI,NCoord),val_zmat(Nat,3),val_zmatTMP(Nat,3)) |

112 |
ALLOCATE(X0(Nat),Y0(Nat),Z0(Nat),Indzmat0(Nat,5)) |

113 |
ALLOCATE(X(Nat),Y(Nat),Z(Nat),Xgeom(NGeomI),XgeomF(NGeomF)) |

114 |
! ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),XyzTmp3(3,Nat)) |

115 |
ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3)) |

116 |
ALLOCATE(IntCoordTmp(NCoord)) |

117 | |

118 |
Do I=1,NGeomI |

119 |
XGeom(I)=FLoat(I)-1.d0 |

120 |
if (Print) THEN |

121 |
WRITE(*,*) "PathCreate - L121 - Initial geometries " |

122 |
DO J=1,Nat |

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

124 |
END DO |

125 |
END IF |

126 |
END DO |

127 | |

128 |
! First iteration of the optimization: |

129 |
IF (First) THEN ! matches at L484 |

130 |
if (debug) WRITE(*,*) "first iteration in PathCreate.90, L93" |

131 | |

132 |
! This is the first call, so we might have to generate a Zmat |

133 |
First=.FALSE. |

134 | |

135 |
IF ( ((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD.EQ."MIXED") & |

136 |
.OR.(COORD.EQ."BAKER")).AND.(IGeomRef.LE.0)) THEN |

137 | |

138 |
! User has not defined the Reference geometry, we have to choose it ourselves! |

139 |
! This is done by counting the number of fragments of each geometry. The |

140 |
! reference geometry is the one with the least fragments. When there are |

141 |
! frozen or cart atoms, they are discarded to count the fragments. |

142 | |

143 |
ALLOCATE(Liaisons(Nat,0:NMaxL)) |

144 |
ALLOCATE(Fragment(nat),NbAtFrag(nat),FragAt(nat,nat),FCart(Nat)) |

145 |
FCart=.FALSE. |

146 |
! FCart(Nat) was under IF below earlier but is being accessed after this loop. |

147 |
IF ((Frozen(1).NE.0).OR.(Cart(1).NE.0)) THEN |

148 |
ALLOCATE(AtCart(Nat),LiaisonsIni(Nat,0:NMaxL)) |

149 |
FCart=.FALSE. |

150 |
NCart=0 |

151 |
Idx=1 |

152 |
DO WHILE (Frozen(Idx).NE.0) |

153 |
IF (Frozen(Idx).LT.0) THEN |

154 |
DO I=Frozen(Idx-1),abs(Frozen(Idx)) |

155 |
IF (.NOT.FCart(I)) THEN |

156 |
NCart=NCart+1 |

157 |
AtCart(NCart)=I |

158 |
FCart(I)=.TRUE. |

159 |
END IF |

160 |
END DO |

161 |
ELSE IF (.NOT.FCart(Frozen(Idx))) THEN |

162 |
FCart(Frozen(Idx))=.TRUE. |

163 |
NCart=NCart+1 |

164 |
AtCart(NCart)=Frozen(Idx) |

165 |
END IF |

166 |
Idx=Idx+1 |

167 |
END DO ! matches DO WHILE (Frozen(Idx).NE.0) |

168 |
NFroz=NCart |

169 |
Idx=1 |

170 |
DO WHILE (Cart(Idx).NE.0) |

171 |
IF (Cart(Idx).LT.0) THEN |

172 |
DO I=Cart(Idx-1),abs(Cart(Idx)) |

173 |
IF (.NOT.FCart(I)) THEN |

174 |
NCart=NCart+1 |

175 |
AtCart(NCart)=I |

176 |
FCart(I)=.TRUE. |

177 |
END IF |

178 |
END DO |

179 |
ELSEIF (.NOT.FCart(Cart(Idx))) THEN |

180 |
FCart(Cart(Idx))=.TRUE. |

181 |
NCart=NCart+1 |

182 |
AtCart(NCart)=Cart(Idx) |

183 |
END IF |

184 |
Idx=Idx+1 |

185 |
END DO !matches DO WHILE (Cart(Idx).NE.0) |

186 |
IF (debug) THEN |

187 |
WRITE(*,'(1X,A,I4,A,I4,A)') "DBG PathCreate - GeomRef :: Found ", & |

188 |
NFroz," frozen atoms, and a total of ",NCart, & |

189 |
" atoms described in cartesian coordinates" |

190 |
WRITE(*,*) "AtCart:",AtCart(1:NCart) |

191 |
END IF |

192 |
END IF ! IF ((Frozen(1).NE.0).OR.(Cart(1).NE.0)) THEN |

193 | |

194 |
NFragRef=2*Nat |

195 |
DO IGeom=1,NGeomI |

196 |
x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |

197 |
y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |

198 |
z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |

199 |
! we calculate the connectivities |

200 |
Call CalcCnct(nat,atome,x,y,z,LIAISONS,r_cov,fact) |

201 |
! if needed, we get rid of links between cart or frozen atoms |

202 |
IF (NCart.GE.2) THEN |

203 |
LiaisonsIni=Liaisons |

204 |
DO I=1,NCart |

205 |
Iat=AtCart(I) |

206 |
if (debug) WRITE(*,'(20(1X,I3))') I,Iat, & |

207 |
(LiaisonsIni(Iat,Itmp),ITmp=0,LiaisonsIni(Iat,0)) |

208 |
DO J=1,LiaisonsIni(Iat,0) |

209 |
Jat=LiaisonsIni(Iat,J) |

210 |
IF (FCart(Jat)) THEN |

211 |
Call Annul(Liaisons,Iat,Jat) |

212 |
Call Annul(Liaisons,Jat,Iat) |

213 |
Call Annul(LiaisonsIni,Jat,Iat) |

214 |
Liaisons(Iat,0)=Liaisons(Iat,0)-1 |

215 |
Liaisons(Jat,0)=Liaisons(Jat,0)-1 |

216 |
LiaisonsIni(Jat,0)=LiaisonsIni(Jat,0)-1 |

217 |
END IF |

218 |
END DO |

219 |
END DO |

220 |
END IF ! matches IF (NCart.GE.2) THEN |

221 |
! we now calculate the number of fragments. |

222 |
Call Decomp_frag(nat,liaisons,.NOT.Fcart,nbfrag,Fragment,NbAtFrag,FragAt) |

223 |
IF (debug) THEN |

224 |
WRITE(*,*) 'Debug PathCreat Line190' |

225 |
WRITE(*,*) 'NbFrag, NbFragRef=',NbFrag,NFragRef |

226 |
END IF |

227 | |

228 |
IF (NbFrag.LT.NFragRef) THEN |

229 |
NFragRef=NbFrag |

230 |
! The reference geometry, IGeomRef, is determined based on the least number |

231 |
! of fragments (here if (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD |

232 |
!.EQ."MIXED")).AND.(IGeomRef.LE.0))), otherwise it still has the value of |

233 |
! IGeomRef=-1. |

234 |
IGeomRef=IGeom |

235 |
END IF |

236 |
END DO ! matches DO IGeom=1,NGeomI |

237 |
DEALLOCATE(FCart) |

238 |
END IF ! matches IF (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD.EQ."MIXED") |

239 |
! .OR.(COORD.EQ."BAKER")).AND.(IGeomRef.LE.0)). |

240 | |

241 |
IF ((COORD.EQ."CART").AND.(IGeomREF.LE.0)) THEN |

242 |
IGeomRef=1 |

243 |
CALL Warning('PathCreate L209','IGeomRef<=0',UNIT=IOOUT) |

244 |
END IF |

245 | |

246 |
if (debug) WRITE(*,*) "DBG PathCreate : IGeomRef= ",IGeomRef |

247 | |

248 |
! we now compute the internal coordinates for this geometry ! |

249 |
! The reference geometry, IGeomRef, is determined based on the least number |

250 |
! of fragments if (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD |

251 |
!.EQ."MIXED")).AND.(IGeomRef.LE.0)), otherwise it has the value of |

252 |
! IGeomRef=-1. |

253 |
x(1:Nat) = XyzGeomI(IGeomRef,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat) |

254 |
y(1:Nat) = XyzGeomI(IGeomRef,2,1:Nat) |

255 |
z(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat) |

256 | |

257 |
IndZmat=0 |

258 |
IndZmat0=0 |

259 | |

260 |
SELECT CASE(COORD) |

261 |
CASE ('BAKER') |

262 |
! atome(Nat), r_cov(0:Max_Z) |

263 |
! Frozen(Nat) contains the indices of frozen atoms |

264 |
Call Calc_baker(atome,r_cov,fact,frozen,IGeomRef) |

265 |
if (debug) THEN |

266 |
WRITE(*,*) "UMat_local after Calc_baker" |

267 |
DO I=1,3*Nat-6 |

268 |
WRITE(*,'(50(1X,F12.8))') UMat_local(:,I) |

269 |
END DO |

270 |
END IF |

271 |
DO IGeom=1,NGeomF |

272 |
UMatF(IGeom,:,:)=UMat_Local(:,:) |

273 |
END DO |

274 | |

275 |
ALLOCATE(Xprimitive_t(NPrim)) |

276 |
! x, y, z corresponds the reference geometry. |

277 |
DO I=1,Nat |

278 |
Order(I)=I |

279 |
OrderInv(I)=I |

280 |
END DO |

281 |
x0=x |

282 |
y0=y |

283 |
z0=z |

284 |
CASE ('ZMAT','HYBRID') |

285 |
! IN case we are using Hybrid or zmat coordinate system, we have to |

286 |
! generate the internal coordinates with a Z-Matrix |

287 |
! IN case we have frozen atoms, we generate a constrained Zmat |

288 |
IF (Frozen(1).NE.0) THEN |

289 |
Renum=.True. |

290 |
call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact,frozen) |

291 |
ELSE |

292 |
!no zmat... we create our own |

293 |
call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact) |

294 |
if (valid('old_zmat')) THEN |

295 |
call Calc_zmat(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact) |

296 |
WRITE(*,*) "DBG PathCreate Calc_Zmat.... STOP" |

297 |
STOP |

298 |
END IF |

299 |
END IF |

300 | |

301 |
if (debug) WRITE(*,*) 'Back from Calc_zmat' |

302 |
IF (.NOT.test_zmat(nat,IndZmat0)) THEN |

303 |
WRITE(IOOUT,*) "I cannot generate a valid Zmat... help" |

304 |
STOP |

305 |
END IF |

306 | |

307 |
! The IndZmat0 generated by Calc_zmat refers to the order of the input x,y,z |

308 |
! However, we need it with the new order so that we can use the zmat |

309 |
! coordinates to generate cartesian coordinates. |

310 |
! So, we have to reorder it... |

311 |
DO I=1,Nat |

312 |
Order(IndZmat0(I,1))=I |

313 |
OrderInv(I)=IndZmat0(I,1) |

314 |
! I let X0 being the Cartesian coordinates of the reference geometry. Before it was |

315 |
! those of the First geometry (which was also the reference one). This might change! |

316 |
X0(i)=X(IndZmat0(i,1)) |

317 |
Y0(i)=Y(IndZmat0(i,1)) |

318 |
Z0(i)=Z(IndZmat0(i,1)) |

319 |
END DO |

320 |
IndZmat(1,1)=Order(IndZmat0(1,1)) |

321 |
IndZmat(2,1)=Order(IndZmat0(2,1)) |

322 |
IndZmat(2,2)=Order(IndZmat0(2,2)) |

323 |
IndZmat(3,1)=Order(IndZmat0(3,1)) |

324 |
IndZmat(3,2)=Order(IndZmat0(3,2)) |

325 |
IndZmat(3,3)=Order(IndZmat0(3,3)) |

326 |
DO I=4,Nat |

327 |
DO J=1,4 |

328 |
IndZmat(I,J)=Order(IndZmat0(I,J)) |

329 |
END DO |

330 |
end do |

331 |
IF (DEBUG) THEN |

332 |
WRITE(IOOUT,*) "DBG PathCreate : IndZmat" |

333 |
DO I=1,Nat |

334 |
WRITE(IOOUT,*) (IndZmat(I,J),J=1,5) |

335 |
END DO |

336 |
WRITE(IOOUT,*) "DBG PathCreate : IndZmat0" |

337 |
DO I=1,Nat |

338 |
WRITE(IOOUT,*) (IndZmat0(I,J),J=1,5) |

339 |
END DO |

340 | |

341 |
END IF |

342 | |

343 |
! It is the first call, we have to create all internal coordinates |

344 |
if (debug) WrITE(*,*) "DBG PathCreate L302: First Call, we generate internal coords" |

345 |
DO igeom=1,NgeomI |

346 |
x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |

347 |
y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |

348 |
z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |

349 | |

350 |
if (debug) WRITE(*,*) "DBG PathCreate L308: we generate internal coords for geom",IGeom |

351 | |

352 |
Call ConvXyzZmat(Nat,x,y,z,IndZmat0,val_zmatTmp) |

353 |
IntCoordI(IGeom,1)=val_zmatTmp(2,1) |

354 |
IntCoordI(IGeom,2)=val_zmatTmp(3,1) |

355 |
IntCoordI(IGeom,3)=val_zmatTmp(3,2)*Pi/180. |

356 |
Idx=4 |

357 |
DO i=4,nat |

358 |
IntCoordI(IGeom,Idx)=val_zmatTmp(i,1) |

359 |
Idx=Idx+1 |

360 |
IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180. |

361 |
Idx=Idx+1 |

362 | |

363 |
! Some tests to check that the dihedral angles are similar... this is |

364 |
! important if they are close to Pi. |

365 |
Delta=0. |

366 | |

367 |
if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN |

368 |
if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN |

369 |
Delta=-1.0d0*TwoPi |

370 |
ELSE |

371 |
Delta=TwoPi |

372 |
END IF |

373 |
if (debug) THEN |

374 |
WRITE(*,*) Nom(atome(i)), & |

375 |
abs(val_zmatTMP(i,3)-val_zmat(i,3)), & |

376 |
val_zmatTMP(i,3),val_zmat(i,3),Delta, & |

377 |
val_zmatTMP(i,3)-Delta |

378 |
END IF |

379 |
END IF |

380 |
IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi |

381 |
Idx=Idx+1 |

382 |
END DO |

383 |
END DO ! End Loop on Igeom |

384 | |

385 |
CASE ('MIXED') |

386 |
! IN case we are using mixed coordinate system, we have to |

387 |
! we generate the internal coordinates with a Z-Matrix |

388 |
Renum=.TRUE. |

389 |
! |

390 |
call Calc_mixed_frag(Nat,atome,IndZmat0,val_zmat,x,y,z, & |

391 |
r_cov,fact,frozen,cart,ncart) |

392 | |

393 |
if (debug) WRITE(*,*) 'Back from Calc_Mixed' |

394 | |

395 |
! How to test this partial zmat ? it must be good as it was automatically generated... |

396 |
! IF (.NOT.test_zmat(nat,IndZmat0)) THEN |

397 |
! WRITE(IOOUT,*) "I cannot generate a valid Zmat... help" |

398 |
! STOP |

399 |
! END IF |

400 | |

401 |
! The IndZmat0 generated by Calc_zmat refers to the order of the input x,y,z |

402 |
! However, we need it with the new order so that we can use the zmat |

403 |
! coordinates to generate cartesian coordinates. |

404 |
! So, we have to reorder it... |

405 |
DO I=1,Nat |

406 |
Order(IndZmat0(I,1))=I |

407 |
OrderInv(I)=IndZmat0(I,1) |

408 |
! I let X0 being the Cartesian coordinates of the reference geometry. Before it was |

409 |
! those of the First geometry (which was also the reference one). This might change! |

410 |
X0(i)=X(IndZmat0(i,1)) |

411 |
Y0(i)=Y(IndZmat0(i,1)) |

412 |
Z0(i)=Z(IndZmat0(i,1)) |

413 |
END DO |

414 |
IndZmat=0 |

415 |
IndZmat(1:NCart,:)=-1 |

416 | |

417 |
DO I=1,Nat |

418 |
IndZmat(I,1)=Order(IndZmat0(I,1)) |

419 |
END DO |

420 | |

421 |
SELECT CASE (NCart) |

422 |
CASE (1) |

423 |
IndZmat(2,2)=Order(IndZmat0(2,2)) |

424 |
IndZmat(3,2)=Order(IndZmat0(3,2)) |

425 |
IndZmat(3,3)=Order(IndZmat0(3,3)) |

426 |
IdX=4 |

427 |
CASE (2) |

428 |
IndZmat(3,2)=Order(IndZmat0(3,2)) |

429 |
IndZmat(3,3)=Order(IndZmat0(3,3)) |

430 |
Idx=4 |

431 |
CASE DEFAULT |

432 |
Idx=NCart+1 |

433 |
END SELECT |

434 | |

435 |
DO I=Idx,Nat |

436 |
DO J=1,4 |

437 |
IndZmat(I,J)=Order(IndZmat0(I,J)) |

438 |
END DO |

439 |
end do |

440 |
IF (DEBUG) THEN |

441 |
WRITE(IOOUT,*) "DBG PathCreate : IndZmat - Mixed" |

442 |
DO I=1,Nat |

443 |
WRITE(IOOUT,*) (IndZmat(I,J),J=1,5) |

444 |
END DO |

445 |
END IF |

446 | |

447 |
! It is the first call, we have to create all internal coordinates |

448 |
! Idx=1 |

449 |
! DO I=1,NCart |

450 |
! IntCoordI(1,Idx)=val_zmat(I,1) |

451 |
! IntCoordI(1,Idx+1)=val_zmat(I,2) |

452 |
! IntCoordI(1,Idx+2)=val_zmat(I,3) |

453 |
! Idx=Idx+3 |

454 |
! END DO |

455 |
! SELECT CASE (NCART) |

456 |
! CASE (1) |

457 |
! IntCoordI(1,Idx)=val_zmat(2,1) |

458 |
! IntCoordI(1,Idx+1)=val_zmat(3,1) |

459 |
! IntCoordI(1,Idx+2)=val_zmat(3,2)*Pi/180. |

460 |
! Idx=Idx+3 |

461 |
! IBeg=4 |

462 |
! CASE (2) |

463 |
! IntCoordI(1,Idx)=val_zmat(3,1) |

464 |
! IntCoordI(1,Idx+1)=val_zmat(3,2)*Pi/180. |

465 |
! Idx=Idx+2 |

466 |
! IBeg=4 |

467 |
! CASE DEFAULT |

468 |
! IBeg=NCart+1 |

469 |
! END SELECT |

470 |
! DO i=iBeg,Nat |

471 |
! IntCoordI(1,Idx)=val_zmat(i,1) |

472 |
! Idx=Idx+1 |

473 |
! DO j=2,3 |

474 |
! IntCoordI(1,Idx)=val_zmat(i,j)*Pi/180. |

475 |
! Idx=Idx+1 |

476 |
! END DO |

477 |
! END DO |

478 | |

479 |
! If (debug) WRITE(*,*) "Idx, NCoord",Idx-1,NCoord |

480 | |

481 |
DO Igeom=1,NgeomI |

482 |
x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |

483 |
y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |

484 |
z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |

485 | |

486 |
Call ConvXyzMixed(Nat,Ncart,x,y,z,IndZmat0,val_zmatTmp(1,1)) |

487 |
Idx=1 |

488 |
DO I=1,NCart |

489 |
IntCoordI(IGeom,Idx)=val_zmatTmp(I,1) |

490 |
IntCoordI(IGeom,Idx+1)=val_zmatTmp(I,2) |

491 |
IntCoordI(IGeom,Idx+2)=val_zmatTmp(I,3) |

492 |
Idx=Idx+3 |

493 |
END DO |

494 | |

495 |
SELECT CASE (NCART) |

496 |
CASE (1) |

497 |
IntCoordI(IGeom,Idx)=val_zmatTmp(2,1) |

498 |
Idx=Idx+1 |

499 |
IntCoordI(IGeom,Idx)=val_zmatTmp(3,1) |

500 |
IntCoordI(IGeom,Idx+1)=val_zmatTmp(3,2)*Pi/180. |

501 |
Idx=Idx+2 |

502 |
IBeg=4 |

503 |
CASE (2) |

504 |
IntCoordI(IGeom,Idx)=val_zmatTmp(3,1) |

505 |
IntCoordI(IGeom,Idx+1)=val_zmatTmp(3,2)*Pi/180. |

506 |
Idx=Idx+2 |

507 |
IBeg=4 |

508 |
CASE DEFAULT |

509 |
Ibeg=Ncart+1 |

510 |
END SELECT |

511 |
DO i=IBeg,Nat |

512 |
IntCoordI(IGeom,Idx)=val_zmatTmp(i,1) |

513 |
Idx=Idx+1 |

514 |
IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180. |

515 |
Idx=Idx+1 |

516 | |

517 |
! Some tests to check that the dihedral angles are similar... this is |

518 |
! important if they are close to Pi. |

519 |
Delta=0. |

520 | |

521 |
if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN |

522 |
if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN |

523 |
Delta=-1.0d0*TwoPi |

524 |
ELSE |

525 |
Delta=TwoPi |

526 |
END IF |

527 |
if (debug) THEN |

528 |
WRITE(*,*) Nom(atome(i)), & |

529 |
abs(val_zmatTMP(i,3)-val_zmat(i,3)), & |

530 |
val_zmatTMP(i,3),val_zmat(i,3),Delta, & |

531 |
val_zmatTMP(i,3)-Delta |

532 |
END IF |

533 |
END IF |

534 |
IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi |

535 |
Idx=Idx+1 |

536 |
END DO |

537 | |

538 |
END DO |

539 |
CASE ('CART') |

540 |
DO I=1,Nat |

541 |
Order(I)=I |

542 |
OrderInv(I)=I |

543 |
END DO |

544 |
X0=X |

545 |
Y0=y |

546 |
Z0=z |

547 |
END SELECT |

548 | |

549 |
ELSE ! First iteration of the optimization. Match at L616 |

550 | |

551 |
IGeom=1 |

552 |
x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |

553 |
y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |

554 |
z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |

555 | |

556 |
SELECT CASE (COORD) |

557 |
CASE ('HYBRID') |

558 |
IndZmat0=IndZmat |

559 |
Call ConvXyzZmat(Nat,x,y,z,IndZmat,val_zmat(1,1)) |

560 |
CASE ('ZMAT') |

561 |
IndZmat0=IndZmat |

562 |
val_zmat=0.d0 |

563 |
val_zmat(2,1)=IntCoordI(1,1) |

564 |
val_zmat(3,1)=IntCoordI(1,2) |

565 |
val_zmat(3,2)=IntCoordI(1,3)*180./Pi |

566 | |

567 |
Idx=4 |

568 |
DO iat=4,Nat |

569 |
val_zmat(iat,1)=IntCoordI(1,idx) |

570 |
Idx=Idx+1 |

571 |
do j=2,3 |

572 |
val_zmat(iat,J)=IntCoordI(1,idx)*180./Pi |

573 |
Idx=Idx+1 |

574 |
END DO |

575 |
END DO |

576 | |

577 |
CASE ('MIXED') |

578 |
IndZmat0=IndZmat |

579 |
val_zmat=0.d0 |

580 |
Idx=1 |

581 |
DO I=1,NCart |

582 |
DO J=1,3 |

583 |
val_zmat(i,j)= IntCoordI(1,Idx) |

584 |
Idx=Idx+1 |

585 |
END DO |

586 |
END DO |

587 |
SELECT CASE (NCart) |

588 |
CASE(1) |

589 |
val_zmat(2,1)=IntCoordI(1,Idx) |

590 |
val_zmat(3,1)=IntCoordI(1,Idx+1) |

591 |
val_zmat(3,2)=IntCoordI(1,Idx+2)*180./Pi |

592 |
Idx=Idx+3 |

593 |
IBeg=4 |

594 |
CASE(2) |

595 |
val_zmat(3,1)=IntCoordI(1,Idx) |

596 |
val_zmat(3,2)=IntCoordI(1,Idx)*180./Pi |

597 |
Idx=Idx+2 |

598 |
IBeg=4 |

599 |
CASE DEFAULT |

600 |
IBeg=NCart+1 |

601 |
END SELECT |

602 | |

603 |
DO Iat=IBeg,Nat |

604 |
val_zmat(iat,1)=IntCoordI(1,idx) |

605 |
Idx=Idx+1 |

606 |
do j=2,3 |

607 |
val_zmat(iat,J)=IntCoordI(1,idx)*180./Pi |

608 |
Idx=Idx+1 |

609 |
END DO |

610 |
END DO |

611 |
CASE ('BAKER') |

612 |
! Nothing to be done here. |

613 |
CASE ('CART') |

614 |
! Nothing to be done here. |

615 |
CASE DEFAULT |

616 |
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate.L554.STOP" |

617 |
STOP |

618 |
END SELECT |

619 |
X0=x |

620 |
y0=y |

621 |
z0=z |

622 |
END IF ! First |

623 | |

624 |
if (Print) THEN |

625 |
WRITE(*,*) "PathCreate - L631 - geometries " |

626 |
DO I=1,NGeomI |

627 |
DO J=1,Nat |

628 |
If (renum) THEN |

629 |
Iat=Order(J) |

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

631 |
ELSE |

632 |
Iat=OrderInv(J) |

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

634 |
END IF |

635 |
END DO |

636 |
END DO |

637 |
END IF |

638 | |

639 | |

640 |
! Now that we have a zmat, we will generate all the IntCoodI corresponding... |

641 |
! First one |

642 |
IF (COORD.EQ.'HYBRID') THEN ! Matches at L680 |

643 |
DO IGeom=1,NGeomI |

644 |
x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |

645 |
y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |

646 |
z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |

647 | |

648 |
Call ConvXyzZmat(Nat,x,y,z,IndZmat0,val_zmatTmp(1,1)) |

649 |
IntCoordI(IGeom,1)=val_zmatTmp(2,1) |

650 |
IntCoordI(IGeom,2)=val_zmatTmp(3,1) |

651 |
IntCoordI(IGeom,3)=val_zmatTmp(3,2)*Pi/180. |

652 |
Idx=4 |

653 |
DO i=4,nat |

654 |
IntCoordI(IGeom,Idx)=val_zmatTmp(i,1) |

655 |
Idx=Idx+1 |

656 |
IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180. |

657 |
Idx=Idx+1 |

658 | |

659 |
! Some tests to check that the dihedral angles are similar... this is |

660 |
! important if they are close to Pi. |

661 |
Delta=0. |

662 | |

663 |
if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN |

664 |
if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN |

665 |
Delta=-1.0d0*TwoPi |

666 |
ELSE |

667 |
Delta=TwoPi |

668 |
END IF |

669 |
if (debug) THEN |

670 |
WRITE(*,*) Nom(atome(i)), & |

671 |
abs(val_zmatTMP(i,3)-val_zmat(i,3)), & |

672 |
val_zmatTMP(i,3),val_zmat(i,3),Delta, & |

673 |
val_zmatTMP(i,3)-Delta |

674 |
END IF |

675 |
END IF |

676 |
IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi |

677 |
Idx=Idx+1 |

678 |
END DO |

679 |
END DO |

680 |
END IF ! Matches Coord=Hybrid |

681 | |

682 |
! PFL 24 Nov 2008 -> |

683 |
! This should not be necessary... |

684 |
! ! if we have frozen atoms, we make sure that they do not move between geometries |

685 |
! IF (IntFroz.NE.0) THEN |

686 |
! if (debug) WRITE(*,*) 'DBG PathCreate L641: Number of frozen coordinates=', IntFroz |

687 |
! DO IGeom=2,NGeomI ! For IOpt .GT. 0, NGeomI is equal to NGeomF |

688 |
! IntCoordI(IGeom,1:IntFroz)=IntCoordI(1,1:IntFroz) |

689 |
! END DO |

690 |
! END IF |

691 |
! -> PFL 24 Nov 2008 |

692 | |

693 |
Idx=min(9,NCoord) |

694 | |

695 |
IF (DEBUG.AND.(COORD.NE.'CART')) THEN |

696 |
WRITE(*,*) "DBG PathCreate. L685 IntCoordI(IGeom,:)" |

697 |
DO I=1,NGeomI |

698 |
WRITE(*,'("Geom:",I5,9(1X,F10.4))') I,IntCoordI(I,1:Idx) |

699 |
END DO |

700 |
! We write it also in terms of Zmat |

701 |
IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'HYBRID')) THEN |

702 |
DO I=1,NGeomI |

703 |
WRITE(*,*) 'Zmat for geom ',I,'/',NGeomI,'Xgeom=',Xgeom(I) |

704 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(1,1))) |

705 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(2,1))), & |

706 |
IndZmat(2,2),IntCoordI(I,1) |

707 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(3,1))), & |

708 |
IndZmat(3,2),IntCoordI(I,2),IndZmat(3,3),IntCoordI(I,3)*180./Pi |

709 |
Idx=4 |

710 |
DO J=4,Nat |

711 |
WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(J,1))), & |

712 |
IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), & |

713 |
IntCoordI(I,Idx+1)*180./Pi, & |

714 |
IndZmat(J,4),IntCoordI(I,Idx+2)*180./Pi |

715 |
Idx=Idx+3 |

716 |
END DO |

717 | |

718 |
XyzTmp2=0. |

719 |
! We convert it into Cartesian coordinates |

720 |
XyzTmp2(2,1)=IntCoordI(I,1) |

721 |
d=IntCoordI(I,2) |

722 |
a_val=IntCoordI(I,3) |

723 |
if (Nat.GE.3) THEN |

724 |
if (IndZmat(3,2).EQ.1) THEN |

725 |
XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val) |

726 |
ELSE |

727 |
XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val) |

728 |
ENDIF |

729 |
XyzTmp2(3,2)=d*sin(a_val) |

730 |
ENDIF |

731 | |

732 |
Idx=4 |

733 |
DO iat=4,Nat |

734 |
val_zmatTMP(iat,1)=IntCoordI(I,idx) |

735 |
Idx=Idx+1 |

736 |
DO j=2,3 |

737 |
val_zmatTMP(iat,J)=IntCoordI(I,idx)*180./Pi |

738 |
Idx=Idx+1 |

739 |
END DO |

740 |
END DO |

741 |
DO iat=4,Nat |

742 |
call ConvertZmat_cart(iat,IndZmat,val_zmatTMP, & |

743 |
XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3)) |

744 |
END DO |

745 |
DO Iat=1,nat |

746 |
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(Iat,1))),& |

747 |
XyzTmp2(iat,1:3) |

748 |
END DO |

749 | |

750 |
END DO |

751 |
END IF |

752 |
IF (COORD.EQ.'MIXED') THEN |

753 |
DO I=1,NGeomI |

754 |
WRITE(*,*) 'Cart+Zmat for geom ',I,'/',NGeomI,'Xgeom=',Xgeom(I) |

755 |
Idx=1 |

756 |
XyzTmp2=0. |

757 |
DO J=1,NCART |

758 |
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(1,1))), & |

759 |
IntCoordI(I,Idx:Idx+2) |

760 |
XyzTmp2(J,1:3)=IntCoordI(I,Idx:Idx+2) |

761 |
Idx=Idx+3 |

762 |
END DO |

763 |
SELECT CASE (NCart) |

764 |
CASE (1) |

765 |
J=2 |

766 |
WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), & |

767 |
IndZmat(J,2),IntCoordI(I,Idx) |

768 |
Idx=Idx+1 |

769 |
J=3 |

770 |
WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), & |

771 |
IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), & |

772 |
IntCoordI(I,Idx+1)*180./Pi |

773 |
Idx=Idx+2 |

774 |
Ibeg=4 |

775 |
CASE (2) |

776 |
J=3 |

777 |
WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), & |

778 |
IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), & |

779 |
IntCoordI(I,Idx+1)*180./Pi |

780 |
Idx=Idx+2 |

781 |
Ibeg=4 |

782 |
CASE DEFAULT |

783 |
IBeg=NCart+1 |

784 |
END SELECT |

785 |
DO J=IBeg,Nat |

786 |
IF (J.EQ.2) & |

787 |
WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), & |

788 |
IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), & |

789 |
IntCoordI(I,Idx+1)*180./Pi, & |

790 |
IndZmat(J,4),IntCoordI(I,Idx+2)*180./Pi |

791 |
Idx=Idx+3 |

792 |
END DO |

793 | |

794 |
! We convert it into Cartesian coordinates |

795 |
IntCoordTmp=IntCoordI(I,1:NCoord) |

796 |
Call Mixed2cart(Nat,IndZmat,IntCoordTmp,XyzTmp2) |

797 | |

798 |
DO Iat=1,nat |

799 |
WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(Iat,1))),& |

800 |
XyzTmp2(iat,1:3) |

801 |
END DO |

802 | |

803 |
END DO |

804 |
END IF |

805 |
END IF ! matches IF (DEBUG.AND.(COORD.NE.'CART')) THEN |

806 | |

807 |
if (debug) WRITE(*,*) "Before interpolation, PathCreate.f90, L740, IOpt=",IOpt, & |

808 |
"ISpline=", ISpline |

809 | |

810 |
if (Print) THEN |

811 |
WRITE(*,*) "PathCreate - L811 - geometries " |

812 |
DO I=1,NGeomI |

813 |
DO J=1,Nat |

814 |
If (renum) THEN |

815 |
Iat=Order(J) |

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

817 |
ELSE |

818 |
Iat=OrderInv(J) |

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

820 |
END IF |

821 |
END DO |

822 |
END DO |

823 |
END IF |

824 | |

825 | |

826 |
! Now comes the Interpolation: |

827 |
IF ((NGeomI>2).AND.(IOpt.GE.ISpline)) THEN |

828 |
Linear=.FALSE. |

829 |
! We have at least three geometries so we can apply cubic spline |

830 |
SELECT CASE (COORD) |

831 |
CASE ('ZMAT','HYBRID','MIXED') |

832 |
DO I=1,NCoord |

833 |
! We compute the spline coefficients |

834 |
call spline(xGeom,IntCoordI(1,I),NGeomI,Inf,Inf, Coef(1,I)) |

835 |
if (debug) THEN |

836 |
WRITE(*,*) 'Coef Spline for coord',I |

837 |
WRITE(*,*) Coef(1:NGeomI,I) |

838 |
END IF |

839 |
END DO |

840 |
CASE ('BAKER') |

841 |
! We compute the spline coefficients |

842 |
! xGeom(NGeomI), IntCoordI(NGeomI,3*Nat-6) declared in Path_module.f90 |

843 |
! can also be used for the Baker's internal coordinates. |

844 |
! we need to get Baker's internal coordinates from cartesian coordinates. |

845 |
! as we have subroutine ConvXyzmat(...) in Z matrix case. |

846 |
! Coef(NGeomI,NCoord) or (NGeomI,N3at) has INTENT(OUT) |

847 |
DO I=1,NCoord ! NCoord=3*Nat-6-symmetry_eliminated. |

848 |
call spline(xGeom,IntCoordI(1,I),NGeomI,Inf,Inf,Coef(1,I)) |

849 |
! IF (debug) THEN |

850 |
! WRITE(*,*) 'Coef Spline for coord',I |

851 |
! WRITE(*,*) Coef(1:NGeomI,I) |

852 |
! END IF |

853 |
END DO |

854 |
CASE ('CART') |

855 | |

856 |
!! PFL 2011 June 22 -> |

857 |
!! The alignment procedure from Extrapol_cart has been moved |

858 |
!! here, before computation of the spline coefficients |

859 | |

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

861 | |

862 |
if (debug) THEN |

863 |
WRITE(*,*) "DBG PathCreate - Calc Spline coeff CART- Initial geometries" |

864 |
DO IGeom=1,NGeomI |

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

866 |
DO I=1,Nat |

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

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

869 |
END DO |

870 |
END DO |

871 |
END IF |

872 | |

873 | |

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

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

876 | |

877 |
DO IGeom=1,NGeomI |

878 | |

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

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

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

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

883 |
! PFL 24 Nov 2008 -> |

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

885 |
! PFL 8 Feb 2011 -> |

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

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

888 |
IF (Align) THEN |

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

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

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

892 |
FrozAtoms,MRot) |

893 |
ELSE |

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

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

896 |
MRot,rmsd,.TRUE.,.TRUE.) |

897 |
END IF |

898 |
! <- PFL 24 Nov 2008 |

899 |
END IF |

900 |
! -> PFL 8 Feb 2011 |

901 |
! END IF |

902 | |

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

904 |
END DO |

905 | |

906 |
if (debug) THEN |

907 |
WRITE(*,*) "DBG PathCreate - Calc Spline coeff CART- Aligned geometries" |

908 |
DO J=1, NGeomI |

909 |
WRITE(IOOUT,*) Nat |

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

911 |
DO i=1,Nat |

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

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

914 |
END DO |

915 |
END DO |

916 |
END IF |

917 | |

918 |
!! <- PFL 2011 June 22 |

919 | |

920 |
DO I=1,Nat |

921 |
! We compute the spline coefficients |

922 |
call spline(xGeom,XyzGeomI(1,1,I),NGeomI,Inf,Inf, Coef(1,3*I-2)) |

923 |
call spline(xGeom,XyzGeomI(1,2,I),NGeomI,Inf,Inf, Coef(1,3*I-1)) |

924 |
call spline(xGeom,XyzGeomI(1,3,I),NGeomI,Inf,Inf, Coef(1,3*I)) |

925 |
if (debug) THEN |

926 |
WRITE(*,*) 'Coef Spline for coord',i |

927 |
WRITE(*,*) Coef(1:NGeomI,I) |

928 |
END IF |

929 |
END DO |

930 |
CASE DEFAULT |

931 |
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate.L783.STOP" |

932 |
STOP |

933 |
END SELECT |

934 |
ELSE ! matches IF ((NGeomI>2).AND.(Iopt.Ge.ISpline)) THEN |

935 |
! On n'a que deux images. On fera une interpolation lineaire, |

936 |
! qui est un cas particulier avec Coef=0 |

937 |
if (debug.AND.(NGeomI.EQ.2)) WRITE(*,*) "DBG : Only 2 starting geometries" |

938 |
if (debug.AND.(Iopt.LE.ISpline)) WRITE(*,*) "DBG : Not enough path -> linear int" |

939 |
Coef=0.d0 |

940 |
Linear=.TRUE. |

941 |
END IF ! matches IF ((NGeomI>2).AND.(Iopt.Ge.ISpline)) THEN |

942 |
! Shall we redistribute the points along the path ? |

943 |
if (debug) WRITE(*,*) "PathCreate.f90, L795, IOpt=", IOpt, "IReparam=", IReparam |

944 | |

945 |
IF (MOD(IOpt,IReparam).EQ.0) THEN |

946 |
! We generate the Path, in two steps: |

947 |
! i) we calculate the length of the mass weighted (MW) path |

948 |
! ii) we sample the path to get the geometries |

949 | |

950 |
NGeomS=NGeomI |

951 |
SELECT CASE (COORD) |

952 |
CASE ('ZMAT','HYBRID') |

953 |
! We call the extrapolation once to get the path length |

954 |
dist=Inf |

955 |
Call ExtraPol_int(IOpt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF) |

956 |
! we have the length, we calculate the distance between two points |

957 | |

958 |
if (s.LE.1e-6) THEN |

959 |
IF (OptGeom.LE.0) THEN |

960 |
WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization." |

961 |
WRITE(*,*) "ERROR !!! Stop !" |

962 |
STOP |

963 |
END IF |

964 |
ELSE |

965 |
! we have the length, we calculate the distance between two points |

966 |
dist=s/Real(NGeomf-1) |

967 |
if (debug) WRITE(*,*) 'Dbg PathCreate ZMAT s, dist:',s,dist |

968 |
! we call it again to obtain the geometries ! |

969 |
Call ExtraPol_int(iopt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF) |

970 |
END IF |

971 |
CASE ('BAKER') |

972 |
! We generate the Path, in two steps: |

973 |
! i) we calculate the length of the mass weighted (MW) path |

974 |
! ii) we sample the path to get the geometries |

975 | |

976 |
! We call the extrapolation first time to get the path length |

977 |
! iopt: Number of the cycles for the optimization. |

978 |
! s: a real number and has INTENT(OUT), probably returns the |

979 |
! length of the path. |

980 |
! dist=Inf=> ExtraPol_baker has to calculate the length of the path. |

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

982 |
! Xgeom: Xgeom(NGeomI) has INTENT(IN) |

983 |
! Coef(NGeomI,NCoord) is obtained from spline interpolation subroutine. |

984 |
! XgeomF: XGeomF(NGeomF) has INTENT(OUT) |

985 |
! Which reference coordinate geometry X0, Y0, Z0 to use??? |

986 |
dist=Inf |

987 |
!if (debug) WRITE(*,*) "Before the call of ExtraPol_baker in redistribution & |

988 |
!of the points, in PathCreate.f90, L836" |

989 |
!Print *, 'Before ExtraPol_baker, PathCreate.90, L848, IntCoordI=' |

990 |
!Print *, IntCoordI |

991 |
Call ExtraPol_baker(s,dist,x0,y0,z0,xgeom,Coef,XgeomF) |

992 |
!Print *, 'After ExtraPol_baker, PathCreate.90, L848, IntCoordI=' |

993 |
!Print *, IntCoordI |

994 |
!if (debug) WRITE(*,*) "After the call of ExtraPol_baker in the redistribution & |

995 |
! of the points, in PathCreate.f90, L843" |

996 | |

997 |
! we have the lenght, we calculate the distance between two points |

998 |
!Print *, 'After first call of ExtraPol_baker in PathCreate.f90, s=', s |

999 | |

1000 |
if (s.LE.1e-6) THEN |

1001 |
IF (OptGeom.LE.0) THEN |

1002 |
WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization." |

1003 |
WRITE(*,*) "ERROR !!! Stop !" |

1004 |
STOP |

1005 |
END IF |

1006 |
ELSE |

1007 | |

1008 |
dist=s/Real(NGeomf-1) |

1009 |
!if (debug) WRITE(*,*) 'Dbg PathCreate ZMAT s, dist:',s,dist |

1010 | |

1011 |
! we call it again to obtain the geometries ! |

1012 |
Call ExtraPol_baker(s,dist,x0,y0,z0,xgeom,Coef,XgeomF) |

1013 |
WRITE(*,*) "At the end of Baker case in the redistribution of the points, in PathCreate.f90, L966" |

1014 |
END IF |

1015 | |

1016 |
DO IGeom=1,NGeomF |

1017 |
WRITE(*,*) "UMatF for IGeom=",IGeom |

1018 |
DO I=1,3*Nat-6 |

1019 |
WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I) |

1020 |
END DO |

1021 |
END DO |

1022 | |

1023 |
CASE ('MIXED') |

1024 |
! WRITE(*,*) "IOIOIOIOIOIOIOI" |

1025 |
! We call the extrapolation once to get the path length |

1026 |
Dist=Inf |

1027 |
Call ExtraPol_Mixed(s,dist,x0,y0,z0,xgeom,Coef) |

1028 | |

1029 |
if (s.LE.1e-6) THEN |

1030 |
IF (OptGeom.LE.0) THEN |

1031 |
WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization." |

1032 |
WRITE(*,*) "ERROR !!! Stop !" |

1033 |
STOP |

1034 |
END IF |

1035 |
ELSE |

1036 | |

1037 |
! we have the length, we calculate the distance between two points |

1038 |
dist=s/Real(NGeomf-1) |

1039 |
if (debug) WRITE(*,*) 'Dbg PathCreate Mixed s, dist:',s,dist |

1040 | |

1041 |
! we call it again to obtain the geometries ! |

1042 |
Call ExtraPol_Mixed(s,dist,x0,y0,z0,xgeom,Coef) |

1043 | |

1044 |
IF (debug) THEN |

1045 |
DO I=1,NGeomF |

1046 |
WRITE(*,'(1X,A,12(1X,F10.5))') "IntTangent:",IntTangent(I,1:NCoord) |

1047 |
END DO |

1048 |
END IF |

1049 |
END IF |

1050 |
CASE ('CART') |

1051 |
! We call the extrapolation once to get the path length |

1052 |
Dist=Inf |

1053 |
Call ExtraPol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef) |

1054 | |

1055 |
if (s.LE.1e-6) THEN |

1056 |
IF (OptGeom.LE.0) THEN |

1057 |
WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization." |

1058 |
WRITE(*,*) "ERROR !!! Stop !" |

1059 |
STOP |

1060 |
END IF |

1061 |
ELSE |

1062 |
! we have the lenght, we calculate the distance between two points |

1063 |
dist=s/Real(NGeomf-1) |

1064 |
if (debug) WRITE(*,*) 'Debug s, dist:',s,dist |

1065 | |

1066 |
! we call it again to obtain the geometries ! |

1067 |
Call ExtraPol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef) |

1068 |
END IF |

1069 |
END SELECT |

1070 |
ELSE |

1071 |
! PFL 3 Jan 2008 -> |

1072 |
! Tangents are not recalculated if the points are not reparameterized |

1073 |
! PFL 6 Apr 2008 -> |

1074 |
! Unless user has asked to reparameterized them ! |

1075 |
IF (mod(Iopt,iReparamT).EQ.0) THEN |

1076 |
XGeomF=XGeom |

1077 |
NGeomS=NGeomI |

1078 |
Call CalcTangent(XgeomF,NGeomS,xgeom,Coef) |

1079 |
! |

1080 |
END IF |

1081 |
! <- PFL 3 Jan 2008 / 6 Apr 2008 |

1082 |
END IF ! if (mod(Iop,iReparam).EQ.0) |

1083 |
if (debug) WRITE(*,*) "PathCreate.f90, L885" |

1084 | |

1085 |
! PFL 22.Nov.2007 -> |

1086 |
! We do not write tangents for baker coordinates because they are hard to |

1087 |
! interpret for a chemist. |

1088 | |

1089 |
IF (debug.AND.(COORD.NE."BAKER")) THEN |

1090 |
! we write the tangents :-) |

1091 |
DO I=1,NGeomF |

1092 |
SELECT CASE(Coord) |

1093 |
CASE ('ZMAT','MIXED') |

1094 |
IntCoordTmp=IntTangent(I,:) |

1095 |
WRITE(*,'(1X,A,12(1X,F10.5))') "IntCoordTmp:",IntCoordTmp(1:NCoord) |

1096 |
WRITE(*,'(1X,A,12(1X,F10.5))') "IntTangent:",IntTangent(I,1:NCoord) |

1097 |
CASE ('CART','HYBRID') |

1098 |
IntCoordTmp=XyzTangent(I,:) |

1099 |
CASE DEFAULT |

1100 |
WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate. STOP" |

1101 |
STOP |

1102 |
END SELECT |

1103 |
WRITE(Title,"('Tangents for geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF |

1104 |
Call PrintGeom(Title,Nat,Ncoord,IntCoordTmp,Coord,6,Atome,Order,OrderInv,IndZmat) |

1105 |
END DO |

1106 |
END IF |

1107 | |

1108 | |

1109 |
IF (DEBUG.AND.(COORD.NE.'CART')) THEN |

1110 |
Idx=min(12,NCoord) |

1111 |
WRITE(*,*) "DBG PathCreate. IntCoordF(IGeom,:)" |

1112 |
DO I=1,NGeomF |

1113 |
WRITE(*,'("Geom:",I5,12(1X,F10.4))') I,IntCoordF(I,1:Idx) |

1114 |
END DO |

1115 |
END IF |

1116 | |

1117 |
DEALLOCATE(Coef,val_zmat,val_zmatTMP) |

1118 |
DEALLOCATE(X0,Y0,Z0,Indzmat0) |

1119 |
DEALLOCATE(X,Y,Z,Xgeom,XgeomF) |

1120 |
DEALLOCATE(XyzTmp,XyzTmp2) |

1121 | |

1122 |
if (debug) Call header("PathCreate over") |

1123 | |

1124 |
END SUBROUTINE PathCreate |

1125 |