| 1 |
subroutine do_electrostatic_new() |
| 2 |
|
| 3 |
ci = ElectrostaticMap(me1)%charge |
| 4 |
|
| 5 |
mui = ElectrostaticMap(me1)%dipole_moment |
| 6 |
Thetai = ElectrostaticMap(me1)%quadrupole_moments |
| 7 |
eFramei = eFrame(atom1) |
| 8 |
uz_i = eFramei(:,3) |
| 9 |
|
| 10 |
|
| 11 |
di = mui * uz_i |
| 12 |
qi = Thetai(:)*eFramei |
| 13 |
|
| 14 |
|
| 15 |
f = electric / dielec |
| 16 |
xji = x(jj) - x(ii) |
| 17 |
yji = y(jj) - y(ii) |
| 18 |
zji = z(jj) - z(ii) |
| 19 |
r2 = xji*xji + yji*yji + zji*zji |
| 20 |
r = sqrt(r2) |
| 21 |
rr1 = 1.0d0 / r |
| 22 |
rr3 = rr1 / r2 |
| 23 |
rr5 = 3.0d0 * rr3 / r2 |
| 24 |
rr7 = 5.0d0 * rr5 / r2 |
| 25 |
rr9 = 7.0d0 * rr7 / r2 |
| 26 |
rr11 = 9.0d0 * rr9 / r2 |
| 27 |
|
| 28 |
! construct auxiliary vectors |
| 29 |
|
| 30 |
qir(1) = qi(1)*xji + qi(4)*yji + qi(7)*zji |
| 31 |
qir(2) = qi(2)*xji + qi(5)*yji + qi(8)*zji |
| 32 |
qir(3) = qi(3)*xji + qi(6)*yji + qi(9)*zji |
| 33 |
qjr(1) = qj(1)*xji + qj(4)*yji + qj(7)*zji |
| 34 |
qjr(2) = qj(2)*xji + qj(5)*yji + qj(8)*zji |
| 35 |
qjr(3) = qj(3)*xji + qj(6)*yji + qj(9)*zji |
| 36 |
|
| 37 |
! get intermediate variables for energy terms |
| 38 |
|
| 39 |
sc(2) = di(1)*dj(1) + di(2)*dj(2) + di(3)*dj(3) |
| 40 |
sc(3) = di(1)*xji + di(2)*yji + di(3)*zji |
| 41 |
sc(4) = dj(1)*xji + dj(2)*yji + dj(3)*zji |
| 42 |
sc(5) = qir(1)*xji + qir(2)*yji + qir(3)*zji |
| 43 |
sc(6) = qjr(1)*xji + qjr(2)*yji + qjr(3)*zji |
| 44 |
sc(7) = qir(1)*dj(1) + qir(2)*dj(2) + qir(3)*dj(3) |
| 45 |
sc(8) = qjr(1)*di(1) + qjr(2)*di(2) + qjr(3)*di(3) |
| 46 |
sc(9) = qir(1)*qjr(1) + qir(2)*qjr(2) + qir(3)*qjr(3) |
| 47 |
sc(10) = qi(1)*qj(1) + qi(2)*qj(2) + qi(3)*qj(3) & |
| 48 |
+ qi(4)*qj(4) + qi(5)*qj(5) + qi(6)*qj(6) & |
| 49 |
+ qi(7)*qj(7) + qi(8)*qj(8) + qi(9)*qj(9) |
| 50 |
|
| 51 |
! calculate the gl functions for potential energy |
| 52 |
|
| 53 |
gl(0) = ci*cj |
| 54 |
gl(1) = cj*sc(3) - ci*sc(4) |
| 55 |
gl(2) = ci*sc(6) + cj*sc(5) - sc(3)*sc(4) |
| 56 |
gl(3) = sc(3)*sc(6) - sc(4)*sc(5) |
| 57 |
gl(4) = sc(5)*sc(6) |
| 58 |
gl(6) = sc(2) |
| 59 |
gl(7) = 2.0d0 * (sc(7)-sc(8)) |
| 60 |
gl(8) = 2.0d0 * sc(10) |
| 61 |
gl(5) = -4.0d0 * sc(9) |
| 62 |
|
| 63 |
! get the permanent multipole energy |
| 64 |
|
| 65 |
e = rr1*gl(0) + rr3*(gl(1)+gl(6)) & |
| 66 |
+ rr5*(gl(2)+gl(7)+gl(8)) & |
| 67 |
+ rr7*(gl(3)+gl(5)) + rr9*gl(4) |
| 68 |
|
| 69 |
! intermediate variables for the multipole terms |
| 70 |
|
| 71 |
gf(1) = rr3*gl(0) + rr5*(gl(1)+gl(6)) & |
| 72 |
+ rr7*(gl(2)+gl(7)+gl(8)) & |
| 73 |
+ rr9*(gl(3)+gl(5)) + rr11*gl(4) |
| 74 |
gf(2) = -cj*rr3 + sc(4)*rr5 - sc(6)*rr7 |
| 75 |
gf(3) = ci*rr3 + sc(3)*rr5 + sc(5)*rr7 |
| 76 |
gf(4) = 2.0d0 * rr5 |
| 77 |
gf(5) = 2.0d0 * (-cj*rr5+sc(4)*rr7-sc(6)*rr9) |
| 78 |
gf(6) = 2.0d0 * (-ci*rr5-sc(3)*rr7-sc(5)*rr9) |
| 79 |
gf(7) = 4.0d0 * rr7 |
| 80 |
|
| 81 |
! get the force |
| 82 |
|
| 83 |
ftm2(1) = gf(1)*xji + gf(2)*di(1) + gf(3)*dj(1) & |
| 84 |
+ gf(4)*(qjdi(1)-qidj(1)) + gf(5)*qir(1) & |
| 85 |
+ gf(6)*qjr(1) + gf(7)*(qiqjr(1)+qjqir(1)) |
| 86 |
ftm2(2) = gf(1)*yji + gf(2)*di(2) + gf(3)*dj(2) & |
| 87 |
+ gf(4)*(qjdi(2)-qidj(2)) + gf(5)*qir(2) & |
| 88 |
+ gf(6)*qjr(2) + gf(7)*(qiqjr(2)+qjqir(2)) |
| 89 |
ftm2(3) = gf(1)*zji + gf(2)*di(3) + gf(3)*dj(3) & |
| 90 |
+ gf(4)*(qjdi(3)-qidj(3)) + gf(5)*qir(3) & |
| 91 |
+ gf(6)*qjr(3) + gf(7)*(qiqjr(3)+qjqir(3)) |
| 92 |
|
| 93 |
! now perform the torque calculation |
| 94 |
|
| 95 |
! get the interaction torques |
| 96 |
|
| 97 |
ttm2(1) = -rr3*dixdj(1) + gf(2)*dixr(1)-gf(5)*rxqir(1) & |
| 98 |
+ gf(4)*(dixqjr(1)+djxqir(1)+rxqidj(1)-2.0d0*qixqj(1))& |
| 99 |
- gf(7)*(rxqijr(1)+qjrxqir(1)) |
| 100 |
ttm2(2) = -rr3*dixdj(2) + gf(2)*dixr(2)-gf(5)*rxqir(2) & |
| 101 |
+ gf(4)*(dixqjr(2)+djxqir(2)+rxqidj(2)-2.0d0*qixqj(2))& |
| 102 |
- gf(7)*(rxqijr(2)+qjrxqir(2)) |
| 103 |
ttm2(3) = -rr3*dixdj(3) + gf(2)*dixr(3)-gf(5)*rxqir(3) & |
| 104 |
+ gf(4)*(dixqjr(3)+djxqir(3)+rxqidj(3)-2.0d0*qixqj(3))& |
| 105 |
- gf(7)*(rxqijr(3)+qjrxqir(3)) |
| 106 |
ttm3(1) = rr3*dixdj(1) + gf(3)*djxr(1) -gf(6)*rxqjr(1) & |
| 107 |
- gf(4)*(dixqjr(1)+djxqir(1)+rxqjdi(1)-2.0d0*qixqj(1))& |
| 108 |
- gf(7)*(rxqjir(1)-qjrxqir(1)) |
| 109 |
ttm3(2) = rr3*dixdj(2) + gf(3)*djxr(2) -gf(6)*rxqjr(2) & |
| 110 |
- gf(4)*(dixqjr(2)+djxqir(2)+rxqjdi(2)-2.0d0*qixqj(2))& |
| 111 |
- gf(7)*(rxqjir(2)-qjrxqir(2)) |
| 112 |
ttm3(3) = rr3*dixdj(3) + gf(3)*djxr(3) -gf(6)*rxqjr(3) & |
| 113 |
- gf(4)*(dixqjr(3)+djxqir(3)+rxqjdi(3)-2.0d0*qixqj(3))& |
| 114 |
- gf(7)*(rxqjir(3)-qjrxqir(3)) |
| 115 |
|
| 116 |
! update force and torque on site j |
| 117 |
|
| 118 |
ftm1(1,j) = ftm1(1,j) + ftm2(1) |
| 119 |
ftm1(2,j) = ftm1(2,j) + ftm2(2) |
| 120 |
ftm1(3,j) = ftm1(3,j) + ftm2(3) |
| 121 |
ttm1(1,j) = ttm1(1,j) + ttm3(1) |
| 122 |
ttm1(2,j) = ttm1(2,j) + ttm3(2) |
| 123 |
ttm1(3,j) = ttm1(3,j) + ttm3(3) |
| 124 |
|
| 125 |
! update force and torque on site i |
| 126 |
|
| 127 |
ftm1(1,i) = ftm1(1,i) - ftm2(1) |
| 128 |
ftm1(2,i) = ftm1(2,i) - ftm2(2) |
| 129 |
ftm1(3,i) = ftm1(3,i) - ftm2(3) |
| 130 |
ttm1(1,i) = ttm1(1,i) + ttm2(1) |
| 131 |
ttm1(2,i) = ttm1(2,i) + ttm2(2) |
| 132 |
ttm1(3,i) = ttm1(3,i) + ttm2(3) |
