Coverage for /builds/debichem-team/python-ase/ase/utils/ff.py : 57.72%
Hot-keys on this page
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
148def get_morse_potential_gradient(atoms, morse):
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
227def get_bond_potential_gradient(atoms, bond):
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
349def get_angle_potential_gradient(atoms, angle):
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
594def get_dihedral_potential_gradient(atoms, dihedral):
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
651 i,j,k,l,g = get_dihedral_potential_gradient(atoms, dihedral)
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)))
667 i,j,k,l,geps = get_dihedral_potential_gradient(a, dihedral_eps)
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
805def get_vdw_potential_gradient(atoms, vdw):
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
860def get_coulomb_potential_gradient(atoms, coulomb):
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