Coverage for /builds/ase/ase/ase/utils/ff.py: 58.86%

1# fmt: off

3# flake8: noqa

4import numpy as np

5from numpy import linalg

7from ase import units

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

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

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

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

13 [0, 0, 0, -1, 0, 0, 1, 0, 0],

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

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

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

17 [0, 1, 0, 0, -1, 0],

18 [0, 0, 1, 0, 0, -1]])

19Mx = Bx

22class Morse:

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

24 self.atomi = atomi

25 self.atomj = atomj

26 self.D = D

27 self.alpha = alpha

28 self.r0 = r0

29 self.r = None

32class Bond:

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

34 alpha=None, rref=None):

35 self.atomi = atomi

36 self.atomj = atomj

37 self.k = k

38 self.b0 = b0

39 self.alpha = alpha

40 self.rref = rref

41 self.b = None

44class Angle:

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

46 alpha=None, rref=None):

47 self.atomi = atomi

48 self.atomj = atomj

49 self.atomk = atomk

50 self.k = k

51 self.a0 = a0

52 self.cos = cos

53 self.alpha = alpha

54 self.rref = rref

55 self.a = None

58class Dihedral:

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

60 alpha=None, rref=None):

61 self.atomi = atomi

62 self.atomj = atomj

63 self.atomk = atomk

64 self.atoml = atoml

65 self.k = k

66 self.d0 = d0

67 self.n = n

68 self.alpha = alpha

69 self.rref = rref

70 self.d = None

73class VdW:

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

75 Aij=None, Bij=None, epsiloni=None, epsilonj=None,

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

77 self.atomi = atomi

78 self.atomj = atomj

79 if epsilonij is not None:

80 if sigmaij is not None:

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

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

83 elif rminij is not None:

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

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

86 else:

87 raise NotImplementedError("not implemented combination"

88 "of vdW parameters.")

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

90 self.Aij = scale * Aij

91 self.Bij = scale * Bij

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

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

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

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

96 self.Bij = (scale * 4.0 * np.sqrt(epsiloni * epsilonj)

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

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

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

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

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

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

103 else:

104 raise NotImplementedError("not implemented combination"

105 "of vdW parameters.")

106 self.r = None

107

108

109class Coulomb:

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

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

112 self.atomi = atomi

113 self.atomj = atomj

114 if chargeij is not None:

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

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

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

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

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

120 else:

121 raise NotImplementedError("not implemented combination"

122 "of Coulomb parameters.")

123 self.r = None

124

125

126def get_morse_potential_eta(atoms, morse):

127 i = morse.atomi

128 j = morse.atomj

129

130 rij = rel_pos_pbc(atoms, i, j)

131 dij = linalg.norm(rij)

132

133 if dij > morse.r0:

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

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

136 else:

137 eta = 1.0

138

139 return eta

140

141

142def get_morse_potential_value(atoms, morse):

143 i = morse.atomi

144 j = morse.atomj

145

146 rij = rel_pos_pbc(atoms, i, j)

147 dij = linalg.norm(rij)

148

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

150

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

153 morse.r = dij

155 return i, j, v

158def get_morse_potential_gradient(atoms, morse):

159 i = morse.atomi

160 j = morse.atomj

161

162 rij = rel_pos_pbc(atoms, i, j)

163 dij = linalg.norm(rij)

164 eij = rij / dij

165

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

167

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

169

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

172 morse.r = dij

174 return i, j, gx

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

178 i = morse.atomi

179 j = morse.atomj

180

181 rij = rel_pos_pbc(atoms, i, j)

182 dij = linalg.norm(rij)

183 eij = rij / dij

184

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

186 Qij = np.eye(3) - Pij

187

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

189

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

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

192

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

194

195 if spectral:

196 eigvals, eigvecs = linalg.eigh(Hx)

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

198 U = eigvecs

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

201 morse.r = dij

203 return i, j, Hx

206def get_morse_potential_reduced_hessian(atoms, morse):

207 i = morse.atomi

208 j = morse.atomj

209

210 rij = rel_pos_pbc(atoms, i, j)

211 dij = linalg.norm(rij)

212 eij = rij / dij

213

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

215

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

217

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

219

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

222 morse.r = dij

224 return i, j, Hx

227def get_bond_potential_value(atoms, bond):

228 i = bond.atomi

229 j = bond.atomj

230

231 rij = rel_pos_pbc(atoms, i, j)

232 dij = linalg.norm(rij)

233

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

236 bond.b = dij

238 return i, j, v

241def get_bond_potential_gradient(atoms, bond):

242 i = bond.atomi

243 j = bond.atomj

244

245 rij = rel_pos_pbc(atoms, i, j)

246 dij = linalg.norm(rij)

247 eij = rij / dij

248

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

250

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

253 bond.b = dij

255 return i, j, gx

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

259 i = bond.atomi

260 j = bond.atomj

261

262 rij = rel_pos_pbc(atoms, i, j)

263 dij = linalg.norm(rij)

264 eij = rij / dij

265

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

267 Qij = np.eye(3) - Pij

268

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

270

271 if bond.alpha is not None:

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

273

274 if morses is not None:

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

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

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

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

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

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

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

282

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

284

285 if spectral:

286 eigvals, eigvecs = linalg.eigh(Hx)

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

288 U = eigvecs

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

291 bond.b = dij

293 return i, j, Hx

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

297 i = bond.atomi

298 j = bond.atomj

299

300 rij = rel_pos_pbc(atoms, i, j)

301 dij = linalg.norm(rij)

302 eij = rij / dij

303

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

305

306 Hr = bond.k * Pij

307

308 if bond.alpha is not None:

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

310

311 if morses is not None:

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

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

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

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

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

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

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

319

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

322 bond.b = dij

324 return i, j, Hx

327def get_bond_potential_reduced_hessian_test(atoms, bond):

328

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

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

331

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

333

334 return i, j, Hx

335

336

337def get_angle_potential_value(atoms, angle):

338

339 i = angle.atomi

340 j = angle.atomj

341 k = angle.atomk

342

343 rij = rel_pos_pbc(atoms, i, j)

344 dij = linalg.norm(rij)

345 eij = rij / dij

346 rkj = rel_pos_pbc(atoms, k, j)

347 dkj = linalg.norm(rkj)

348 ekj = rkj / dkj

349 eijekj = np.dot(eij, ekj)

350 if np.abs(eijekj) > 1.0:

351 eijekj = np.sign(eijekj)

352

353 a = np.arccos(eijekj)

354

355 if angle.cos:

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

357 else:

358 da = a - angle.a0

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

360

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

362

363 angle.a = a

364

365 return i, j, k, v

366

367

368def get_angle_potential_gradient(atoms, angle):

369 i = angle.atomi

370 j = angle.atomj

371 k = angle.atomk

372

373 rij = rel_pos_pbc(atoms, i, j)

374 dij = linalg.norm(rij)

375 eij = rij / dij

376 rkj = rel_pos_pbc(atoms, k, j)

377 dkj = linalg.norm(rkj)

378 ekj = rkj / dkj

379 eijekj = np.dot(eij, ekj)

380 if np.abs(eijekj) > 1.0:

381 eijekj = np.sign(eijekj)

382

383 a = np.arccos(eijekj)

384 if angle.cos:

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

386 else:

387 da = a - angle.a0

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

389 sina = np.sin(a)

390

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

392 Qij = np.eye(3) - Pij

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

394 Qkj = np.eye(3) - Pkj

395

396 gr = np.zeros(6)

397 if angle.cos:

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

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

400 elif np.abs(sina) > 0.001:

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

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

403

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

405

406 angle.a = a

407

408 return i, j, k, gx

409

410

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

412 i = angle.atomi

413 j = angle.atomj

414 k = angle.atomk

415

416 rij = rel_pos_pbc(atoms, i, j)

417 dij = linalg.norm(rij)

418 dij2 = dij * dij

419 eij = rij / dij

420 rkj = rel_pos_pbc(atoms, k, j)

421 dkj = linalg.norm(rkj)

422 dkj2 = dkj * dkj

423 ekj = rkj / dkj

424 dijdkj = dij * dkj

425 eijekj = np.dot(eij, ekj)

426 if np.abs(eijekj) > 1.0:

427 eijekj = np.sign(eijekj)

428

429 a = np.arccos(eijekj)

430 if angle.cos:

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

432 cosa0 = np.cos(angle.a0)

433 else:

434 da = a - angle.a0

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

436 sina = np.sin(a)

437 cosa = np.cos(a)

438 ctga = cosa / sina

439

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

441 Qij = np.eye(3) - Pij

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

443 Qkj = np.eye(3) - Pkj

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

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

446 P = np.eye(3) * eijekj

447

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

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

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

451

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

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

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

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

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

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

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

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

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

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

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

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

464 + np.dot(Qkj, Pik) -

465 np.dot(Pkj, Pik)

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

467 elif np.abs(sina) > 0.001:

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

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

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

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

472 + da * (-ctga * QijPkiQkj / sina

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

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

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

476 + da * (-ctga * QkjPijQkj / sina

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

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

479

480 if angle.alpha is not None:

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

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

483

484 if morses is not None:

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

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

487 morses[m].atomi == j or

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

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

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

491 morses[m].atomj == j or

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

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

494

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

496

497 if spectral:

498 eigvals, eigvecs = linalg.eigh(Hx)

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

500 U = eigvecs

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

502

503 angle.a = a

504

505 return i, j, k, Hx

506

507

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

509 i = angle.atomi

510 j = angle.atomj

511 k = angle.atomk

512

513 rij = rel_pos_pbc(atoms, i, j)

514 dij = linalg.norm(rij)

515 dij2 = dij * dij

516 eij = rij / dij

517 rkj = rel_pos_pbc(atoms, k, j)

518 dkj = linalg.norm(rkj)

519 dkj2 = dkj * dkj

520 ekj = rkj / dkj

521 dijdkj = dij * dkj

522 eijekj = np.dot(eij, ekj)

523 if np.abs(eijekj) > 1.0:

524 eijekj = np.sign(eijekj)

525

526 a = np.arccos(eijekj)

527 sina = np.sin(a)

528 sina2 = sina * sina

529

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

531 Qij = np.eye(3) - Pij

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

533 Qkj = np.eye(3) - Pkj

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

535

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

537 if np.abs(sina) > 0.001:

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

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

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

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

542

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

544 cosa = np.cos(a)

545 cosa0 = np.cos(angle.a0)

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

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

548 elif np.abs(sina) > 0.001:

549 Hr = Hr * angle.k / sina2

550

551 if angle.alpha is not None:

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

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

554

555 if morses is not None:

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

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

558 morses[m].atomi == j or

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

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

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

562 morses[m].atomj == j or

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

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

565

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

567

568 angle.a = a

569

570 return i, j, k, Hx

571

572

573def get_angle_potential_reduced_hessian_test(atoms, angle):

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

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

576

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

578

579 return i, j, k, Hx

580

581

582def get_dihedral_potential_value(atoms, dihedral):

583 i = dihedral.atomi

584 j = dihedral.atomj

585 k = dihedral.atomk

586 l = dihedral.atoml

587

588 rij = rel_pos_pbc(atoms, i, j)

589 rkj = rel_pos_pbc(atoms, k, j)

590 rkl = rel_pos_pbc(atoms, k, l)

591

592 rmj = np.cross(rij, rkj)

593 dmj = linalg.norm(rmj)

594 emj = rmj / dmj

595 rnk = np.cross(rkj, rkl)

596 dnk = linalg.norm(rnk)

597 enk = rnk / dnk

598 emjenk = np.dot(emj, enk)

599 if np.abs(emjenk) > 1.0:

600 emjenk = np.sign(emjenk)

601

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

603

604 if dihedral.d0 is None:

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

606 else:

607 dd = d - dihedral.d0

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

609 if dihedral.n is None:

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

611 else:

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

613

614 dihedral.d = d

615

616 return i, j, k, l, v

617

618

619def get_dihedral_potential_gradient(atoms, dihedral):

620 i = dihedral.atomi

621 j = dihedral.atomj

622 k = dihedral.atomk

623 l = dihedral.atoml

624

625 rij = rel_pos_pbc(atoms, i, j)

626 rkj = rel_pos_pbc(atoms, k, j)

627 dkj = linalg.norm(rkj)

628 dkj2 = dkj * dkj

629 rkl = rel_pos_pbc(atoms, k, l)

630

631 rijrkj = np.dot(rij, rkj)

632 rkjrkl = np.dot(rkj, rkl)

633

634 rmj = np.cross(rij, rkj)

635 dmj = linalg.norm(rmj)

636 dmj2 = dmj * dmj

637 emj = rmj / dmj

638 rnk = np.cross(rkj, rkl)

639 dnk = linalg.norm(rnk)

640 dnk2 = dnk * dnk

641 enk = rnk / dnk

642 emjenk = np.dot(emj, enk)

643 if np.abs(emjenk) > 1.0:

644 emjenk = np.sign(emjenk)

645

646 dddri = dkj / dmj2 * rmj

647 dddrl = -dkj / dnk2 * rnk

648

649 gx = np.zeros(12)

650

651 gx[0:3] = dddri

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

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

654 gx[9:12] = dddrl

655

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

657

658 if dihedral.d0 is None:

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

660 else:

661 dd = d - dihedral.d0

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

663 if dihedral.n is None:

664 gx *= dihedral.k * dd

665 else:

666 gx *= -dihedral.k * dihedral.n * \

667 np.sin(dihedral.n * d - dihedral.d0)

668

669 dihedral.d = d

670

671 return i, j, k, l, gx

672

673

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

675 spectral=False):

676 eps = 0.000001

677

678 i, j, k, l, g = get_dihedral_potential_gradient(atoms, dihedral)

679

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

681

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

683 dihedral.atomk, dihedral.atoml,

684 dihedral.k, dihedral.d0, dihedral.n)

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

686 3 * j, 3 * j + 1, 3 * j + 2,

687 3 * k, 3 * k + 1, 3 * k + 2,

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

689 for x in range(12):

690 a = atoms.copy()

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

692 positions[indx[x]] += eps

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

694 i, j, k, l, geps = get_dihedral_potential_gradient(a, dihedral_eps)

695 for y in range(12):

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

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

698

699 if dihedral.alpha is not None:

700 rij = rel_pos_pbc(atoms, i, j)

701 dij = linalg.norm(rij)

702 rkj = rel_pos_pbc(atoms, k, j)

703 dkj = linalg.norm(rkj)

704 rkl = rel_pos_pbc(atoms, k, l)

705 dkl = linalg.norm(rkl)

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

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

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

709

710 if morses is not None:

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

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

713 morses[m].atomi == j or

714 morses[m].atomi == k or

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

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

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

718 morses[m].atomj == j or

719 morses[m].atomj == k or

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

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

722

723 if spectral:

724 eigvals, eigvecs = linalg.eigh(Hx)

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

726 U = eigvecs

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

728

729 return i, j, k, l, Hx

730

731

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

733 i = dihedral.atomi

734 j = dihedral.atomj

735 k = dihedral.atomk

736 l = dihedral.atoml

737

738 rij = rel_pos_pbc(atoms, i, j)

739 rkj = rel_pos_pbc(atoms, k, j)

740 dkj = linalg.norm(rkj)

741 dkj2 = dkj * dkj

742 rkl = rel_pos_pbc(atoms, k, l)

743

744 rijrkj = np.dot(rij, rkj)

745 rkjrkl = np.dot(rkj, rkl)

746

747 rmj = np.cross(rij, rkj)

748 dmj = linalg.norm(rmj)

749 dmj2 = dmj * dmj

750 emj = rmj / dmj

751 rnk = np.cross(rkj, rkl)

752 dnk = linalg.norm(rnk)

753 dnk2 = dnk * dnk

754 enk = rnk / dnk

755 emjenk = np.dot(emj, enk)

756 if np.abs(emjenk) > 1.0:

757 emjenk = np.sign(emjenk)

758

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

760

761 dddri = dkj / dmj2 * rmj

762 dddrl = -dkj / dnk2 * rnk

763

764 gx = np.zeros(12)

765

766 gx[0:3] = dddri

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

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

769 gx[9:12] = dddrl

770

771 if dihedral.d0 is None:

772 Hx = np.abs(2.0 * dihedral.k * np.cos(2.0 * d)) * \

773 np.tensordot(gx, gx, axes=0)

774 if dihedral.n is None:

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

776 else:

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

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

779

780 if dihedral.alpha is not None:

781 rij = rel_pos_pbc(atoms, i, j)

782 dij = linalg.norm(rij)

783 rkj = rel_pos_pbc(atoms, k, j)

784 dkj = linalg.norm(rkj)

785 rkl = rel_pos_pbc(atoms, k, l)

786 dkl = linalg.norm(rkl)

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

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

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

790

791 if morses is not None:

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

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

794 morses[m].atomi == j or

795 morses[m].atomi == k or

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

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

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

799 morses[m].atomj == j or

800 morses[m].atomj == k or

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

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

803

804 dihedral.d = d

805

806 return i, j, k, l, Hx

807

808

809def get_dihedral_potential_reduced_hessian_test(atoms, dihedral):

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

811

812 if dihedral.n is None:

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

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

815 else:

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

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

818 * np.cos(arg))

819

820 return i, j, k, l, Hx

821

822

823def get_vdw_potential_value(atoms, vdw):

824 i = vdw.atomi

825 j = vdw.atomj

826

827 rij = rel_pos_pbc(atoms, i, j)

828 dij = linalg.norm(rij)

829

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

832 vdw.r = dij

834 return i, j, v

837def get_vdw_potential_gradient(atoms, vdw):

838 i = vdw.atomi

839 j = vdw.atomj

840

841 rij = rel_pos_pbc(atoms, i, j)

842 dij = linalg.norm(rij)

843 eij = rij / dij

844

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

846

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

849 vdw.r = dij

851 return i, j, gx

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

855 i = vdw.atomi

856 j = vdw.atomj

857

858 rij = rel_pos_pbc(atoms, i, j)

859 dij = linalg.norm(rij)

860 eij = rij / dij

861

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

863 Qij = np.eye(3) - Pij

864

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

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

867

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

869

870 if spectral:

871 eigvals, eigvecs = linalg.eigh(Hx)

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

873 U = eigvecs

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

876 vdw.r = dij

878 return i, j, Hx

881def get_coulomb_potential_value(atoms, coulomb):

882 i = coulomb.atomi

883 j = coulomb.atomj

884

885 rij = rel_pos_pbc(atoms, i, j)

886 dij = linalg.norm(rij)

887

888 v = coulomb.chargeij / dij

890 coulomb.r = dij

892 return i, j, v

895def get_coulomb_potential_gradient(atoms, coulomb):

896 i = coulomb.atomi

897 j = coulomb.atomj

898

899 rij = rel_pos_pbc(atoms, i, j)

900 dij = linalg.norm(rij)

901 eij = rij / dij

902

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

904

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

907 coulomb.r = dij

909 return i, j, gx

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

913 i = coulomb.atomi

914 j = coulomb.atomj

915

916 rij = rel_pos_pbc(atoms, i, j)

917 dij = linalg.norm(rij)

918 eij = rij / dij

919

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

921 Qij = np.eye(3) - Pij

922

923 Hr = (2.0 * coulomb.chargeij / dij**3) * Pij + \

924 (-coulomb.chargeij / dij / dij) / dij * Qij

925

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

927

928 if spectral:

929 eigvals, eigvecs = linalg.eigh(Hx)

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

931 U = eigvecs

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

934 coulomb.r = dij

936 return i, j, Hx

939def rel_pos_pbc(atoms, i, j):

940 """

941 Return difference between two atomic positions,

942 correcting for jumps across PBC

943 """

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

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

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

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

948 return d