## root / src / ConvertBakerInternal_cart.f90 @ 8

Historique | Voir | Annoter | Télécharger (23,95 ko)

1 |
SUBROUTINE ConvertBakerInternal_cart(IntCoord_k,IntCoord,x_k,y_k,z_k,x,y,z,XPrim,XPrimRef) |
---|---|

2 | |

3 |
!================================================================ |

4 |
! Converts the positions in Baker Coordinates to cartesian coordinates |

5 |
!================================================================ |

6 |
! |

7 |
! Input: |

8 |
! Incoord_k(NCoord) REAL: Starting internal coordinates |

9 |
! IntCoord(NCoord) REAL: Coordinate to convert to cartesian coordinates |

10 |
! x_k,y_k,z_k(Nat) REAL: Starting cartesian geometry |

11 |
! XPrimRef(NPrim) REAL: Reference values for Primitives coordinates (mainly for dihedral) |

12 |
! |

13 |
! Output: |

14 |
! x,y,z(Nat) REAL: Cartesian coordinates corresponding to IntCoord |

15 |
! XPrim(NPrim) REAL: Primitives for the cartesian coordinates. (This is a by-product) |

16 |
! |

17 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |

18 | |

19 |
! BTransInv,UMat should be passed as arguments of this subroutine. |

20 | |

21 | |

22 |
! IdxGeom: geometry index. |

23 |
! Internal coordinates IntCoord(NCoord) is to be converted |

24 |
! into Cartesian coordinates x(Nat),y(Nat),z(Nat). |

25 | |

26 |
use Path_module, only : Nat,AtName,BMat_BakerT,NPrim,& |

27 |
BBT,BBT_inv,BprimT,ScanCoord,Coordinate,& |

28 |
BTransInv_local,Xprimitive_t,& |

29 |
NCoord,UMat_local |

30 |
! UMat_local is allocated in Calc_baker.f90 |

31 |
! XyzGeomI(NGeomI,3,Nat),BMat_BakerT(3*Nat,NCoord),a0=0.529177249d0 |

32 | |

33 |
Use Io_Module |

34 |
use VarTypes |

35 |
IMPLICIT NONE |

36 | |

37 |
REAL(KREAL), INTENT(OUT) :: x(Nat),y(Nat),z(Nat) |

38 |
REAL(KREAL), INTENT(IN) :: x_k(Nat),y_k(Nat),z_k(Nat) |

39 |
REAL(KREAL), INTENT(IN) :: IntCoord(NCoord) |

40 |
REAL(KREAL), INTENT(INOUT) :: IntCoord_k(NCoord) |

41 |
REAL(KREAL), INTENT(IN) :: XPrimRef(NPrim) |

42 |
REAL(KREAL), INTENT(OUT) :: XPrim(NPrim) |

43 | |

44 | |

45 |
! IdxGeom: geometry index. |

46 |
Integer(KINT) :: I, J, Counter |

47 |
REAL(KREAL) :: normDeltaX_int |

48 | |

49 |
REAL(KREAL), ALLOCATABLE :: DeltaX_int(:) !DeltaX_int(NPrim,3*Nat-6), 3*Nat-6?? |

50 |
REAL(KREAL), ALLOCATABLE :: DeltaX_cart(:) !DeltaX_cart(3*Nat) |

51 |
REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat) |

52 |
REAL(KREAL), ALLOCATABLE :: X_cart_k(:,:) ! (3,Nat) |

53 |
REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim) |

54 |
REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim) |

55 |
REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:) |

56 | |

57 |
REAL(KREAL) :: Angle, angle_d |

58 |
REAL(KREAL), PARAMETER :: Crit=1e-08 |

59 | |

60 |
EXTERNAL Angle, angle_d |

61 |
LOGICAL :: debug, FAIL |

62 | |

63 |
INTEGER(KINT), PARAMETER :: NbItMax=50 |

64 | |

65 | |

66 |
INTERFACE |

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

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

69 |
logical :: isValid |

70 |
END function VALID |

71 | |

72 | |

73 | |

74 |
SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef) |

75 |
! |

76 |
! This subroutine uses the description of a list of Coordinates |

77 |
! to compute the values of the coordinates for a given geometry. |

78 |
! |

79 |
!!!!!!!!!! |

80 |
! Input: |

81 |
! Na: INTEGER, Number of atoms |

82 |
! x,y,z(Na): REAL, cartesian coordinates of the considered geometry |

83 |
! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes |

84 |
! NPrim, INTEGER: Number of coordinates to compute |

85 |
! |

86 |
! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for |

87 |
! a former geometry. Useful for Dihedral angles evolution... |

88 | |

89 |
!!!!!!!!!!! |

90 |
! Output: |

91 |
! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates |

92 |
! |

93 |
!!!!!!!!! |

94 | |

95 |
Use VarTypes |

96 |
Use Io_module |

97 |
Use Path_module, only : pi |

98 | |

99 |
IMPLICIT NONE |

100 | |

101 |
Type (ListCoord), POINTER :: Coordinate |

102 |
INTEGER(KINT), INTENT(IN) :: Nat,NPrim |

103 |
REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat) |

104 |
REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) |

105 |
REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim) |

106 | |

107 |
END SUBROUTINE CALC_XPRIM |

108 |
END INTERFACE |

109 | |

110 | |

111 |
debug=valid('ConvertBakerInternal_cart') |

112 |
if (debug) WRITE(*,*) '=======Entering ConvertBakerInternal_cart=======' |

113 | |

114 |
FAIL = .FALSE. |

115 | |

116 |
!!!!!! |

117 |
! Allocation phase |

118 | |

119 |
ALLOCATE(DeltaX_int(NCoord),DeltaX_cart(3*Nat)) |

120 |
ALLOCATE(UMat_tmp(NPrim,NCoord)) |

121 |
ALLOCATE(X_cart_k(3,Nat)) |

122 |
ALLOCATE(Geom(3,Nat),BprimT(3*Nat,NPrim)) |

123 |
ALLOCATE(BBT(NCoord,NCoord)) |

124 |
ALLOCATE(BBT_inv(NCoord,NCoord)) |

125 |
ALLOCATE(Gmat(NPrim,NPrim)) |

126 |
ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim)) |

127 | |

128 | |

129 |
if (debug) THEN |

130 |
WRITE(*,*) "Checking purposes...." |

131 |
WRITE(*,*) "Trying ot convert Xint initial (IntCoord)" |

132 |
WRITE(*,'(50(1X,F12.8))') IntCoord(:) |

133 |
WRITE(*,*) "BTransInv_local" |

134 |
DO J=1,3*Nat |

135 |
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) |

136 |
END DO |

137 |
WRITE(*,*) "Xint initial (IntCoord_k)" |

138 |
WRITE(*,'(50(1X,F12.8))') IntCoord_k(:) |

139 |
WRITE(*,*) "Cartesian coordinates" |

140 |
WRITE(*,*) Nat |

141 |
WRITE(*,*) 'GeomCart in ConvertBakerInternal_cart' |

142 |
DO I=1,Nat |

143 |
WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,X_k(I),Y_k(i),Z_k(i) |

144 |
END DO |

145 |
WRITE(*,*) "Calculating actual values using x_k,y_k,z_k" |

146 | |

147 |
! WRITE(*,*) "First, with Calc_XPrim" |

148 |
! WRITE(*,*) "Coordinate:",associated(Coordinate) |

149 |
Call Calc_Xprim(nat,x_k,y_k,z_k,Coordinate,NPrim,XPrimitive_t,XPrimRef) |

150 | |

151 |
! I=0 ! index for the NPrim (NPrim is the number of ! primitive coordinates). |

152 |
! ScanCoord=>Coordinate |

153 |
! DO WHILE (Associated(ScanCoord%next)) |

154 |
! I=I+1 |

155 |
! SELECT CASE (ScanCoord%Type) |

156 |
! CASE ('BOND') |

157 |
! Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2, '=', Xprimitive_t(I) |

158 |
! CASE ('ANGLE') |

159 |
! Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2,'-',ScanCoord%At3, '=', Xprimitive_t(I)*180./Pi |

160 |
! CASE ('DIHEDRAL') |

161 |
! Print *, I, ':', ScanCoord%At1, '-',ScanCoord%At2,'-',ScanCoord%At3,'-',& |

162 |
! ScanCoord%At4,'=',Xprimitive_t(I)*180./Pi |

163 |
! END SELECT |

164 |
! ScanCoord => ScanCoord%next |

165 |
! END DO ! matches DO WHILE (Associated(ScanCoord%next)) |

166 | |

167 |
! WRITE(*,*) "Second with normal proc" |

168 |
! I=0 ! index for the NPrim (NPrim is the number of ! primitive coordinates). |

169 |
! Xprimitive_t=0.d0 |

170 |
! ScanCoord=>Coordinate |

171 |
! DO WHILE (Associated(ScanCoord%next)) |

172 |
! I=I+1 |

173 |
! SELECT CASE (ScanCoord%Type) |

174 |
! CASE ('BOND') |

175 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x_k,y_k,z_k,vx2,vy2,vz2,Norm2) |

176 |
! Xprimitive_t(I) = Norm2 |

177 |
! Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2, '=', Xprimitive_t(I) |

178 |
! CASE ('ANGLE') |

179 |
! Call vecteur(ScanCoord%At2,ScanCoord%At3,x_k,y_k,z_k,vx1,vy1,vz1,Norm1) |

180 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x_k,y_k,z_k,vx2,vy2,vz2,Norm2) |

181 |
! Xprimitive_t(I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |

182 |
! Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2,'-',ScanCoord%At3, '=', Xprimitive_t(I)*180./Pi |

183 |
! CASE ('DIHEDRAL') |

184 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x_k,y_k,z_k,vx1,vy1,vz1,Norm1) |

185 |
! Call vecteur(ScanCoord%At2,ScanCoord%At3,x_k,y_k,z_k,vx2,vy2,vz2,Norm2) |

186 |
! Call vecteur(ScanCoord%At3,ScanCoord%At4,x_k,y_k,z_k,vx3,vy3,vz3,Norm3) |

187 |
! Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5) |

188 |
! Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4) |

189 |
! !!! PFL 28 Feb 2008 Test to be sure that angle does not change by more than Pi |

190 |
! ! Xprimitive_t(I)=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)*Pi/180. |

191 |
! DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |

192 |
! vx2,vy2,vz2,norm2) |

193 |
! Xprimitive_t(I) = DiheTmp*Pi/180. |

194 |
! ! We treat large dihedral angles differently as +180=-180 mathematically and physically |

195 |
! ! but this causes lots of troubles when converting baker to cart. |

196 |
! ! So we ensure that large dihedral angles always have the same sign |

197 |
! if (abs(ScanCoord%SignDihedral).NE.1) THEN |

198 |
! ScanCoord%SignDihedral=Int(Sign(1.,DiheTmp)) |

199 |
! ELSE |

200 |
! If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN |

201 |
! Xprimitive_t(I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi |

202 |
! END IF |

203 |
! END IF |

204 | |

205 |
! IF (abs(Xprimitive_t(I)-XPrimRef(I)).GE.Pi) THEN |

206 |
! if ((Xprimitive_t(I)-XPrimRef(I)).GE.Pi) THEN |

207 |
! Xprimitive_t(I)=Xprimitive_t(I)-2.d0*Pi |

208 |
! ELSE |

209 |
! Xprimitive_t(I)=Xprimitive_t(I)+2.d0*Pi |

210 |
! END IF |

211 |
! END IF |

212 | |

213 |
! Print *, I, ':', ScanCoord%At1, '-',ScanCoord%At2,'-',ScanCoord%At3,'-',& |

214 |
! ScanCoord%At4,'=',Xprimitive_t(I)*180./Pi |

215 |
! END SELECT |

216 |
! ScanCoord => ScanCoord%next |

217 |
! END DO ! matches DO WHILE (Associated(ScanCoord%next)) |

218 | |

219 |
IntCoord_k=0.d0 |

220 |
DO I=1, NPrim |

221 |
! Transpose of UMat is needed below, that is why UMat(I,:). |

222 |
IntCoord_k(:)=IntCoord_k(:)+UMat_local(I,:)*Xprimitive_t(I) |

223 |
END DO |

224 |
WRITE(*,*) "Xint (intCoord_k) cooresponding to x_k,y_k,z_k" |

225 |
WRITE(*,'(50(1X,F12.8))') IntCoord_k(:) |

226 |
WRITE(*,*) "UMat_local" |

227 |
DO I=1, NPrim |

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

229 |
END DO |

230 | |

231 |
! Geom(1,:)=X_k(1:Nat)/a0 |

232 |
! Geom(2,:)=y_k(1:Nat)/a0 |

233 |
! Geom(3,:)=z_k(1:Nat)/a0 |

234 | |

235 |
Geom(1,:)=X_k(1:Nat) |

236 |
Geom(2,:)=y_k(1:Nat) |

237 |
Geom(3,:)=z_k(1:Nat) |

238 | |

239 |
! Computing B^prim: |

240 |
BprimT=0.d0 |

241 |
ScanCoord=>Coordinate |

242 |
I=0 |

243 |
DO WHILE (Associated(ScanCoord%next)) |

244 |
I=I+1 |

245 |
SELECT CASE (ScanCoord%Type) |

246 |
CASE ('BOND') |

247 |
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,Geom,BprimT(1,I)) |

248 |
CASE ('ANGLE') |

249 |
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,Geom,BprimT(1,I)) |

250 |
CASE ('DIHEDRAL') |

251 |
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I)) |

252 |
END SELECT |

253 |
ScanCoord => ScanCoord%next |

254 |
END DO |

255 | |

256 |
! if (debug) THEN |

257 |
! WRITE(*,*) "PFL DBG, I,NPrim,BPrimT",I,NPrim |

258 |
! DO I=1,NPrim |

259 |
! WRITE(*,'(50(1X,F12.6))') BPrimT(:,I) |

260 |
! END DO |

261 |
! END IF |

262 | |

263 |
BMat_BakerT = 0.d0 |

264 |
DO I=1,NCoord |

265 |
DO J=1,NPrim |

266 |
!BprimT is transpose of B^prim. |

267 |
!B = UMat^T B^prim, B^T = (B^prim)^T UMat |

268 |
BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I) |

269 |
END DO |

270 |
END DO |

271 | |

272 |
BBT = 0.d0 |

273 |
DO I=1,NCoord |

274 |
DO J=1,3*Nat ! BBT(:,I) forms BB^T |

275 |
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) |

276 |
END DO |

277 |
END DO |

278 | |

279 |
BBT_inv = 0.d0 |

280 | |

281 |
!Print *, 'NCoord=', NCoord |

282 |
!Print *, 'BBT=' |

283 |
DO J=1,NCoord |

284 |
! WRITE(*,'(50(1X,F12.6))') BBT(:,J) |

285 |
END DO |

286 |
!Print *, 'End of BBT' |

287 |
! GenInv is in Mat_util.f90 |

288 |
Call GenInv(NCoord,BBT,BBT_inv,NCoord) |

289 |
! ALLOCATE(V(NCoord,NCoord)) |

290 |
! tol = 1e-12 |

291 |
! ! BBT is destroyed by GINVSE |

292 |
! Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT) |

293 |
! DEALLOCATE(V) |

294 |
! IF(Fail) Then |

295 |
! WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90" |

296 |
! STOP |

297 |
! END IF |

298 | |

299 |
!Print *, 'BBT_inv=' |

300 |
DO J=1,NCoord |

301 |
! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J) |

302 |
! Print *, BBT_inv(:,J) |

303 |
END DO |

304 |
!Print *, 'End of BBT_inv' |

305 | |

306 |
! Calculation of (B^T)^-1 = (BB^T)^-1B: |

307 |
BTransInv_local = 0.d0 |

308 |
DO I=1,3*Nat |

309 |
DO J=1,NCoord |

310 |
BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J) |

311 |
END DO |

312 |
END DO |

313 | |

314 |
!Print *, 'BMat_BakerT=' |

315 |
DO J=1,3*Nat |

316 |
!WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J) |

317 |
END DO |

318 |
!Print *, 'End of BMat_BakerT' |

319 | |

320 |
WRITE(*,*) "BTransInv_local Trans corresponding to x_k,y_k,z_k" |

321 |
DO J=1,3*Nat |

322 |
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) |

323 |
END DO |

324 | |

325 | |

326 |
END IF |

327 | |

328 |
DeltaX_int(:) = IntCoord(:)-IntCoord_k(:) |

329 | |

330 |
X_cart_k(1,:) = x_k(:) |

331 |
X_cart_k(2,:) = y_k(:) |

332 |
X_cart_k(3,:) = z_k(:) |

333 | |

334 |
normDeltaX_int=0.d0 ! This is to be implemented. |

335 |
DO I=1, NCoord |

336 |
normDeltaX_int=normDeltaX_int+DeltaX_int(I)*DeltaX_int(I) |

337 |
END DO |

338 |
normDeltaX_int = sqrt(normDeltaX_int) |

339 |
! I just need initial value of normDeltaX_int to be .GT. 1d-11, |

340 |
! that is why it is initialized to 1.d0 |

341 |
!normDeltaX_int = 1.d0 |

342 | |

343 |
if (debug) WRITE(*,*) "Entering the loop" |

344 |
Counter = 0 |

345 |
DO While ((normDeltaX_int .GT. Crit) .AND. (Counter .LT. NbItMax)) |

346 |
!DO While (normDeltaX_int .GT. 1d-6) |

347 |
Counter = Counter + 1 |

348 |
DeltaX_cart = 0.d0 |

349 |
! if (normDeltaX_int.LE.1e-4) THEN |

350 |
! DeltaX_int=DeltaX_int*1e4 |

351 |
! END IF |

352 |
DO I=1,NCoord |

353 |
! below BTransInv_local(I,:) is used because actually we need Transpose of BTransInv_local. |

354 |
DeltaX_cart(:)=DeltaX_cart(:)+BTransInv_local(I,:)*DeltaX_int(I) |

355 |
END DO |

356 |
! if (normDeltaX_int.LE.1e-4) THEN |

357 |
! DeltaX_int=DeltaX_int*1e-4 |

358 |
! DeltaX_cart=DeltaX_cart*1e-4 |

359 |
! END IF |

360 | |

361 | |

362 |
if (debug) THEN |

363 |
WRITE(*,*) "Info for iteration:",counter |

364 |
Print *, 'DeltaX_int=' |

365 |
WRITE(*,'(50(1X,F12.8))') DeltaX_int |

366 |
Print *, 'DeltaX_cart,norm',sqrt(dot_product(DeltaX_cart,DeltaX_cart)) |

367 |
DO J=1,Nat |

368 |
WRITE(*,'(1X,A4,3(1X,F15.11),3(1X,D25.15))') AtName(J),DeltaX_cart(3*J-2:3*J),DeltaX_cart(3*J-2:3*J) |

369 |
END DO |

370 |
Print *, 'BTransInv_local Trans=' |

371 |
DO J=1,3*Nat |

372 |
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) |

373 |
END DO |

374 |
Print *, 'DeltaX_int=' |

375 |
WRITE(*,'(50(1X,F12.8))') DeltaX_int |

376 |
END IF |

377 | |

378 |
! Finite precision error correction step: |

379 |
! DO I=1,3*Nat |

380 |
! IF (DeltaX_cart(I) .NE. 0.d0) Then |

381 |
! IF(abs(DeltaX_cart(I)) .LT. 1d-10) Then |

382 |
! !Print *, 'Correcting DeltaX_cart(',I,')=', DeltaX_cart(I) |

383 |
! DeltaX_cart(I) = 0.d0 |

384 |
! END IF |

385 |
! END IF |

386 |
! END DO |

387 | |

388 |
! new cartesian coordinates: |

389 |
DO I=1, Nat |

390 |
X_cart_k(1,I) = X_cart_k(1,I) + DeltaX_cart(3*I-2) |

391 |
X_cart_k(2,I) = X_cart_k(2,I) + DeltaX_cart(3*I-1) |

392 |
X_cart_k(3,I) = X_cart_k(3,I) + DeltaX_cart(3*I) |

393 |
END DO |

394 | |

395 |
! Geom(1,:)=X_cart_k(1,1:Nat)/a0 |

396 |
! Geom(2,:)=X_cart_k(2,1:Nat)/a0 |

397 |
! Geom(3,:)=X_cart_k(3,1:Nat)/a0 |

398 |
Geom(1,:)=X_cart_k(1,1:Nat) |

399 |
Geom(2,:)=X_cart_k(2,1:Nat) |

400 |
Geom(3,:)=X_cart_k(3,1:Nat) |

401 | |

402 |
if (debug) THEN |

403 |
WRITE(*,*) 'GeomCart used to compute BPrim' |

404 |
DO I=1,Nat |

405 |
WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,Geom(1,I),Geom(2,i),Geom(3,i) |

406 |
END DO |

407 |
END IF |

408 | |

409 |
! Computing B^prim: |

410 |
BprimT=0.d0 |

411 |
ScanCoord=>Coordinate |

412 |
I=0 |

413 |
DO WHILE (Associated(ScanCoord%next)) |

414 |
I=I+1 |

415 |
SELECT CASE (ScanCoord%Type) |

416 |
CASE ('BOND') |

417 |
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,Geom,BprimT(1,I)) |

418 |
CASE ('ANGLE') |

419 |
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,Geom,BprimT(1,I)) |

420 |
CASE ('DIHEDRAL') |

421 |
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I)) |

422 |
END SELECT |

423 |
ScanCoord => ScanCoord%next |

424 |
END DO |

425 | |

426 |
if (debug) THEN |

427 |
WRITE(*,*) "Bprim " |

428 |
DO J=1,3*Nat |

429 |
WRITE(*,'(50(1X,F12.6))') BprimT(J,:) |

430 |
END DO |

431 |
END IF |

432 | |

433 | |

434 |
BMat_BakerT = 0.d0 |

435 |
DO I=1,NCoord |

436 |
DO J=1,NPrim |

437 |
!BprimT is transpose of B^prim. |

438 |
!B = UMat^T B^prim, B^T = (B^prim)^T UMat |

439 |
BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I) |

440 |
END DO |

441 |
END DO |

442 | |

443 |
! ADDED FROM HERE: |

444 |
! We now compute G=B(BT) matrix |

445 |
!IF (IOptt .EQ. 5000) Then |

446 |
!GMat=0.d0 |

447 |
!DO I=1,NPrim |

448 |
! DO J=1,3*Nat |

449 |
! GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0))) |

450 |
!END DO |

451 |
!END DO |

452 | |

453 |
! We diagonalize G |

454 |
!Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim) |

455 |
!Call Trie(NPrim,EigVal,EigVec,NPrim) |

456 | |

457 |
! We construct the new DzDc(b=baker) matrix |

458 |
! UMat is nonredundant vector set, i.e. set of eigenvectors of |

459 |
! BB^T corresponding to eigenvalues > zero. |

460 | |

461 |
!UMat=0.d0 |

462 |
!UMat_tmp=0.d0 |

463 |
!BMat_BakerT = 0.d0 |

464 |
!J=0 |

465 |
!DO I=1,NPrim |

466 |
! IF (abs(eigval(I))>=1e-6) THEN |

467 |
! J=J+1 |

468 |
! DO K=1,NPrim |

469 |
! BprimT is transpose of B^prim. |

470 |
! B = UMat^T B^prim, B^T = (B^prim)^T UMat |

471 |
!BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I) |

472 |
! END DO |

473 |
!IF(J .GT. 3*Nat-6) THEN |

474 |
!WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. & |

475 |
! 1e-6 (=UMat) (=' ,J,') exceeded 3*Nat-6=',3*Nat-6, & |

476 |
! 'Stopping the calculation.' |

477 |
! STOP |

478 |
! END IF |

479 |
! UMat_tmp(:,J) = Eigvec(:,I) |

480 |
! END IF |

481 |
!END DO |

482 | |

483 |
! THIS IS TO BE CORRECTED, A special case for debugging C2H3F case: |

484 |
!J=0 |

485 |
!DO I=1, 3*Nat-6 |

486 |
! IF ((I .NE. 6) .AND. (I .NE. 7) .AND. (I .NE. 9)) Then |

487 |
! J=J+1 |

488 |
! Print *, 'J=', J, ' I=', I |

489 |
!UMat(:,J) = UMat_tmp(:,I) |

490 |
! END IF |

491 |
! END DO |

492 |
!BMat_BakerT=0.d0 |

493 |
! DO I=1,NCoord |

494 |
! DO J=1,NPrim |

495 |
! B = UMat^T B^prim, B^T = (B^prim)^T UMat |

496 |
! BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat(J,I) |

497 |
! END DO |

498 |
!END DO |

499 |
! TO BE CORRECTED, ENDS HERE. |

500 | |

501 |
!END IF |

502 |
! END of the new changes. |

503 | |

504 |
BBT = 0.d0 |

505 |
DO I=1,NCoord |

506 |
DO J=1,3*Nat ! BBT(:,I) forms BB^T |

507 |
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) |

508 |
END DO |

509 |
END DO |

510 | |

511 |
BBT_inv = 0.d0 |

512 | |

513 |
!Print *, 'NCoord=', NCoord |

514 |
!Print *, 'BBT=' |

515 |
DO J=1,NCoord |

516 |
! WRITE(*,'(50(1X,F12.6))') BBT(:,J) |

517 |
END DO |

518 |
!Print *, 'End of BBT' |

519 |
! GenInv is in Mat_util.f90 |

520 |
Call GenInv(NCoord,BBT,BBT_inv,NCoord) |

521 |
! if (debug) THEN |

522 |
! WRITE(*,*) 'BBT_inv by GenInv' |

523 |
! DO J=1,NCoord |

524 |
! WRITE(*,'(50(1X,F12.6))') BBT_inv(:,J) |

525 |
! END DO |

526 |
! END IF |

527 | |

528 |
! ALLOCATE(V(NCoord,NCoord)) |

529 |
! tol = 1e-12 |

530 |
! ! BBT is destroyed by GINVSE |

531 |
! Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT) |

532 |
! DEALLOCATE(V) |

533 |
! IF(Fail) Then |

534 |
! WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90" |

535 |
! STOP |

536 |
! END IF |

537 | |

538 |
! if (debug) THEN |

539 |
! WRITE(*,*) 'BBT_inv by Genvse' |

540 |
! DO J=1,NCoord |

541 |
! WRITE(*,'(50(1X,F12.6))') BBT_inv(:,J) |

542 |
! END DO |

543 |
! END IF |

544 | |

545 |
!Print *, 'End of BBT_inv' |

546 | |

547 |
! Calculation of (B^T)^-1 = (BB^T)^-1B: |

548 |
BTransInv_local = 0.d0 |

549 |
DO I=1,3*Nat |

550 |
DO J=1,NCoord |

551 |
BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J) |

552 |
END DO |

553 |
END DO |

554 | |

555 |
! if (debug) THEN |

556 |
! Print *, 'BMat_BakerT=' |

557 |
! DO J=1,3*Nat |

558 |
! WRITE(*,'(50(1X,F12.6))') BMat_BakerT(:,J) |

559 |
! END DO |

560 |
! Print *, 'End of BMat_BakerT' |

561 |
! END IF |

562 | |

563 |
! if (debug) THEN |

564 |
! Print *, 'BTransInv_local Trans=' |

565 |
! DO J=1,3*Nat |

566 |
! WRITE(*,'(50(1X,F16.6))') BTransInv_local(:,J) |

567 |
! END DO |

568 |
! Print *, 'End of BTransInv_local Trans' |

569 |
! END IF |

570 | |

571 |
! New cartesian cooordinates: |

572 |
x(:) = X_cart_k(1,:) |

573 |
y(:) = X_cart_k(2,:) |

574 |
z(:) = X_cart_k(3,:) |

575 | |

576 | |

577 |
Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive_t,XPrimRef) |

578 | |

579 |
! I=0 ! index for the NPrim (NPrim is the number of |

580 |
! ! primitive coordinates). |

581 |
! Xprimitive_t=0.d0 |

582 |
! ScanCoord=>Coordinate |

583 |
! DO WHILE (Associated(ScanCoord%next)) |

584 |
! I=I+1 |

585 |
! SELECT CASE (ScanCoord%Type) |

586 |
! CASE ('BOND') |

587 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) |

588 |
! Xprimitive_t(I) = Norm2 |

589 |
! !Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2, '=', Xprimitive_t(I) |

590 |
! CASE ('ANGLE') |

591 |
! Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1) |

592 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) |

593 |
! Xprimitive_t(I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |

594 |
! !Print *, I, ':', ScanCoord%At1, '-', ScanCoord%At2,'-',ScanCoord%At3, '=', Xprimitive_t(I)*180./Pi |

595 |
! CASE ('DIHEDRAL') |

596 |
! Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2) |

597 |
! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1) |

598 |
! Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3) |

599 |
! Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5) |

600 |
! Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4) |

601 |
! !!! PFL 28 Feb 2008 Test to be sure that angle does not change by more than Pi |

602 |
! ! Xprimitive_t(I)=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)*Pi/180. |

603 |
! DiheTmp=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2) |

604 |
! Xprimitive_t(I) = DiheTmp*Pi/180. |

605 |
! ! We treat large dihedral angles differently as +180=-180 mathematically and physically |

606 |
! ! but this causes lots of troubles when converting baker to cart. |

607 |
! ! So we ensure that large dihedral angles always have the same sign |

608 |
! if (abs(ScanCoord%SignDihedral).NE.1) THEN |

609 |
! ScanCoord%SignDihedral=Int(Sign(1.d0,DiheTmp)) |

610 |
! ELSE |

611 |
! If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN |

612 |
! Xprimitive_t(I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi |

613 |
! END IF |

614 |
! END IF |

615 | |

616 | |

617 | |

618 |
! !Print *, I, ':', ScanCoord%At1, '-',ScanCoord%At2,'-',ScanCoord%At3,'-',& |

619 |
! !ScanCoord%At4,'=',Xprimitive_t(I)*180./Pi |

620 |
! END SELECT |

621 |
! ScanCoord => ScanCoord%next |

622 |
! END DO ! matches DO WHILE (Associated(ScanCoord%next)) |

623 | |

624 |
! if (debug) THEN |

625 |
! WRITE(*,*) "Xprimitive_t" |

626 |
! WRITE(*,'(50(1X,F12.6))') Xprimitive_t |

627 |
! END IF |

628 | |

629 |
IntCoord_k=0.d0 |

630 |
DO I=1, NPrim |

631 |
! Transpose of UMat is needed below, that is why UMat(I,:). |

632 |
IntCoord_k(:)=IntCoord_k(:)+UMat_local(I,:)*Xprimitive_t(I) |

633 |
END DO |

634 | |

635 |
if (debug) THEN |

636 |
WRITE(*,*) "New Xint (IntCoord_k)" |

637 |
WRITE(*,'(50(1X,F12.8))') IntCoord_k |

638 |
WRITE(*,*) "Target (IntCoord)" |

639 |
WRITE(*,'(50(1X,F12.8))') IntCoord |

640 | |

641 |
END IF |

642 | |

643 | |

644 |
! Computing DeltaX_int |

645 | |

646 |
DeltaX_int(:) = IntCoord(:)-IntCoord_k(:) |

647 | |

648 | |

649 | |

650 | |

651 |
! norm2 of DeltaX_int is being calculated here: |

652 |
normDeltaX_int=0.d0 |

653 |
DO I=1, NCoord |

654 |
normDeltaX_int=normDeltaX_int+DeltaX_int(I)*DeltaX_int(I) |

655 |
END DO |

656 |
normDeltaX_int = sqrt(normDeltaX_int) |

657 | |

658 |
if (debug) THEN |

659 |
WRITE(*,*) "New Delta_Xint (deltaX_int)" |

660 |
WRITE(*,'(50(1X,F12.6))') DeltaX_int |

661 |
WRITE(*,*) "Norm:",normDeltaX_int |

662 |
END IF |

663 | |

664 |
!IF (Counter .GE. 400) THEN |

665 |
!Print *, 'Counter=', Counter, 'normDeltaX_int=', normDeltaX_int |

666 |
!Print *, 'DeltaX_int=' |

667 |
!WRITE(*,'(50(1X,F8.2))') DeltaX_int |

668 |
!END IF |

669 | |

670 |
IF (Mod(Counter,1).EQ.0) THEN |

671 |
!WRITE(*,*) Nat |

672 |
!WRITE(*,*) 'GeomCart in ConvertBakerInternal_cart' |

673 |
DO I=1,Nat |

674 |
! WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,X_Cart_k(1:3,I) |

675 |
END DO |

676 |
END IF |

677 | |

678 |
!call two_norm(DeltaX_int,normDeltaX_int,3*Nat) |

679 |
END DO !matches DO While (normDeltaX_int .GT. 1d-11) |

680 | |

681 |
IF (debug) THEN |

682 |
WRITE(*,*) 'GeomCart in ConvertBakerInternal_cart' |

683 |
DO I=1,Nat |

684 |
WRITE(*,'(1X,A,I4,3(1X,F12.4))') AtName(I),I,X_Cart_k(1:3,I) |

685 |
END DO |

686 |
END IF |

687 | |

688 |
if (debug) THEN |

689 |
WRITE(*,*) "Bmat_bakerT" |

690 |
DO J=1,NCoord |

691 |
WRITE(*,'(50(1X,F12.6))') BMat_BakerT(:,J) |

692 |
END DO |

693 |
END IF |

694 | |

695 | |

696 | |

697 |
!IF (IOptt .GE. 2) THEN |

698 |
! Print *, 'Counter=', Counter |

699 |
!END IF |

700 |
IF (Counter .GE. NbItMax) Then |

701 |
Print *, 'Counter GE 500, Stopping, normDeltaX_int=', normDeltaX_int |

702 |
STOP |

703 |
END IF |

704 | |

705 | |

706 |
x(:) = X_cart_k(1,:) |

707 |
y(:) = X_cart_k(2,:) |

708 |
z(:) = X_cart_k(3,:) |

709 | |

710 | |

711 |
! call CalcRmsd(Nat,x_k,y_k,z_k,x,y,z,MRot,rmsd,.TRUE.,.TRUE.) |

712 | |

713 | |

714 |
if (debug) WRITE(*,*) "COnverted in ",counter," iterations" |

715 | |

716 |
DEALLOCATE(Geom,DeltaX_int,DeltaX_cart) |

717 |
DEALLOCATE(BprimT,BBT,BBT_inv,UMat_tmp) |

718 |
DEALLOCATE(X_cart_k) |

719 |
DEALLOCATE(GMat,EigVal,EigVec) |

720 | |

721 |
XPrim=Xprimitive_t |

722 | |

723 |
if (debug) WRITE(*,*) '=====ConvertBakerInternal_cart Over=====' |

724 | |

725 |
END SUBROUTINE ConvertBakerInternal_cart |