1 |
chrisfen |
2919 |
subroutine doMultipolePair(atom1, atom2, d, rij, r2, rcut, sw, & |
2 |
|
|
vpair, fpair, pot, eFrame, f, t, do_pot) |
3 |
|
|
|
4 |
|
|
!*************************************************************** |
5 |
|
|
! doMultipolePair evaluates the potential, forces, and torques |
6 |
|
|
! between point multipoles. It is based on the ewald2 real |
7 |
|
|
! space routine in the MDMULP code written by W. Smith and |
8 |
|
|
! available from the CCP5 website. |
9 |
|
|
! |
10 |
|
|
! We're using the damped real space portion of the Ewald sum |
11 |
|
|
! as the entire interaction. Details on this portion of the |
12 |
|
|
! sum can be found in: |
13 |
|
|
! |
14 |
|
|
! "Point Multipoles in the Ewald Summation (Revisited)," by |
15 |
|
|
! W. Smith, CCP5 Newsletter, 46, pp. 18-30 (1998). |
16 |
|
|
! |
17 |
|
|
!************************************************************** |
18 |
|
|
|
19 |
|
|
logical, intent(in) :: do_pot |
20 |
|
|
|
21 |
|
|
integer, intent(in) :: atom1, atom2 |
22 |
|
|
real(kind=dp), intent(in) :: rij, r2, sw, rcut |
23 |
|
|
real(kind=dp), intent(in), dimension(3) :: d |
24 |
|
|
real(kind=dp), intent(inout) :: vpair |
25 |
|
|
real(kind=dp), intent(inout), dimension(3) :: fpair |
26 |
|
|
|
27 |
|
|
real( kind = dp ) :: pot |
28 |
|
|
real( kind = dp ), dimension(9,nLocal) :: eFrame |
29 |
|
|
real( kind = dp ), dimension(3,nLocal) :: f |
30 |
|
|
real( kind = dp ), dimension(3,nLocal) :: t |
31 |
|
|
logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole |
32 |
|
|
logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole |
33 |
|
|
integer :: me1, me2, id1, id2 |
34 |
|
|
real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j |
35 |
|
|
real (kind=dp) :: qxx_i, qyy_i, qzz_i |
36 |
|
|
real (kind=dp) :: qxx_j, qyy_j, qzz_j |
37 |
|
|
|
38 |
|
|
real (kind=dp) :: epot, qii, alsq2, alsq2n, exp2a, ralpha |
39 |
|
|
real (kind=dp), dimension(3) :: ftm, ttm_i, ttm_j |
40 |
|
|
real (kind=dp), dimension(3) :: qir, qjr, qidj, qjdi, qiqjr, qjqir |
41 |
|
|
real (kind=dp), dimension(3) :: dixdj, qixqj, dixr, djxr, rxqir, rxqjr |
42 |
|
|
real (kind=dp), dimension(3) :: dixqjr, djxqir, rxqidj, rxqjdi, rxqijr |
43 |
|
|
real (kind=dp), dimension(3) :: rxqjir, qjrxqir |
44 |
|
|
real (kind=dp), dimension(6) :: bn |
45 |
|
|
real (kind=dp), dimension(10) :: sc |
46 |
|
|
real (kind=dp), dimension(9) :: gl |
47 |
|
|
|
48 |
|
|
real(kind=dp) :: a1, a2, a3, a4, a5, p |
49 |
|
|
|
50 |
|
|
DATA A1,A2,A3/0.254829592,-0.284496736,1.421413741/ |
51 |
|
|
DATA A4,A5,P/-1.453152027,1.061405429,0.3275911/ |
52 |
|
|
|
53 |
|
|
if (.not.summationMethodChecked) then |
54 |
|
|
call checkSummationMethod() |
55 |
|
|
endif |
56 |
|
|
|
57 |
|
|
#ifdef IS_MPI |
58 |
|
|
me1 = atid_Row(atom1) |
59 |
|
|
me2 = atid_Col(atom2) |
60 |
|
|
#else |
61 |
|
|
me1 = atid(atom1) |
62 |
|
|
me2 = atid(atom2) |
63 |
|
|
#endif |
64 |
|
|
|
65 |
|
|
! logicals |
66 |
|
|
i_is_Charge = ElectrostaticMap(me1)%is_Charge |
67 |
|
|
i_is_Dipole = ElectrostaticMap(me1)%is_Dipole |
68 |
|
|
i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole |
69 |
|
|
|
70 |
|
|
j_is_Charge = ElectrostaticMap(me2)%is_Charge |
71 |
|
|
j_is_Dipole = ElectrostaticMap(me2)%is_Dipole |
72 |
|
|
j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole |
73 |
|
|
|
74 |
|
|
if (i_is_Charge) then |
75 |
|
|
q_i = ElectrostaticMap(me1)%charge |
76 |
|
|
else |
77 |
|
|
q_i = 0.0_dp |
78 |
|
|
endif |
79 |
|
|
|
80 |
|
|
if (j_is_Charge) then |
81 |
|
|
q_j = ElectrostaticMap(me2)%charge |
82 |
|
|
else |
83 |
|
|
q_j = 0.0_dp |
84 |
|
|
endif |
85 |
|
|
|
86 |
|
|
if (i_is_Dipole) then |
87 |
|
|
mu_i = ElectrostaticMap(me1)%dipole_moment |
88 |
|
|
#ifdef IS_MPI |
89 |
|
|
d_i(1) = eFrame_Row(3,atom1) |
90 |
|
|
d_i(2) = eFrame_Row(6,atom1) |
91 |
|
|
d_i(3) = eFrame_Row(9,atom1) |
92 |
|
|
#else |
93 |
|
|
d_i(1) = eFrame(3,atom1) |
94 |
|
|
d_i(2) = eFrame(6,atom1) |
95 |
|
|
d_i(3) = eFrame(9,atom1) |
96 |
|
|
#endif |
97 |
|
|
d_i = d_i * mu_i |
98 |
|
|
else |
99 |
|
|
d_i = 0.0_dp |
100 |
|
|
endif |
101 |
|
|
|
102 |
|
|
if (j_is_Dipole) then |
103 |
|
|
mu_j = ElectrostaticMap(me2)%dipole_moment |
104 |
|
|
#ifdef IS_MPI |
105 |
|
|
d_j(1) = eFrame_Col(3,atom2) |
106 |
|
|
d_j(2) = eFrame_Col(6,atom2) |
107 |
|
|
d_j(3) = eFrame_Col(9,atom2) |
108 |
|
|
#else |
109 |
|
|
d_j(1) = eFrame(3,atom2) |
110 |
|
|
d_j(2) = eFrame(6,atom2) |
111 |
|
|
d_j(3) = eFrame(9,atom2) |
112 |
|
|
#endif |
113 |
|
|
d_j = d_j * mu_j |
114 |
|
|
else |
115 |
|
|
d_j = 0.0_dp |
116 |
|
|
endif |
117 |
|
|
|
118 |
|
|
if (i_is_Quadrupole) then |
119 |
|
|
qxx_i = ElectrostaticMap(me1)%quadrupole_moments(1) |
120 |
|
|
qyy_i = ElectrostaticMap(me1)%quadrupole_moments(2) |
121 |
|
|
qzz_i = ElectrostaticMap(me1)%quadrupole_moments(3) |
122 |
|
|
#ifdef IS_MPI |
123 |
|
|
qpole_i(1) = qxx_i * eFrame_Row(1,atom1) |
124 |
|
|
qpole_i(2) = qxx_i * eFrame_Row(4,atom1) |
125 |
|
|
qpole_i(3) = qxx_i * eFrame_Row(7,atom1) |
126 |
|
|
qpole_i(4) = qyy_i * eFrame_Row(2,atom1) |
127 |
|
|
qpole_i(5) = qyy_i * eFrame_Row(5,atom1) |
128 |
|
|
qpole_i(6) = qyy_i * eFrame_Row(8,atom1) |
129 |
|
|
qpole_i(7) = qzz_i * eFrame_Row(3,atom1) |
130 |
|
|
qpole_i(8) = qzz_i * eFrame_Row(6,atom1) |
131 |
|
|
qpole_i(9) = qzz_i * eFrame_Row(9,atom1) |
132 |
|
|
#else |
133 |
|
|
qpole_i(1) = qxx_i * eFrame(1,atom1) |
134 |
|
|
qpole_i(2) = qxx_i * eFrame(4,atom1) |
135 |
|
|
qpole_i(3) = qxx_i * eFrame(7,atom1) |
136 |
|
|
qpole_i(4) = qyy_i * eFrame(2,atom1) |
137 |
|
|
qpole_i(5) = qyy_i * eFrame(5,atom1) |
138 |
|
|
qpole_i(6) = qyy_i * eFrame(8,atom1) |
139 |
|
|
qpole_i(7) = qzz_i * eFrame(3,atom1) |
140 |
|
|
qpole_i(8) = qzz_i * eFrame(6,atom1) |
141 |
|
|
qpole_i(9) = qzz_i * eFrame(9,atom1) |
142 |
|
|
#endif |
143 |
|
|
else |
144 |
|
|
qpole_i = 0.0_dp |
145 |
|
|
endif |
146 |
|
|
|
147 |
|
|
if (j_is_Quadrupole) then |
148 |
|
|
qxx_j = ElectrostaticMap(me2)%quadrupole_moments(1) |
149 |
|
|
qyy_j = ElectrostaticMap(me2)%quadrupole_moments(2) |
150 |
|
|
qzz_j = ElectrostaticMap(me2)%quadrupole_moments(3) |
151 |
|
|
#ifdef IS_MPI |
152 |
|
|
qpole_j(1) = qxx_j * eFrame_Col(1,atom2) |
153 |
|
|
qpole_j(2) = qxx_j * eFrame_Col(4,atom2) |
154 |
|
|
qpole_j(3) = qxx_j * eFrame_Col(7,atom2) |
155 |
|
|
qpole_j(4) = qyy_j * eFrame_Col(2,atom2) |
156 |
|
|
qpole_j(5) = qyy_j * eFrame_Col(5,atom2) |
157 |
|
|
qpole_j(6) = qyy_j * eFrame_Col(8,atom2) |
158 |
|
|
qpole_j(7) = qzz_j * eFrame_Col(3,atom2) |
159 |
|
|
qpole_j(8) = qzz_j * eFrame_Col(6,atom2) |
160 |
|
|
qpole_j(9) = qzz_j * eFrame_Col(9,atom2) |
161 |
|
|
#else |
162 |
|
|
qpole_j(1) = qxx_j * eFrame(1,atom2) |
163 |
|
|
qpole_j(2) = qxx_j * eFrame(4,atom2) |
164 |
|
|
qpole_j(3) = qxx_j * eFrame(7,atom2) |
165 |
|
|
qpole_j(4) = qyy_j * eFrame(2,atom2) |
166 |
|
|
qpole_j(5) = qyy_j * eFrame(5,atom2) |
167 |
|
|
qpole_j(6) = qyy_j * eFrame(8,atom2) |
168 |
|
|
qpole_j(7) = qzz_j * eFrame(3,atom2) |
169 |
|
|
qpole_j(8) = qzz_j * eFrame(6,atom2) |
170 |
|
|
qpole_j(9) = qzz_j * eFrame(9,atom2) |
171 |
|
|
#endif |
172 |
|
|
else |
173 |
|
|
qpole_j = 0.0_dp |
174 |
|
|
endif |
175 |
|
|
|
176 |
|
|
! |
177 |
|
|
! INITIALISE TEMPORARY FORCE AND TORQUE ACCUMULATORS |
178 |
|
|
! |
179 |
|
|
|
180 |
|
|
ftm = 0.0 |
181 |
|
|
ttm_i = 0.0 |
182 |
|
|
ttm_j = 0.0 |
183 |
|
|
|
184 |
|
|
ri = 1.0_dp / rij |
185 |
|
|
ri2 = ri * ri |
186 |
|
|
|
187 |
|
|
! |
188 |
|
|
! CALCULATE BN FUNCTIONS |
189 |
|
|
! |
190 |
|
|
if (screeningMethod .eq. DAMPED) then |
191 |
|
|
! assemble the damping variables |
192 |
|
|
call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal) |
193 |
|
|
else |
194 |
|
|
dampingAlpha = 0.0_dp |
195 |
|
|
erfcVal = 1.0_dp |
196 |
|
|
endif |
197 |
|
|
|
198 |
|
|
rrtpi = 1.0 / sqrt(pi) |
199 |
|
|
alsq2 = 2.0 * dampingAlpha**2 |
200 |
|
|
ralpi = 0.0 |
201 |
|
|
if(dampingAlpha .gt. 0.0) ralpi = rrtpi / dampingAlpha |
202 |
|
|
|
203 |
|
|
alphar = dampingAlpha * rij |
204 |
|
|
t = 1.0 / (1.0 + p*alphar) |
205 |
|
|
exp2a = exp(-alphar**2) |
206 |
|
|
bn(1) = ((((a5*t+a4)*t+a3)*t+a2)*t+a1) * t * exp2a * ri |
207 |
|
|
alsq2n = ralpi |
208 |
|
|
|
209 |
|
|
do n = 1, 5 |
210 |
|
|
bfac = real(n+n-1, kind=dp) |
211 |
|
|
alsq2n = alsq2*alsq2n |
212 |
|
|
bn(n+1) = ri2 * (bfac*bn(n) + alsq2n*exp2a) |
213 |
|
|
enddo |
214 |
|
|
|
215 |
|
|
ralpha = dampingAlpha * rij |
216 |
|
|
bn(0) = erfcVal * ri |
217 |
|
|
alsq2 = 2.0d0 * dampingAlpha**2 |
218 |
|
|
alsq2n = 0.0d0 |
219 |
|
|
if (dampingAlpha .gt. 0.0_dp) alsq2n = 1.0_dp / (sqrtpi*dampingAlpha) |
220 |
|
|
exp2a = exp(-ralpha**2) |
221 |
|
|
do m = 1, 5 |
222 |
|
|
bfac = real(m+m-1, kind=dp) |
223 |
|
|
alsq2n = alsq2 * alsq2n |
224 |
|
|
bn(m) = (bfac*bn(m-1)+alsq2n*exp2a) / r2 |
225 |
|
|
end do |
226 |
|
|
|
227 |
|
|
! |
228 |
|
|
! CONSTRUCT AUXILIARY VECTORS |
229 |
|
|
! |
230 |
|
|
|
231 |
|
|
dixdj(1) = d_i(2)*d_j(3) - d_i(3)*d_j(2) |
232 |
|
|
dixdj(2) = d_i(3)*d_j(1) - d_i(1)*d_j(3) |
233 |
|
|
dixdj(3) = d_i(1)*d_j(2) - d_i(2)*d_j(1) |
234 |
|
|
|
235 |
|
|
dixr(1) = d_i(2)*d(3) - d_i(3)*d(2) |
236 |
|
|
dixr(2) = d_i(3)*d(1) - d_i(1)*d(3) |
237 |
|
|
dixr(3) = d_i(1)*d(2) - d_i(2)*d(1) |
238 |
|
|
|
239 |
|
|
djxr(1) = d_j(2)*d(3) - d_j(3)*d(2) |
240 |
|
|
djxr(2) = d_j(3)*d(1) - d_j(1)*d(3) |
241 |
|
|
djxr(3) = d_j(1)*d(2) - d_j(2)*d(1) |
242 |
|
|
|
243 |
|
|
qir(1) = qpole_i(1)*d(1) + qpole_i(4)*d(2) + qpole_i(7)*d(3) |
244 |
|
|
qir(2) = qpole_i(2)*d(1) + qpole_i(5)*d(2) + qpole_i(8)*d(3) |
245 |
|
|
qir(3) = qpole_i(3)*d(1) + qpole_i(6)*d(2) + qpole_i(9)*d(3) |
246 |
|
|
|
247 |
|
|
qjr(1) = qpole_j(1)*d(1) + qpole_j(4)*d(2) + qpole_j(7)*d(3) |
248 |
|
|
qjr(2) = qpole_j(2)*d(1) + qpole_j(5)*d(2) + qpole_j(8)*d(3) |
249 |
|
|
qjr(3) = qpole_j(3)*d(1) + qpole_j(6)*d(2) + qpole_j(9)*d(3) |
250 |
|
|
|
251 |
|
|
qiqjr(1) = qpole_i(1)*qjr(1) + qpole_i(4)*qjr(2) + qpole_i(7)*qjr(3) |
252 |
|
|
qiqjr(2) = qpole_i(2)*qjr(1) + qpole_i(5)*qjr(2) + qpole_i(8)*qjr(3) |
253 |
|
|
qiqjr(3) = qpole_i(3)*qjr(1) + qpole_i(6)*qjr(2) + qpole_i(9)*qjr(3) |
254 |
|
|
|
255 |
|
|
|
256 |
|
|
qjqir(1) = qpole_j(1)*QIR(1) + qpole_j(4)*QIR(2) + qpole_j(7)*QIR(3) |
257 |
|
|
qjqir(2) = qpole_j(2)*QIR(1) + qpole_j(5)*QIR(2) + qpole_j(8)*QIR(3) |
258 |
|
|
qjqir(3) = qpole_j(3)*QIR(1) + qpole_j(6)*QIR(2) + qpole_j(9)*QIR(3) |
259 |
|
|
|
260 |
|
|
qixqj(1) = qpole_i(2)*qpole_j(3) + qpole_i(5)*qpole_j(6) + & |
261 |
|
|
qpole_i(8)*qpole_j(9) - qpole_i(3)*qpole_j(2) - & |
262 |
|
|
qpole_i(6)*qpole_j(5) - qpole_i(9)*qpole_j(8) |
263 |
|
|
|
264 |
|
|
qixqj(2) = qpole_i(3)*qpole_j(1) + qpole_i(6)*qpole_j(4) + & |
265 |
|
|
qpole_i(9)*qpole_j(7) - qpole_i(1)*qpole_j(3) - & |
266 |
|
|
qpole_i(4)*qpole_j(6) - qpole_i(7)*qpole_j(9) |
267 |
|
|
|
268 |
|
|
qixqj(3) = qpole_i(1)*qpole_j(2) + qpole_i(4)*qpole_j(5) + & |
269 |
|
|
qpole_i(7)*qpole_j(8) - qpole_i(2)*qpole_j(1) - & |
270 |
|
|
qpole_i(5)*qpole_j(4) - qpole_i(8)*qpole_j(7) |
271 |
|
|
|
272 |
|
|
rxqir(1) = d(2)*qir(3) - d(3)*qir(2) |
273 |
|
|
rxqir(2) = d(3)*qir(1) - d(1)*qir(3) |
274 |
|
|
rxqir(3) = d(1)*qir(2) - d(2)*qir(1) |
275 |
|
|
|
276 |
|
|
rxqjr(1) = d(2)*qjr(3) - d(3)*qjr(2) |
277 |
|
|
rxqjr(2) = d(3)*qjr(1) - d(1)*qjr(3) |
278 |
|
|
rxqjr(3) = d(1)*qjr(2) - d(2)*qjr(1) |
279 |
|
|
|
280 |
|
|
rxqijr(1) = d(2)*qiqjr(3) - d(3)*qiqjr(2) |
281 |
|
|
rxqijr(2) = d(3)*qiqjr(1) - d(1)*qiqjr(3) |
282 |
|
|
rxqijr(3) = d(1)*qiqjr(2) - d(2)*qiqjr(1) |
283 |
|
|
|
284 |
|
|
rxqjir(1) = d(2)*qjqir(3) - d(3)*qjqir(2) |
285 |
|
|
rxqjir(2) = d(3)*qjqir(1) - d(1)*qjqir(3) |
286 |
|
|
rxqjir(3) = d(1)*qjqir(2) - d(2)*qjqir(1) |
287 |
|
|
|
288 |
|
|
qjrxqir(1) = qjr(2)*qir(3) - qjr(3)*qir(2) |
289 |
|
|
qjrxqir(2) = qjr(3)*qir(1) - qjr(1)*qir(3) |
290 |
|
|
qjrxqir(3) = qjr(1)*qir(2) - qjr(2)*qir(1) |
291 |
|
|
|
292 |
|
|
qidj(1) = qpole_i(1)*d_j(1) + qpole_i(4)*d_j(2) + qpole_i(7)*d_j(3) |
293 |
|
|
qidj(2) = qpole_i(2)*d_j(1) + qpole_i(5)*d_j(2) + qpole_i(8)*d_j(3) |
294 |
|
|
qidj(3) = qpole_i(3)*d_j(1) + qpole_i(6)*d_j(2) + qpole_i(9)*d_j(3) |
295 |
|
|
|
296 |
|
|
qjdi(1) = qpole_j(1)*d_i(1) + qpole_j(4)*d_i(2) + qpole_j(7)*d_i(3) |
297 |
|
|
qjdi(2) = qpole_j(2)*d_i(1) + qpole_j(5)*d_i(2) + qpole_j(8)*d_i(3) |
298 |
|
|
qjdi(3) = qpole_j(3)*d_i(1) + qpole_j(6)*d_i(2) + qpole_j(9)*d_i(3) |
299 |
|
|
|
300 |
|
|
dixqjr(1) = d_i(2)*qjr(3) - d_i(3)*qjr(2) |
301 |
|
|
dixqjr(2) = d_i(3)*qjr(1) - d_i(1)*qjr(3) |
302 |
|
|
dixqjr(3) = d_i(1)*qjr(2) - d_i(2)*qjr(1) |
303 |
|
|
|
304 |
|
|
djxqir(1) = d_j(2)*qir(3) - d_j(3)*qir(2) |
305 |
|
|
djxqir(2) = d_j(3)*qir(1) - d_j(1)*qir(3) |
306 |
|
|
djxqir(3) = d_j(1)*qir(2) - d_j(2)*qir(1) |
307 |
|
|
|
308 |
|
|
rxqidj(1) = d(2)*qidj(3) - d(3)*qidj(2) |
309 |
|
|
rxqidj(2) = d(3)*qidj(1) - d(1)*qidj(3) |
310 |
|
|
rxqidj(3) = d(1)*qidj(2) - d(2)*qidj(1) |
311 |
|
|
|
312 |
|
|
rxqjdi(1) = d(2)*qjdi(3) - d(3)*qjdi(2) |
313 |
|
|
rxqjdi(2) = d(3)*qjdi(1) - d(1)*qjdi(3) |
314 |
|
|
rxqjdi(3) = d(1)*qjdi(2) - d(2)*qjdi(1) |
315 |
|
|
|
316 |
|
|
! |
317 |
|
|
! CALCULATE SCALAR PRODUCTS |
318 |
|
|
! |
319 |
|
|
|
320 |
|
|
qii = qpole_i(1) + qpole_i(5) + qpole_i(9) |
321 |
|
|
|
322 |
|
|
sc(1) = qpole_j(1) + qpole_j(5) + qpole_j(9) |
323 |
|
|
sc(2) = d_i(1)*d_j(1) + d_i(2)*d_j(2) + d_i(3)*d_j(3) |
324 |
|
|
sc(3) = d_i(1)*d(1) + d_i(2)*d(2) + d_i(3)*d(3) |
325 |
|
|
sc(4) = d_j(1)*d(1) + d_j(2)*d(2) + d_j(3)*d(3) |
326 |
|
|
sc(5) = qir(1)*d(1) + qir(2)*d(2) + qir(3)*d(3) |
327 |
|
|
sc(6) = qjr(1)*d(1) + qjr(2)*d(2) + qjr(3)*d(3) |
328 |
|
|
sc(7) = qir(1)*d_j(1) + qir(2)*d_j(2) + qir(3)*d_j(3) |
329 |
|
|
sc(8) = qjr(1)*d_i(1) + qjr(2)*d_i(2) + qjr(3)*d_i(3) |
330 |
|
|
sc(9) = qir(1)*qjr(1) + qir(2)*qjr(2) + qir(3)*qjr(3) |
331 |
|
|
sc(10) = qpole_i(1)*qpole_j(1) + qpole_i(2)*qpole_j(2) + & |
332 |
|
|
qpole_i(3)*qpole_j(3) + qpole_i(4)*qpole_j(4) + & |
333 |
|
|
qpole_i(5)*qpole_j(5) + qpole_i(6)*qpole_j(6) + & |
334 |
|
|
qpole_i(7)*qpole_j(7) + qpole_i(8)*qpole_j(8) + & |
335 |
|
|
qpole_i(9)*qpole_j(9) |
336 |
|
|
|
337 |
|
|
! |
338 |
|
|
! CALCULATE GL FUNCTIONS FOR POTENTIAL |
339 |
|
|
! |
340 |
|
|
|
341 |
|
|
gl(1) = q_i*q_j |
342 |
|
|
gl(2) = q_j*sc(3) - q_i*sc(4) |
343 |
|
|
gl(3) = q_i*sc(6) + q_j*sc(5) - sc(3)*sc(4) |
344 |
|
|
gl(4) = sc(3)*sc(6) - sc(4)*sc(5) |
345 |
|
|
gl(5) = sc(5)*sc(6) |
346 |
|
|
gl(7) = sc(2) - q_j*qii - q_i*sc(1) |
347 |
|
|
gl(8) = sc(4)*qii - sc(1)*sc(3) + 2.0*(sc(7) - sc(8)) |
348 |
|
|
gl(9) = 2.0*sc(10) + sc(1)*qii |
349 |
|
|
gl(6) = -(sc(1)*sc(5) + sc(6)*qii + 4.0*sc(9)) |
350 |
|
|
|
351 |
|
|
! |
352 |
|
|
! CALCULATE POTENTIAL AND VIRIAL |
353 |
|
|
! |
354 |
|
|
|
355 |
|
|
epot = bn(1)*gl(1) + bn(2)*(gl(2) + gl(7)) + bn(3)*(gl(3) & |
356 |
|
|
+ gl(8) + gl(9)) + bn(4)*(gl(4) + gl(6)) + bn(5)*gl(5) |
357 |
|
|
|
358 |
|
|
vpair = vpair + epot |
359 |
|
|
|
360 |
|
|
if (do_pot) then |
361 |
|
|
#ifdef IS_MPI |
362 |
|
|
pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + & |
363 |
|
|
0.5_dp*epot*sw |
364 |
|
|
pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + & |
365 |
|
|
0.5_dp*epot*sw |
366 |
|
|
#else |
367 |
|
|
pot = pot + epot*sw |
368 |
|
|
#endif |
369 |
|
|
endif |
370 |
|
|
|
371 |
|
|
! |
372 |
|
|
! CALCULATE FORCE AND TORQUE COEFFICIENTS |
373 |
|
|
! |
374 |
|
|
|
375 |
|
|
gl(1) = bn(2)*gl(1) + bn(3)*(gl(2) + gl(7)) + bn(4)*(gl(3) & |
376 |
|
|
+ gl(8) + gl(9)) + bn(5)*(gl(4) + gl(6)) + bn(6)*gl(5) |
377 |
|
|
gl(2) = -q_j*bn(2) + (sc(4) + sc(1))*bn(3) - sc(6)*bn(4) |
378 |
|
|
gl(3) = q_i*bn(2) + (sc(3) - qii)*bn(3) + sc(5)*bn(4) |
379 |
|
|
gl(4) = 2.0*bn(3) |
380 |
|
|
gl(5) = 2.0*( - q_j*bn(3) + (sc(1) + sc(4))*bn(4) - sc(6)* bn(5)) |
381 |
|
|
gl(6) = 2.0*( - q_i*bn(3) + (qii - sc(3))*bn(4) - sc(5)*bn(5)) |
382 |
|
|
gl(7) = 4.0*bn(4) |
383 |
|
|
|
384 |
|
|
! |
385 |
|
|
! CALCULATE FORCES BETWEEN MULTIPOLES I AND J |
386 |
|
|
! |
387 |
|
|
|
388 |
|
|
ftm(1) = gl(1)*d(1) + gl(2)*d_i(1) + gl(3)*d_j(1) + & |
389 |
|
|
gl(4)*(qjdi(1) - qidj(1)) + gl(5)*qir(1) + & |
390 |
|
|
gl(6)*qjr(1) + gl(7)*(qiqjr(1) + qjqir(1)) |
391 |
|
|
ftm(2) = gl(1)*d(2) + gl(2)*d_i(2) + gl(3)*d_j(2) + & |
392 |
|
|
gl(4)*(qjdi(2) - qidj(2)) + gl(5)*qir(2) + & |
393 |
|
|
gl(6)*qjr(2) + gl(7)*(qiqjr(2) + qjqir(2)) |
394 |
|
|
ftm(3) = gl(1)*d(3) + gl(2)*d_i(3) + gl(3)*d_j(3) + & |
395 |
|
|
gl(4)*(qjdi(3) - qidj(3)) + gl(5)*qir(3) + & |
396 |
|
|
gl(6)*qjr(3) + gl(7)*(qiqjr(3) + qjqir(3)) |
397 |
|
|
|
398 |
|
|
! |
399 |
|
|
! CALCULATE TORQUES BETWEEN MULTIPOLES I AND J |
400 |
|
|
! |
401 |
|
|
|
402 |
|
|
ttm_i(1) = - bn(2)*dixdj(1) + gl(2)*dixr(1) + gl(4)* & |
403 |
|
|
(dixqjr(1) + djxqir(1) + rxqidj(1) - 2.0*qixqj(1)) - & |
404 |
|
|
gl(5)*rxqir(1) - gl(7)*(rxqijr(1) + qjrxqir(1)) |
405 |
|
|
ttm_i(2) = - bn(2)*dixdj(2) + gl(2)*dixr(2) + gl(4)* & |
406 |
|
|
(dixqjr(2) + djxqir(2) + rxqidj(2) - 2.0*qixqj(2)) - & |
407 |
|
|
gl(5)*rxqir(2) - gl(7)*(rxqijr(2) + qjrxqir(2)) |
408 |
|
|
ttm_i(3) = - bn(2)*dixdj(3) + gl(2)*dixr(3) + gl(4)* & |
409 |
|
|
(dixqjr(3) + djxqir(3) + rxqidj(3) - 2.0*qixqj(3)) - & |
410 |
|
|
gl(5)*rxqir(3) - gl(7)*(rxqijr(3) + qjrxqir(3)) |
411 |
|
|
ttm_j(1) = bn(2)*dixdj(1) + gl(3)*djxr(1) - gl(4)* & |
412 |
|
|
(dixqjr(1) + djxqir(1) + rxqjdi(1) - 2.0*qixqj(1)) - & |
413 |
|
|
gl(6)*rxqjr(1) - gl(7)*(rxqjir(1) - qjrxqir(1)) |
414 |
|
|
ttm_j(2) = bn(2)*dixdj(2) + gl(3)*djxr(2) - gl(4)* & |
415 |
|
|
(dixqjr(2) + djxqir(2) + rxqjdi(2) - 2.0*qixqj(2)) - & |
416 |
|
|
gl(6)*rxqjr(2) - gl(7)*(rxqjir(2) - qjrxqir(2)) |
417 |
|
|
ttm_j(3) = bn(2)*dixdj(3) + gl(3)*djxr(3) - gl(4)* & |
418 |
|
|
(dixqjr(3) + djxqir(3) + rxqjdi(3) - 2.0*qixqj(3)) - & |
419 |
|
|
gl(6)*rxqjr(3) - gl(7)*(rxqjir(3) - qjrxqir(3)) |
420 |
|
|
|
421 |
|
|
ftm = ftm*sw |
422 |
|
|
ttm_i = ttm_i*sw |
423 |
|
|
ttm_j = ttm_j*sw |
424 |
|
|
|
425 |
|
|
#ifdef IS_MPI |
426 |
|
|
f_Row(1,atom1) = f_Row(1,atom1) + ftm(1) |
427 |
|
|
f_Row(2,atom1) = f_Row(2,atom1) + ftm(2) |
428 |
|
|
f_Row(3,atom1) = f_Row(3,atom1) + ftm(3) |
429 |
|
|
|
430 |
|
|
f_Col(1,atom2) = f_Col(1,atom2) - ftm(1) |
431 |
|
|
f_Col(2,atom2) = f_Col(2,atom2) - ftm(2) |
432 |
|
|
f_Col(3,atom2) = f_Col(3,atom2) - ftm(3) |
433 |
|
|
|
434 |
|
|
t_Row(1,atom1) = t_Row(1,atom1) + ttm_i(1) |
435 |
|
|
t_Row(2,atom1) = t_Row(2,atom1) + ttm_i(2) |
436 |
|
|
t_Row(3,atom1) = t_Row(3,atom1) + ttm_i(3) |
437 |
|
|
|
438 |
|
|
t_Col(1,atom2) = t_Col(1,atom2) + ttm_j(1) |
439 |
|
|
t_Col(2,atom2) = t_Col(2,atom2) + ttm_j(2) |
440 |
|
|
t_Col(3,atom2) = t_Col(3,atom2) + ttm_j(3) |
441 |
|
|
#else |
442 |
|
|
f(1,atom1) = f(1,atom1) + ftm(1) |
443 |
|
|
f(2,atom1) = f(2,atom1) + ftm(2) |
444 |
|
|
f(3,atom1) = f(3,atom1) + ftm(3) |
445 |
|
|
|
446 |
|
|
f(1,atom2) = f(1,atom2) - ftm(1) |
447 |
|
|
f(2,atom2) = f(2,atom2) - ftm(2) |
448 |
|
|
f(3,atom2) = f(3,atom2) - ftm(3) |
449 |
|
|
|
450 |
|
|
t(1,atom1) = t(1,atom1) + ttm_i(1) |
451 |
|
|
t(2,atom1) = t(2,atom1) + ttm_i(2) |
452 |
|
|
t(3,atom1) = t(3,atom1) + ttm_i(3) |
453 |
|
|
|
454 |
|
|
t(1,atom2) = t(1,atom2) + ttm_j(1) |
455 |
|
|
t(2,atom2) = t(2,atom2) + ttm_j(2) |
456 |
|
|
t(3,atom2) = t(3,atom2) + ttm_j(3) |
457 |
|
|
#endif |
458 |
|
|
|
459 |
|
|
#ifdef IS_MPI |
460 |
|
|
id1 = AtomRowToGlobal(atom1) |
461 |
|
|
id2 = AtomColToGlobal(atom2) |
462 |
|
|
#else |
463 |
|
|
id1 = atom1 |
464 |
|
|
id2 = atom2 |
465 |
|
|
#endif |
466 |
|
|
|
467 |
|
|
if (molMembershipList(id1) .ne. molMembershipList(id2)) then |
468 |
|
|
|
469 |
|
|
fpair(1) = fpair(1) + ftm(1) |
470 |
|
|
fpair(2) = fpair(2) + ftm(2) |
471 |
|
|
fpair(3) = fpair(3) + ftm(3) |
472 |
|
|
|
473 |
|
|
endif |
474 |
|
|
|
475 |
|
|
return |
476 |
|
|
|
477 |
|
|
end subroutine domultipolepair |