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

1# flake8: noqa

2import numpy as np

3from numpy import linalg

5from ase import units

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

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

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

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

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

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

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

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

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

16 [0, 0, 1, 0, 0, -1]])

17Mx = 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

30class Bond:

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

32 alpha=None, rref=None):

33 self.atomi = atomi

34 self.atomj = atomj

35 self.k = k

36 self.b0 = b0

37 self.alpha = alpha

38 self.rref = rref

39 self.b = None

42class Angle:

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

44 alpha=None, rref=None):

45 self.atomi = atomi

46 self.atomj = atomj

47 self.atomk = atomk

48 self.k = k

49 self.a0 = a0

50 self.cos = cos

51 self.alpha = alpha

52 self.rref = rref

53 self.a = None

56class Dihedral:

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

58 alpha=None, rref=None):

59 self.atomi = atomi

60 self.atomj = atomj

61 self.atomk = atomk

62 self.atoml = atoml

63 self.k = k

64 self.d0 = d0

65 self.n = n

66 self.alpha = alpha

67 self.rref = rref

68 self.d = None

71class VdW:

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

73 Aij=None, Bij=None, epsiloni=None, epsilonj=None,

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

75 self.atomi = atomi

76 self.atomj = atomj

77 if epsilonij is not None:

78 if sigmaij is not None:

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

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

81 elif rminij is not None:

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

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

84 else:

85 raise NotImplementedError("not implemented combination"

86 "of vdW parameters.")

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

88 self.Aij = scale * Aij

89 self.Bij = scale * Bij

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

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

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

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

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

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

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

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

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

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

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

101 else:

102 raise NotImplementedError("not implemented combination"

103 "of vdW parameters.")

104 self.r = None

105

106

107class Coulomb:

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

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

110 self.atomi = atomi

111 self.atomj = atomj

112 if chargeij is not None:

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

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

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

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

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

118 else:

119 raise NotImplementedError("not implemented combination"

120 "of Coulomb parameters.")

121 self.r = None

122

123

124def get_morse_potential_eta(atoms, morse):

125 i = morse.atomi

126 j = morse.atomj

127

128 rij = rel_pos_pbc(atoms, i, j)

129 dij = linalg.norm(rij)

130

131 if dij > morse.r0:

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

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

134 else:

135 eta = 1.0

136

137 return eta

138

139

140def get_morse_potential_value(atoms, morse):

141 i = morse.atomi

142 j = morse.atomj

143

144 rij = rel_pos_pbc(atoms, i, j)

145 dij = linalg.norm(rij)

146

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

148

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

151 morse.r = dij

153 return i, j, v

156def get_morse_potential_gradient(atoms, morse):

157 i = morse.atomi

158 j = morse.atomj

159

160 rij = rel_pos_pbc(atoms, i, j)

161 dij = linalg.norm(rij)

162 eij = rij / dij

163

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

165

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

167

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

170 morse.r = dij

172 return i, j, gx

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

176 i = morse.atomi

177 j = morse.atomj

178

179 rij = rel_pos_pbc(atoms, i, j)

180 dij = linalg.norm(rij)

181 eij = rij / dij

182

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

184 Qij = np.eye(3) - Pij

185

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

187

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

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

190

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

192

193 if spectral:

194 eigvals, eigvecs = linalg.eigh(Hx)

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

196 U = eigvecs

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

199 morse.r = dij

201 return i, j, Hx

204def get_morse_potential_reduced_hessian(atoms, morse):

205 i = morse.atomi

206 j = morse.atomj

207

208 rij = rel_pos_pbc(atoms, i, j)

209 dij = linalg.norm(rij)

210 eij = rij / dij

211

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

213

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

215

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

217

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

220 morse.r = dij

222 return i, j, Hx

225def get_bond_potential_value(atoms, bond):

226 i = bond.atomi

227 j = bond.atomj

228

229 rij = rel_pos_pbc(atoms, i, j)

230 dij = linalg.norm(rij)

231

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

234 bond.b = dij

236 return i, j, v

239def get_bond_potential_gradient(atoms, bond):

240 i = bond.atomi

241 j = bond.atomj

242

243 rij = rel_pos_pbc(atoms, i, j)

244 dij = linalg.norm(rij)

245 eij = rij / dij

246

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

248

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

251 bond.b = dij

253 return i, j, gx

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

257 i = bond.atomi

258 j = bond.atomj

259

260 rij = rel_pos_pbc(atoms, i, j)

261 dij = linalg.norm(rij)

262 eij = rij / dij

263

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

265 Qij = np.eye(3) - Pij

266

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

268

269 if bond.alpha is not None:

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

271

272 if morses is not None:

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

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

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

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

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

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

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

280

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

282

283 if spectral:

284 eigvals, eigvecs = linalg.eigh(Hx)

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

286 U = eigvecs

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

289 bond.b = dij

291 return i, j, Hx

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

295 i = bond.atomi

296 j = bond.atomj

297

298 rij = rel_pos_pbc(atoms, i, j)

299 dij = linalg.norm(rij)

300 eij = rij / dij

301

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

303

304 Hr = bond.k * Pij

305

306 if bond.alpha is not None:

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

308

309 if morses is not None:

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

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

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

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

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

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

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

317

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

320 bond.b = dij

322 return i, j, Hx

325def get_bond_potential_reduced_hessian_test(atoms, bond):

326

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

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

329

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

331

332 return i, j, Hx

333

334

335def get_angle_potential_value(atoms, angle):

336

337 i = angle.atomi

338 j = angle.atomj

339 k = angle.atomk

340

341 rij = rel_pos_pbc(atoms, i, j)

342 dij = linalg.norm(rij)

343 eij = rij / dij

344 rkj = rel_pos_pbc(atoms, k, j)

345 dkj = linalg.norm(rkj)

346 ekj = rkj / dkj

347 eijekj = np.dot(eij, ekj)

348 if np.abs(eijekj) > 1.0:

349 eijekj = np.sign(eijekj)

350

351 a = np.arccos(eijekj)

352

353 if angle.cos:

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

355 else:

356 da = a - angle.a0

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

358

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

360

361 angle.a = a

362

363 return i, j, k, v

364

365

366def get_angle_potential_gradient(atoms, angle):

367 i = angle.atomi

368 j = angle.atomj

369 k = angle.atomk

370

371 rij = rel_pos_pbc(atoms, i, j)

372 dij = linalg.norm(rij)

373 eij = rij / dij

374 rkj = rel_pos_pbc(atoms, k, j)

375 dkj = linalg.norm(rkj)

376 ekj = rkj / dkj

377 eijekj = np.dot(eij, ekj)

378 if np.abs(eijekj) > 1.0:

379 eijekj = np.sign(eijekj)

380

381 a = np.arccos(eijekj)

382 if angle.cos:

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

384 else:

385 da = a - angle.a0

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

387 sina = np.sin(a)

388

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

390 Qij = np.eye(3) - Pij

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

392 Qkj = np.eye(3) - Pkj

393

394 gr = np.zeros(6)

395 if angle.cos:

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

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

398 elif np.abs(sina) > 0.001:

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

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

401

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

403

404 angle.a = a

405

406 return i, j, k, gx

407

408

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

410 i = angle.atomi

411 j = angle.atomj

412 k = angle.atomk

413

414 rij = rel_pos_pbc(atoms, i, j)

415 dij = linalg.norm(rij)

416 dij2 = dij * dij

417 eij = rij / dij

418 rkj = rel_pos_pbc(atoms, k, j)

419 dkj = linalg.norm(rkj)

420 dkj2 = dkj * dkj

421 ekj = rkj / dkj

422 dijdkj = dij * dkj

423 eijekj = np.dot(eij, ekj)

424 if np.abs(eijekj) > 1.0:

425 eijekj = np.sign(eijekj)

426

427 a = np.arccos(eijekj)

428 if angle.cos:

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

430 cosa0 = np.cos(angle.a0)

431 else:

432 da = a - angle.a0

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

434 sina = np.sin(a)

435 cosa = np.cos(a)

436 ctga = cosa / sina

437

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

439 Qij = np.eye(3) - Pij

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

441 Qkj = np.eye(3) - Pkj

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

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

444 P = np.eye(3) * eijekj

445

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

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

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

449

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

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

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

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

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

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

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

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

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

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

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

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

462 + np.dot(Qkj, Pik) -

463 np.dot(Pkj, Pik)

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

465 elif np.abs(sina) > 0.001:

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

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

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

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

470 + da * (-ctga * QijPkiQkj / sina

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

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

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

474 + da * (-ctga * QkjPijQkj / sina

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

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

477

478 if angle.alpha is not None:

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

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

481

482 if morses is not None:

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

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

485 morses[m].atomi == j or

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

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

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

489 morses[m].atomj == j or

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

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

492

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

494

495 if spectral:

496 eigvals, eigvecs = linalg.eigh(Hx)

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

498 U = eigvecs

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

500

501 angle.a = a

502

503 return i, j, k, Hx

504

505

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

507 i = angle.atomi

508 j = angle.atomj

509 k = angle.atomk

510

511 rij = rel_pos_pbc(atoms, i, j)

512 dij = linalg.norm(rij)

513 dij2 = dij * dij

514 eij = rij / dij

515 rkj = rel_pos_pbc(atoms, k, j)

516 dkj = linalg.norm(rkj)

517 dkj2 = dkj * dkj

518 ekj = rkj / dkj

519 dijdkj = dij * dkj

520 eijekj = np.dot(eij, ekj)

521 if np.abs(eijekj) > 1.0:

522 eijekj = np.sign(eijekj)

523

524 a = np.arccos(eijekj)

525 sina = np.sin(a)

526 sina2 = sina * sina

527

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

529 Qij = np.eye(3) - Pij

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

531 Qkj = np.eye(3) - Pkj

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

533

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

535 if np.abs(sina) > 0.001:

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

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

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

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

540

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

542 cosa = np.cos(a)

543 cosa0 = np.cos(angle.a0)

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

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

546 elif np.abs(sina) > 0.001:

547 Hr = Hr * angle.k / sina2

548

549 if angle.alpha is not None:

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

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

552

553 if morses is not None:

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

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

556 morses[m].atomi == j or

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

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

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

560 morses[m].atomj == j or

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

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

563

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

565

566 angle.a = a

567

568 return i, j, k, Hx

569

570

571def get_angle_potential_reduced_hessian_test(atoms, angle):

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

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

574

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

576

577 return i, j, k, Hx

578

579

580def get_dihedral_potential_value(atoms, dihedral):

581 i = dihedral.atomi

582 j = dihedral.atomj

583 k = dihedral.atomk

584 l = dihedral.atoml

585

586 rij = rel_pos_pbc(atoms, i, j)

587 rkj = rel_pos_pbc(atoms, k, j)

588 rkl = rel_pos_pbc(atoms, k, l)

589

590 rmj = np.cross(rij, rkj)

591 dmj = linalg.norm(rmj)

592 emj = rmj / dmj

593 rnk = np.cross(rkj, rkl)

594 dnk = linalg.norm(rnk)

595 enk = rnk / dnk

596 emjenk = np.dot(emj, enk)

597 if np.abs(emjenk) > 1.0:

598 emjenk = np.sign(emjenk)

599

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

601

602 if dihedral.d0 is None:

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

604 else:

605 dd = d - dihedral.d0

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

607 if dihedral.n is None:

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

609 else:

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

611

612 dihedral.d = d

613

614 return i, j, k, l, v

615

616

617def get_dihedral_potential_gradient(atoms, dihedral):

618 i = dihedral.atomi

619 j = dihedral.atomj

620 k = dihedral.atomk

621 l = dihedral.atoml

622

623 rij = rel_pos_pbc(atoms, i, j)

624 rkj = rel_pos_pbc(atoms, k, j)

625 dkj = linalg.norm(rkj)

626 dkj2 = dkj * dkj

627 rkl = rel_pos_pbc(atoms, k, l)

628

629 rijrkj = np.dot(rij, rkj)

630 rkjrkl = np.dot(rkj, rkl)

631

632 rmj = np.cross(rij, rkj)

633 dmj = linalg.norm(rmj)

634 dmj2 = dmj * dmj

635 emj = rmj / dmj

636 rnk = np.cross(rkj, rkl)

637 dnk = linalg.norm(rnk)

638 dnk2 = dnk * dnk

639 enk = rnk / dnk

640 emjenk = np.dot(emj, enk)

641 if np.abs(emjenk) > 1.0:

642 emjenk = np.sign(emjenk)

643

644 dddri = dkj / dmj2 * rmj

645 dddrl = -dkj / dnk2 * rnk

646

647 gx = np.zeros(12)

648

649 gx[0:3] = dddri

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

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

652 gx[9:12] = dddrl

653

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

655

656 if dihedral.d0 is None:

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

658 else:

659 dd = d - dihedral.d0

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

661 if dihedral.n is None:

662 gx *= dihedral.k * dd

663 else:

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

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

666

667 dihedral.d = d

668

669 return i, j, k, l, gx

670

671

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

673 spectral=False):

674 eps = 0.000001

675

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

677

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

679

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

681 dihedral.atomk, dihedral.atoml,

682 dihedral.k, dihedral.d0, dihedral.n)

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

684 3 * j, 3 * j + 1, 3 * j + 2,

685 3 * k, 3 * k + 1, 3 * k + 2,

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

687 for x in range(12):

688 a = atoms.copy()

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

690 positions[indx[x]] += eps

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

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

693 for y in range(12):

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

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

696

697 if dihedral.alpha is not None:

698 rij = rel_pos_pbc(atoms, i, j)

699 dij = linalg.norm(rij)

700 rkj = rel_pos_pbc(atoms, k, j)

701 dkj = linalg.norm(rkj)

702 rkl = rel_pos_pbc(atoms, k, l)

703 dkl = linalg.norm(rkl)

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

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

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

707

708 if morses is not None:

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

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

711 morses[m].atomi == j or

712 morses[m].atomi == k or

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

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

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

716 morses[m].atomj == j or

717 morses[m].atomj == k or

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

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

720

721 if spectral:

722 eigvals, eigvecs = linalg.eigh(Hx)

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

724 U = eigvecs

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

726

727 return i, j, k, l, Hx

728

729

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

731 i = dihedral.atomi

732 j = dihedral.atomj

733 k = dihedral.atomk

734 l = dihedral.atoml

735

736 rij = rel_pos_pbc(atoms, i, j)

737 rkj = rel_pos_pbc(atoms, k, j)

738 dkj = linalg.norm(rkj)

739 dkj2 = dkj * dkj

740 rkl = rel_pos_pbc(atoms, k, l)

741

742 rijrkj = np.dot(rij, rkj)

743 rkjrkl = np.dot(rkj, rkl)

744

745 rmj = np.cross(rij, rkj)

746 dmj = linalg.norm(rmj)

747 dmj2 = dmj * dmj

748 emj = rmj / dmj

749 rnk = np.cross(rkj, rkl)

750 dnk = linalg.norm(rnk)

751 dnk2 = dnk * dnk

752 enk = rnk / dnk

753 emjenk = np.dot(emj, enk)

754 if np.abs(emjenk) > 1.0:

755 emjenk = np.sign(emjenk)

756

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

758

759 dddri = dkj / dmj2 * rmj

760 dddrl = -dkj / dnk2 * rnk

761

762 gx = np.zeros(12)

763

764 gx[0:3] = dddri

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

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

767 gx[9:12] = dddrl

768

769 if dihedral.d0 is None:

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

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

772 if dihedral.n is None:

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

774 else:

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

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

777

778 if dihedral.alpha is not None:

779 rij = rel_pos_pbc(atoms, i, j)

780 dij = linalg.norm(rij)

781 rkj = rel_pos_pbc(atoms, k, j)

782 dkj = linalg.norm(rkj)

783 rkl = rel_pos_pbc(atoms, k, l)

784 dkl = linalg.norm(rkl)

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

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

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

788

789 if morses is not None:

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

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

792 morses[m].atomi == j or

793 morses[m].atomi == k or

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

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

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

797 morses[m].atomj == j or

798 morses[m].atomj == k or

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

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

801

802 dihedral.d = d

803

804 return i, j, k, l, Hx

805

806

807def get_dihedral_potential_reduced_hessian_test(atoms, dihedral):

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

809

810 if dihedral.n is None:

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

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

813 else:

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

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

816 * np.cos(arg))

817

818 return i, j, k, l, Hx

819

820

821def get_vdw_potential_value(atoms, vdw):

822 i = vdw.atomi

823 j = vdw.atomj

824

825 rij = rel_pos_pbc(atoms, i, j)

826 dij = linalg.norm(rij)

827

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

830 vdw.r = dij

832 return i, j, v

835def get_vdw_potential_gradient(atoms, vdw):

836 i = vdw.atomi

837 j = vdw.atomj

838

839 rij = rel_pos_pbc(atoms, i, j)

840 dij = linalg.norm(rij)

841 eij = rij / dij

842

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

844

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

847 vdw.r = dij

849 return i, j, gx

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

853 i = vdw.atomi

854 j = vdw.atomj

855

856 rij = rel_pos_pbc(atoms, i, j)

857 dij = linalg.norm(rij)

858 eij = rij / dij

859

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

861 Qij = np.eye(3) - Pij

862

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

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

865

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

867

868 if spectral:

869 eigvals, eigvecs = linalg.eigh(Hx)

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

871 U = eigvecs

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

874 vdw.r = dij

876 return i, j, Hx

879def get_coulomb_potential_value(atoms, coulomb):

880 i = coulomb.atomi

881 j = coulomb.atomj

882

883 rij = rel_pos_pbc(atoms, i, j)

884 dij = linalg.norm(rij)

885

886 v = coulomb.chargeij / dij

888 coulomb.r = dij

890 return i, j, v

893def get_coulomb_potential_gradient(atoms, coulomb):

894 i = coulomb.atomi

895 j = coulomb.atomj

896

897 rij = rel_pos_pbc(atoms, i, j)

898 dij = linalg.norm(rij)

899 eij = rij / dij

900

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

902

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

905 coulomb.r = dij

907 return i, j, gx

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

911 i = coulomb.atomi

912 j = coulomb.atomj

913

914 rij = rel_pos_pbc(atoms, i, j)

915 dij = linalg.norm(rij)

916 eij = rij / dij

917

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

919 Qij = np.eye(3) - Pij

920

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

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

923

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

925

926 if spectral:

927 eigvals, eigvecs = linalg.eigh(Hx)

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

929 U = eigvecs

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

932 coulomb.r = dij

934 return i, j, Hx

937def rel_pos_pbc(atoms, i, j):

938 """

939 Return difference between two atomic positions,

940 correcting for jumps across PBC

941 """

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

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

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

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

946 return d