r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# flake8: noqa

2import numpy as np

3from numpy import linalg

4from ase import units

6# Three variables extracted from what used to be endless repetitions below.

7Ax = np.array([[1, 0, 0, -1, 0, 0, 0, 0, 0],

8 [0, 1, 0, 0, -1, 0, 0, 0, 0],

9 [0, 0, 1, 0, 0, -1, 0, 0, 0],

10 [0, 0, 0, -1, 0, 0, 1, 0, 0],

11 [0, 0, 0, 0, -1, 0, 0, 1, 0],

12 [0, 0, 0, 0, 0, -1, 0, 0, 1]])

13Bx = np.array([[1, 0, 0, -1, 0, 0],

14 [0, 1, 0, 0, -1, 0],

15 [0, 0, 1, 0, 0, -1]])

16Mx = Bx

20class Morse:

21 def __init__(self, atomi, atomj, D, alpha, r0):

22 self.atomi = atomi

23 self.atomj = atomj

24 self.D = D

25 self.alpha = alpha

26 self.r0 = r0

27 self.r = None

29class Bond:

30 def __init__(self, atomi, atomj, k, b0,

31 alpha=None, rref=None):

32 self.atomi = atomi

33 self.atomj = atomj

34 self.k = k

35 self.b0 = b0

36 self.alpha = alpha

37 self.rref = rref

38 self.b = None

40class Angle:

41 def __init__(self, atomi, atomj, atomk, k, a0, cos=False,

42 alpha=None, rref=None):

43 self.atomi = atomi

44 self.atomj = atomj

45 self.atomk = atomk

46 self.k = k

47 self.a0 = a0

48 self.cos = cos

49 self.alpha = alpha

50 self.rref = rref

51 self.a = None

53class Dihedral:

54 def __init__(self, atomi, atomj, atomk, atoml, k, d0=None, n=None,

55 alpha=None, rref=None):

56 self.atomi = atomi

57 self.atomj = atomj

58 self.atomk = atomk

59 self.atoml = atoml

60 self.k = k

61 self.d0 = d0

62 self.n = n

63 self.alpha = alpha

64 self.rref = rref

65 self.d = None

67class VdW:

68 def __init__(self, atomi, atomj, epsilonij=None, sigmaij=None, rminij=None,

69 Aij=None, Bij=None, epsiloni=None, epsilonj=None,

70 sigmai=None, sigmaj=None, rmini=None, rminj=None, scale=1.0):

71 self.atomi = atomi

72 self.atomj = atomj

73 if epsilonij is not None:

74 if sigmaij is not None:

75 self.Aij = scale * 4.0 * epsilonij * sigmaij**12

76 self.Bij = scale * 4.0 * epsilonij * sigmaij**6 * scale

77 elif rminij is not None:

78 self.Aij = scale * epsilonij * rminij**12

79 self.Bij = scale * 2.0 * epsilonij * rminij**6

80 else:

81 raise NotImplementedError("not implemented combination"

82 "of vdW parameters.")

83 elif Aij is not None and Bij is not None:

84 self.Aij = scale * Aij

85 self.Bij = scale * Bij

86 elif epsiloni is not None and epsilonj is not None:

87 if sigmai is not None and sigmaj is not None:

88 self.Aij = ( scale * 4.0 * np.sqrt(epsiloni * epsilonj)

89 * ((sigmai + sigmaj) / 2.0)**12 )

90 self.Bij = ( scale * 2.0 * np.sqrt(epsiloni * epsilonj)

91 * ((sigmai + sigmaj) / 2.0)**6 )

92 elif rmini is not None and rminj is not None:

93 self.Aij = ( scale * np.sqrt(epsiloni * epsilonj)

94 * ((rmini + rminj) / 2.0)**12 )

95 self.Bij = ( scale * 2.0 * np.sqrt(epsiloni * epsilonj)

96 * ((rmini + rminj) / 2.0)**6 )

97 else:

98 raise NotImplementedError("not implemented combination"

99 "of vdW parameters.")

100 self.r = None

102class Coulomb:

103 def __init__(self, atomi, atomj, chargeij=None,

104 chargei=None, chargej=None, scale=1.0):

105 self.atomi = atomi

106 self.atomj = atomj

107 if chargeij is not None:

108 self.chargeij = ( scale * chargeij * 8.9875517873681764e9

109 * units.m * units.J / units.C / units.C )

110 elif chargei is not None and chargej is not None:

111 self.chargeij = ( scale * chargei * chargej * 8.9875517873681764e9

112 * units.m * units.J / units.C / units.C )

113 else:

114 raise NotImplementedError("not implemented combination"

115 "of Coulomb parameters.")

116 self.r = None

118def get_morse_potential_eta(atoms, morse):

119 i = morse.atomi

120 j = morse.atomj

122 rij = rel_pos_pbc(atoms, i, j)

123 dij = linalg.norm(rij)

125 if dij > morse.r0:

126 exp = np.exp(-morse.alpha*(dij-morse.r0))

127 eta = 1.0 - (1.0 - exp)**2

128 else:

129 eta = 1.0

131 return eta

133def get_morse_potential_value(atoms, morse):

134 i = morse.atomi

135 j = morse.atomj

137 rij = rel_pos_pbc(atoms, i, j)

138 dij = linalg.norm(rij)

140 exp = np.exp(-morse.alpha*(dij-morse.r0))

142 v = morse.D*(1.0-exp)**2

144 morse.r = dij

146 return i, j, v

149 i = morse.atomi

150 j = morse.atomj

152 rij = rel_pos_pbc(atoms, i, j)

153 dij = linalg.norm(rij)

154 eij = rij/dij

156 exp = np.exp(-morse.alpha*(dij-morse.r0))

158 gr = 2.0*morse.D*morse.alpha*exp*(1.0-exp)*eij

160 gx = np.dot(Mx.T, gr)

162 morse.r = dij

164 return i, j, gx

166def get_morse_potential_hessian(atoms, morse, spectral=False):

167 i = morse.atomi

168 j = morse.atomj

170 rij = rel_pos_pbc(atoms, i, j)

171 dij = linalg.norm(rij)

172 eij = rij/dij

174 Pij = np.tensordot(eij,eij,axes=0)

175 Qij = np.eye(3)-Pij

177 exp = np.exp(-morse.alpha*(dij-morse.r0))

179 Hr = ( 2.0*morse.D*morse.alpha*exp*(morse.alpha*(2.0*exp-1.0)*Pij

180 + (1.0-exp)/dij*Qij) )

182 Hx = np.dot(Mx.T, np.dot(Hr, Mx))

184 if spectral:

185 eigvals, eigvecs = linalg.eigh(Hx)

186 D = np.diag(np.abs(eigvals))

187 U = eigvecs

188 Hx = np.dot(U,np.dot(D,np.transpose(U)))

190 morse.r = dij

192 return i, j, Hx

194def get_morse_potential_reduced_hessian(atoms, morse):

195 i = morse.atomi

196 j = morse.atomj

198 rij = rel_pos_pbc(atoms, i, j)

199 dij = linalg.norm(rij)

200 eij = rij/dij

202 Pij = np.tensordot(eij,eij,axes=0)

204 exp = np.exp(-morse.alpha*(dij-morse.r0))

206 Hr = np.abs(2.0*morse.D*morse.alpha**2*exp*(2.0*exp-1.0))*Pij

208 Hx = np.dot(Mx.T, np.dot(Hr, Mx))

210 morse.r = dij

212 return i, j, Hx

214def get_bond_potential_value(atoms, bond):

215 i = bond.atomi

216 j = bond.atomj

218 rij = rel_pos_pbc(atoms, i, j)

219 dij = linalg.norm(rij)

221 v = 0.5*bond.k*(dij-bond.b0)**2

223 bond.b = dij

225 return i, j, v

228 i = bond.atomi

229 j = bond.atomj

231 rij = rel_pos_pbc(atoms, i, j)

232 dij = linalg.norm(rij)

233 eij = rij/dij

235 gr = bond.k*(dij-bond.b0)*eij

237 gx = np.dot(Bx.T, gr)

239 bond.b = dij

241 return i, j, gx

243def get_bond_potential_hessian(atoms, bond, morses=None, spectral=False):

244 i = bond.atomi

245 j = bond.atomj

247 rij = rel_pos_pbc(atoms, i, j)

248 dij = linalg.norm(rij)

249 eij = rij/dij

251 Pij = np.tensordot(eij,eij,axes=0)

252 Qij = np.eye(3)-Pij

254 Hr = bond.k*Pij+bond.k*(dij-bond.b0)/dij*Qij

256 if bond.alpha is not None:

257 Hr *= np.exp(bond.alpha[0]*(bond.rref[0]**2-dij**2))

259 if morses is not None:

260 for m in range(len(morses)):

261 if ( morses[m].atomi == i or

262 morses[m].atomi == j ):

263 Hr *= get_morse_potential_eta(atoms, morses[m])

264 elif ( morses[m].atomj == i or

265 morses[m].atomj == j ):

266 Hr *= get_morse_potential_eta(atoms, morses[m])

268 Hx = np.dot(Bx.T, np.dot(Hr, Bx))

270 if spectral:

271 eigvals, eigvecs = linalg.eigh(Hx)

272 D = np.diag(np.abs(eigvals))

273 U = eigvecs

274 Hx = np.dot(U,np.dot(D,np.transpose(U)))

276 bond.b = dij

278 return i, j, Hx

280def get_bond_potential_reduced_hessian(atoms, bond, morses=None):

281 i = bond.atomi

282 j = bond.atomj

284 rij = rel_pos_pbc(atoms, i, j)

285 dij = linalg.norm(rij)

286 eij = rij/dij

288 Pij = np.tensordot(eij,eij,axes=0)

290 Hr = bond.k*Pij

292 if bond.alpha is not None:

293 Hr *= np.exp(bond.alpha[0]*(bond.rref[0]**2-dij**2))

295 if morses is not None:

296 for m in range(len(morses)):

297 if ( morses[m].atomi == i or

298 morses[m].atomi == j ):

299 Hr *= get_morse_potential_eta(atoms, morses[m])

300 elif ( morses[m].atomj == i or

301 morses[m].atomj == j ):

302 Hr *= get_morse_potential_eta(atoms, morses[m])

304 Hx = np.dot(Bx.T, np.dot(Hr, Bx))

306 bond.b = dij

308 return i, j, Hx

310def get_bond_potential_reduced_hessian_test(atoms, bond):

312 i, j, v = get_bond_potential_value(atoms, bond)

313 i, j, gx = get_bond_potential_gradient(atoms, bond)

315 Hx = np.tensordot(gx,gx,axes=0)/v/2.0

317 return i, j, Hx

319def get_angle_potential_value(atoms, angle):

321 i = angle.atomi

322 j = angle.atomj

323 k = angle.atomk

325 rij = rel_pos_pbc(atoms, i, j)

326 dij = linalg.norm(rij)

327 eij = rij/dij

328 rkj = rel_pos_pbc(atoms, k, j)

329 dkj = linalg.norm(rkj)

330 ekj = rkj/dkj

331 eijekj = np.dot(eij, ekj)

332 if np.abs(eijekj) > 1.0:

333 eijekj = np.sign(eijekj)

335 a = np.arccos(eijekj)

337 if angle.cos:

338 da = np.cos(a)-np.cos(angle.a0)

339 else:

340 da = a-angle.a0

341 da = da - np.around(da / np.pi) * np.pi

343 v = 0.5*angle.k*da**2

345 angle.a = a

347 return i, j, k, v

350 i = angle.atomi

351 j = angle.atomj

352 k = angle.atomk

354 rij = rel_pos_pbc(atoms, i, j)

355 dij = linalg.norm(rij)

356 eij = rij/dij

357 rkj = rel_pos_pbc(atoms, k, j)

358 dkj = linalg.norm(rkj)

359 ekj = rkj/dkj

360 eijekj = np.dot(eij, ekj)

361 if np.abs(eijekj) > 1.0:

362 eijekj = np.sign(eijekj)

364 a = np.arccos(eijekj)

365 if angle.cos:

366 da = np.cos(a)-np.cos(angle.a0)

367 else:

368 da = a-angle.a0

369 da = da - np.around(da / np.pi) * np.pi

370 sina = np.sin(a)

372 Pij = np.tensordot(eij,eij,axes=0)

373 Qij = np.eye(3)-Pij

374 Pkj = np.tensordot(ekj,ekj,axes=0)

375 Qkj = np.eye(3)-Pkj

377 gr = np.zeros(6)

378 if angle.cos:

379 gr[0:3] = angle.k*da/dij*np.dot(Qij,ekj)

380 gr[3:6] = angle.k*da/dkj*np.dot(Qkj,eij)

381 elif np.abs(sina) > 0.001:

382 gr[0:3] = -angle.k*da/sina/dij*np.dot(Qij,ekj)

383 gr[3:6] = -angle.k*da/sina/dkj*np.dot(Qkj,eij)

385 gx = np.dot(Ax.T, gr)

387 angle.a = a

389 return i, j, k, gx

391def get_angle_potential_hessian(atoms, angle, morses=None, spectral=False):

392 i = angle.atomi

393 j = angle.atomj

394 k = angle.atomk

396 rij = rel_pos_pbc(atoms, i, j)

397 dij = linalg.norm(rij)

398 dij2 = dij*dij

399 eij = rij/dij

400 rkj = rel_pos_pbc(atoms, k, j)

401 dkj = linalg.norm(rkj)

402 dkj2 = dkj*dkj

403 ekj = rkj/dkj

404 dijdkj = dij*dkj

405 eijekj = np.dot(eij, ekj)

406 if np.abs(eijekj) > 1.0:

407 eijekj = np.sign(eijekj)

409 a = np.arccos(eijekj)

410 if angle.cos:

411 da = np.cos(a)-np.cos(angle.a0)

412 cosa0 = np.cos(angle.a0)

413 else:

414 da = a-angle.a0

415 da = da - np.around(da / np.pi) * np.pi

416 sina = np.sin(a)

417 cosa = np.cos(a)

418 ctga = cosa/sina

420 Pij = np.tensordot(eij,eij,axes=0)

421 Qij = np.eye(3)-Pij

422 Pkj = np.tensordot(ekj,ekj,axes=0)

423 Qkj = np.eye(3)-Pkj

424 Pik = np.tensordot(eij,ekj,axes=0)

425 Pki = np.tensordot(ekj,eij,axes=0)

426 P = np.eye(3)*eijekj

428 QijPkjQij = np.dot(Qij, np.dot(Pkj, Qij))

429 QijPkiQkj = np.dot(Qij, np.dot(Pki, Qkj))

430 QkjPijQkj = np.dot(Qkj, np.dot(Pij, Qkj))

432 Hr = np.zeros((6,6))

433 if angle.cos and np.abs(sina) > 0.001:

434 factor = 1.0-2.0*cosa*cosa+cosa*cosa0

435 Hr[0:3,0:3] = ( angle.k*(factor*QijPkjQij/sina

436 - sina*da*(-ctga*QijPkjQij/sina+np.dot(Qij, Pki)

437 -np.dot(Pij, Pki)*2.0+(Pik+P)))/sina/dij2 )

438 Hr[0:3,3:6] = ( angle.k*(factor*QijPkiQkj/sina

439 - sina*da*(-ctga*QijPkiQkj/sina

440 -np.dot(Qij, Qkj)))/sina/dijdkj )

441 Hr[3:6,0:3] = Hr[0:3,3:6].T

442 Hr[3:6,3:6] = ( angle.k*(factor*QkjPijQkj/sina

443 - sina*da*(-ctga*QkjPijQkj/sina

444 +np.dot(Qkj, Pik)-np.dot(Pkj, Pik)

445 *2.0+(Pki+P)))/sina/dkj2 )

446 elif np.abs(sina) > 0.001:

447 Hr[0:3,0:3] = ( angle.k*(QijPkjQij/sina

448 + da*(-ctga*QijPkjQij/sina+np.dot(Qij, Pki)

449 -np.dot(Pij, Pki)*2.0+(Pik+P)))/sina/dij2 )

450 Hr[0:3,3:6] = ( angle.k*(QijPkiQkj/sina

451 + da*(-ctga*QijPkiQkj/sina

452 -np.dot(Qij, Qkj)))/sina/dijdkj )

453 Hr[3:6,0:3] = Hr[0:3,3:6].T

454 Hr[3:6,3:6] = ( angle.k*(QkjPijQkj/sina

455 + da*(-ctga*QkjPijQkj/sina

456 +np.dot(Qkj, Pik)-np.dot(Pkj, Pik)

457 *2.0+(Pki+P)))/sina/dkj2 )

459 if angle.alpha is not None:

460 Hr *= ( np.exp(angle.alpha[0]*(angle.rref[0]**2-dij**2))

461 *np.exp(angle.alpha[1]*(angle.rref[1]**2-dkj**2)) )

463 if morses is not None:

464 for m in range(len(morses)):

465 if ( morses[m].atomi == i or

466 morses[m].atomi == j or

467 morses[m].atomi == k ):

468 Hr *= get_morse_potential_eta(atoms, morses[m])

469 elif ( morses[m].atomj == i or

470 morses[m].atomj == j or

471 morses[m].atomj == k ):

472 Hr *= get_morse_potential_eta(atoms, morses[m])

474 Hx = np.dot(Ax.T, np.dot(Hr, Ax))

476 if spectral:

477 eigvals, eigvecs = linalg.eigh(Hx)

478 D = np.diag(np.abs(eigvals))

479 U = eigvecs

480 Hx = np.dot(U,np.dot(D,np.transpose(U)))

482 angle.a = a

484 return i, j, k, Hx

486def get_angle_potential_reduced_hessian(atoms, angle, morses=None):

487 i = angle.atomi

488 j = angle.atomj

489 k = angle.atomk

491 rij = rel_pos_pbc(atoms, i, j)

492 dij = linalg.norm(rij)

493 dij2 = dij*dij

494 eij = rij/dij

495 rkj = rel_pos_pbc(atoms, k, j)

496 dkj = linalg.norm(rkj)

497 dkj2 = dkj*dkj

498 ekj = rkj/dkj

499 dijdkj = dij*dkj

500 eijekj = np.dot(eij, ekj)

501 if np.abs(eijekj) > 1.0:

502 eijekj = np.sign(eijekj)

504 a = np.arccos(eijekj)

505 sina = np.sin(a)

506 sina2 = sina*sina

508 Pij = np.tensordot(eij,eij,axes=0)

509 Qij = np.eye(3)-Pij

510 Pkj = np.tensordot(ekj,ekj,axes=0)

511 Qkj = np.eye(3)-Pkj

512 Pki = np.tensordot(ekj,eij,axes=0)

514 Hr = np.zeros((6,6))

515 if np.abs(sina) > 0.001:

516 Hr[0:3,0:3] = np.dot(Qij, np.dot(Pkj, Qij))/dij2

517 Hr[0:3,3:6] = np.dot(Qij, np.dot(Pki, Qkj))/dijdkj

518 Hr[3:6,0:3] = Hr[0:3,3:6].T

519 Hr[3:6,3:6] = np.dot(Qkj, np.dot(Pij, Qkj))/dkj2

521 if angle.cos and np.abs(sina) > 0.001:

522 cosa = np.cos(a)

523 cosa0 = np.cos(angle.a0)

524 factor = np.abs(1.0-2.0*cosa*cosa+cosa*cosa0)

525 Hr = Hr*factor*angle.k/sina2

526 elif np.abs(sina) > 0.001:

527 Hr = Hr*angle.k/sina2

529 if angle.alpha is not None:

530 Hr *= ( np.exp(angle.alpha[0]*(angle.rref[0]**2-dij**2))

531 *np.exp(angle.alpha[1]*(angle.rref[1]**2-dkj**2)) )

533 if morses is not None:

534 for m in range(len(morses)):

535 if ( morses[m].atomi == i or

536 morses[m].atomi == j or

537 morses[m].atomi == k ):

538 Hr *= get_morse_potential_eta(atoms, morses[m])

539 elif ( morses[m].atomj == i or

540 morses[m].atomj == j or

541 morses[m].atomj == k ):

542 Hr *= get_morse_potential_eta(atoms, morses[m])

544 Hx=np.dot(Ax.T, np.dot(Hr, Ax))

546 angle.a = a

548 return i, j, k, Hx

550def get_angle_potential_reduced_hessian_test(atoms, angle):

551 i, j, k, v = get_angle_potential_value(atoms, angle)

552 i, j, k, gx = get_angle_potential_gradient(atoms, angle)

554 Hx = np.tensordot(gx,gx,axes=0)/v/2.0

556 return i, j, k, Hx

558def get_dihedral_potential_value(atoms, dihedral):

559 i = dihedral.atomi

560 j = dihedral.atomj

561 k = dihedral.atomk

562 l = dihedral.atoml

564 rij = rel_pos_pbc(atoms, i, j)

565 rkj = rel_pos_pbc(atoms, k, j)

566 rkl = rel_pos_pbc(atoms, k, l)

568 rmj = np.cross(rij, rkj)

569 dmj = linalg.norm(rmj)

570 emj = rmj/dmj

571 rnk = np.cross(rkj, rkl)

572 dnk = linalg.norm(rnk)

573 enk = rnk/dnk

574 emjenk = np.dot(emj, enk)

575 if np.abs(emjenk) > 1.0:

576 emjenk = np.sign(emjenk)

578 d = np.sign(np.dot(rkj, np.cross(rmj, rnk)))*np.arccos(emjenk)

580 if dihedral.d0 is None:

581 v = 0.5*dihedral.k*(1.0 - np.cos(2.0 * d))

582 else:

583 dd = d-dihedral.d0

584 dd = dd - np.around(dd / np.pi / 2.0) * np.pi * 2.0

585 if dihedral.n is None:

586 v = 0.5*dihedral.k*dd**2

587 else:

588 v = dihedral.k*(1.0 + np.cos(dihedral.n*d - dihedral.d0))

590 dihedral.d = d

592 return i, j, k, l, v

595 i = dihedral.atomi

596 j = dihedral.atomj

597 k = dihedral.atomk

598 l = dihedral.atoml

600 rij = rel_pos_pbc(atoms, i, j)

601 rkj = rel_pos_pbc(atoms, k, j)

602 dkj = linalg.norm(rkj)

603 dkj2 = dkj*dkj

604 rkl = rel_pos_pbc(atoms, k, l)

606 rijrkj = np.dot(rij, rkj)

607 rkjrkl = np.dot(rkj, rkl)

609 rmj = np.cross(rij, rkj)

610 dmj = linalg.norm(rmj)

611 dmj2 = dmj*dmj

612 emj = rmj/dmj

613 rnk = np.cross(rkj, rkl)

614 dnk = linalg.norm(rnk)

615 dnk2 = dnk*dnk

616 enk = rnk/dnk

617 emjenk = np.dot(emj, enk)

618 if np.abs(emjenk) > 1.0:

619 emjenk = np.sign(emjenk)

621 dddri = dkj/dmj2*rmj

622 dddrl = -dkj/dnk2*rnk

624 gx = np.zeros(12)

626 gx[0:3] = dddri

627 gx[3:6] = (rijrkj/dkj2-1.0)*dddri-rkjrkl/dkj2*dddrl

628 gx[6:9] = (rkjrkl/dkj2-1.0)*dddrl-rijrkj/dkj2*dddri

629 gx[9:12] = dddrl

631 d = np.sign(np.dot(rkj, np.cross(rmj, rnk)))*np.arccos(emjenk)

633 if dihedral.d0 is None:

634 gx *= dihedral.k*np.sin(2.0 * d)

635 else:

636 dd = d-dihedral.d0

637 dd = dd - np.around(dd / np.pi / 2.0) * np.pi * 2.0

638 if dihedral.n is None:

639 gx *= dihedral.k*dd

640 else:

641 gx *= -dihedral.k*dihedral.n*np.sin(dihedral.n*d - dihedral.d0)

643 dihedral.d = d

645 return i, j, k, l, gx

647def get_dihedral_potential_hessian(atoms, dihedral, morses=None,

648 spectral=False):

649 eps = 0.000001

653 Hx = np.zeros((12,12))

655 dihedral_eps = Dihedral(dihedral.atomi, dihedral.atomj,

656 dihedral.atomk, dihedral.atoml,

657 dihedral.k, dihedral.d0, dihedral.n)

658 indx = [3*i, 3*i+1, 3*i+2,

659 3*j, 3*j+1, 3*j+2,

660 3*k, 3*k+1, 3*k+2,

661 3*l, 3*l+1, 3*l+2]

662 for x in range(12):

663 a = atoms.copy()

664 positions = np.reshape(a.get_positions(),-1)

665 positions[indx[x]] += eps

666 a.set_positions(np.reshape(positions, (len(a),3)))

668 for y in range(12):

669 Hx[x,y] += 0.5*(geps[y]-g[y])/eps

670 Hx[y,x] += 0.5*(geps[y]-g[y])/eps

672 if dihedral.alpha is not None:

673 rij = rel_pos_pbc(atoms, i, j)

674 dij = linalg.norm(rij)

675 rkj = rel_pos_pbc(atoms, k, j)

676 dkj = linalg.norm(rkj)

677 rkl = rel_pos_pbc(atoms, k, l)

678 dkl = linalg.norm(rkl)

679 Hx *= ( np.exp(dihedral.alpha[0]*(dihedral.rref[0]**2-dij**2))

680 *np.exp(dihedral.alpha[1]*(dihedral.rref[1]**2-dkj**2))

681 *np.exp(dihedral.alpha[2]*(dihedral.rref[2]**2-dkl**2)) )

683 if morses is not None:

684 for m in range(len(morses)):

685 if ( morses[m].atomi == i or

686 morses[m].atomi == j or

687 morses[m].atomi == k or

688 morses[m].atomi == l ):

689 Hx *= get_morse_potential_eta(atoms, morses[m])

690 elif ( morses[m].atomj == i or

691 morses[m].atomj == j or

692 morses[m].atomj == k or

693 morses[m].atomj == l ):

694 Hx *= get_morse_potential_eta(atoms, morses[m])

696 if spectral:

697 eigvals, eigvecs = linalg.eigh(Hx)

698 D = np.diag(np.abs(eigvals))

699 U = eigvecs

700 Hx = np.dot(U,np.dot(D,np.transpose(U)))

702 return i, j, k, l, Hx

704def get_dihedral_potential_reduced_hessian(atoms, dihedral, morses=None):

705 i = dihedral.atomi

706 j = dihedral.atomj

707 k = dihedral.atomk

708 l = dihedral.atoml

710 rij = rel_pos_pbc(atoms, i, j)

711 rkj = rel_pos_pbc(atoms, k, j)

712 dkj = linalg.norm(rkj)

713 dkj2 = dkj*dkj

714 rkl = rel_pos_pbc(atoms, k, l)

716 rijrkj = np.dot(rij, rkj)

717 rkjrkl = np.dot(rkj, rkl)

719 rmj = np.cross(rij, rkj)

720 dmj = linalg.norm(rmj)

721 dmj2 = dmj*dmj

722 emj = rmj/dmj

723 rnk = np.cross(rkj, rkl)

724 dnk = linalg.norm(rnk)

725 dnk2 = dnk*dnk

726 enk = rnk/dnk

727 emjenk = np.dot(emj, enk)

728 if np.abs(emjenk) > 1.0:

729 emjenk = np.sign(emjenk)

731 d = np.sign(np.dot(rkj, np.cross(rmj, rnk)))*np.arccos(emjenk)

733 dddri = dkj/dmj2*rmj

734 dddrl = -dkj/dnk2*rnk

736 gx = np.zeros(12)

738 gx[0:3] = dddri

739 gx[3:6] = (rijrkj/dkj2-1.0)*dddri-rkjrkl/dkj2*dddrl

740 gx[6:9] = (rkjrkl/dkj2-1.0)*dddrl-rijrkj/dkj2*dddri

741 gx[9:12] = dddrl

743 if dihedral.d0 is None:

744 Hx = np.abs(2.0*dihedral.k*np.cos(2.0 * d))*np.tensordot(gx,gx,axes=0)

745 if dihedral.n is None:

746 Hx = dihedral.k*np.tensordot(gx,gx,axes=0)

747 else:

748 Hx = ( np.abs(-dihedral.k*dihedral.n**2

749 *np.cos(dihedral.n*d-dihedral.d0))*np.tensordot(gx,gx,axes=0) )

751 if dihedral.alpha is not None:

752 rij = rel_pos_pbc(atoms, i, j)

753 dij = linalg.norm(rij)

754 rkj = rel_pos_pbc(atoms, k, j)

755 dkj = linalg.norm(rkj)

756 rkl = rel_pos_pbc(atoms, k, l)

757 dkl = linalg.norm(rkl)

758 Hx *= ( np.exp(dihedral.alpha[0]*(dihedral.rref[0]**2-dij**2))

759 *np.exp(dihedral.alpha[1]*(dihedral.rref[1]**2-dkj**2))

760 *np.exp(dihedral.alpha[2]*(dihedral.rref[2]**2-dkl**2)) )

762 if morses is not None:

763 for m in range(len(morses)):

764 if ( morses[m].atomi == i or

765 morses[m].atomi == j or

766 morses[m].atomi == k or

767 morses[m].atomi == l ):

768 Hx *= get_morse_potential_eta(atoms, morses[m])

769 elif ( morses[m].atomj == i or

770 morses[m].atomj == j or

771 morses[m].atomj == k or

772 morses[m].atomj == l ):

773 Hx *= get_morse_potential_eta(atoms, morses[m])

775 dihedral.d = d

777 return i, j, k, l, Hx

779def get_dihedral_potential_reduced_hessian_test(atoms, dihedral):

780 i, j, k, l, gx = get_dihedral_potential_gradient(atoms, dihedral)

782 if dihedral.n is None:

783 i, j, k, l, v = get_dihedral_potential_value(atoms, dihedral)

784 Hx = np.tensordot(gx,gx,axes=0)/v/2.0

785 else:

786 arg = dihedral.n*dihedral.d - dihedral.d0

787 Hx = ( np.tensordot(gx,gx,axes=0)/dihedral.k/np.sin(arg)/np.sin(arg)

788 *np.cos(arg) )

790 return i, j, k, l, Hx

792def get_vdw_potential_value(atoms, vdw):

793 i = vdw.atomi

794 j = vdw.atomj

796 rij = rel_pos_pbc(atoms, i, j)

797 dij = linalg.norm(rij)

799 v = vdw.Aij/dij**12 - vdw.Bij/dij**6

801 vdw.r = dij

803 return i, j, v

806 i = vdw.atomi

807 j = vdw.atomj

809 rij = rel_pos_pbc(atoms, i, j)

810 dij = linalg.norm(rij)

811 eij = rij/dij

813 gr = (-12.0*vdw.Aij/dij**13+6.0*vdw.Bij/dij**7)*eij

815 gx = np.dot(Bx.T, gr)

817 vdw.r = dij

819 return i, j, gx

821def get_vdw_potential_hessian(atoms, vdw, spectral=False):

822 i = vdw.atomi

823 j = vdw.atomj

825 rij = rel_pos_pbc(atoms, i, j)

826 dij = linalg.norm(rij)

827 eij = rij/dij

829 Pij = np.tensordot(eij,eij,axes=0)

830 Qij = np.eye(3)-Pij

832 Hr = ( (156.0*vdw.Aij/dij**14-42.0*vdw.Bij/dij**8)*Pij

833 +(-12.0*vdw.Aij/dij**13+6.0*vdw.Bij/dij**7)/dij*Qij )

835 Hx = np.dot(Bx.T, np.dot(Hr, Bx))

837 if spectral:

838 eigvals, eigvecs = linalg.eigh(Hx)

839 D = np.diag(np.abs(eigvals))

840 U = eigvecs

841 Hx = np.dot(U,np.dot(D,np.transpose(U)))

843 vdw.r = dij

845 return i, j, Hx

847def get_coulomb_potential_value(atoms, coulomb):

848 i = coulomb.atomi

849 j = coulomb.atomj

851 rij = rel_pos_pbc(atoms, i, j)

852 dij = linalg.norm(rij)

854 v = coulomb.chargeij/dij

856 coulomb.r = dij

858 return i, j, v

861 i = coulomb.atomi

862 j = coulomb.atomj

864 rij = rel_pos_pbc(atoms, i, j)

865 dij = linalg.norm(rij)

866 eij = rij/dij

868 gr = -coulomb.chargeij/dij/dij*eij

870 gx = np.dot(Bx.T, gr)

872 coulomb.r = dij

874 return i, j, gx

876def get_coulomb_potential_hessian(atoms, coulomb, spectral=False):

877 i = coulomb.atomi

878 j = coulomb.atomj

880 rij = rel_pos_pbc(atoms, i, j)

881 dij = linalg.norm(rij)

882 eij = rij/dij

884 Pij = np.tensordot(eij,eij,axes=0)

885 Qij = np.eye(3)-Pij

887 Hr = (2.0*coulomb.chargeij/dij**3)*Pij+(-coulomb.chargeij/dij/dij)/dij*Qij

889 Hx = np.dot(Bx.T, np.dot(Hr, Bx))

891 if spectral:

892 eigvals, eigvecs = linalg.eigh(Hx)

893 D = np.diag(np.abs(eigvals))

894 U = eigvecs

895 Hx = np.dot(U,np.dot(D,np.transpose(U)))

897 coulomb.r = dij

899 return i, j, Hx

901def rel_pos_pbc(atoms, i, j):

902 """

903 Return difference between two atomic positions,

904 correcting for jumps across PBC

905 """

906 d = atoms.get_positions()[i,:]-atoms.get_positions()[j,:]

907 g = linalg.inv(atoms.get_cell().T)

908 f = np.floor(np.dot(g, d.T) + 0.5)

909 d -= np.dot(atoms.get_cell().T, f).T

910 return d