ls1-MarDyn
ls1-MarDyn molecular dynamics code
potforce.h
1#ifndef POTFORCE_H_
2#define POTFORCE_H_
3
6#include "molecules/Comp2Param.h"
7#include "molecules/Molecule.h"
8
18inline void PotForceLJ(const double dr[3], const double& dr2,
19 const double& eps24, const double& sig2,
20 double f[3], double& u6)
21{
22 double invdr2 = 1. / dr2;
23 double lj6 = sig2 * invdr2; lj6 = lj6 * lj6 * lj6;
24 double lj12 = lj6 * lj6;
25 double lj12m6 = lj12 - lj6;
26 u6 = eps24 * lj12m6;
27 double fac = eps24 * (lj12 + lj12m6) * invdr2;
28 for (unsigned short d = 0; d < 3; ++d)
29 f[d] = fac * dr[d];
30}
31
32
36inline void PotForce2Dipole(const double dr[3], const double& dr2, const double* eii, const double* ejj,
37 const double& my2, const double& rffac,
38 double f[3], double m1[3], double m2[3], double& u, double& MyRF)
39{
40 double invdr2 = 1. / dr2;
41 double invdr1 = sqrt(invdr2);
42 double myfac = my2 * invdr2 * invdr1;
43 double costi = 0.;
44 double costj = 0.;
45 double cosgij = 0.;
46 for (unsigned short d = 0; d < 3; ++d) {
47 const double& drd = dr[d];
48 const double& eiid = eii[d];
49 const double& ejjd = ejj[d];
50 costi += eiid * drd;
51 costj += ejjd * drd;
52 cosgij += eiid * ejjd;
53 }
54 costi *= invdr1;
55 costj *= invdr1;
56 u = myfac * (cosgij - 3. * costi * costj);
57 MyRF -= rffac * cosgij;
58 const double partialRijInvdr1 = -3. * u * invdr2;
59 const double partialTiInvdr1 = -myfac * 3. * costj * invdr1;
60 const double partialTjInvdr1 = -myfac * 3. * costi * invdr1;
61 const double& partialGij = myfac;
62 const double fac = -partialRijInvdr1 + (costi*partialTiInvdr1 + costj*partialTjInvdr1) * invdr1;
63 for (unsigned short d = 0; d < 3; ++d)
64 f[d] = fac * dr[d] - partialTiInvdr1*eii[d] - partialTjInvdr1*ejj[d];
65
66 double eiXej[3], eXrij[3];
67 eiXej[0] = eii[1] * ejj[2] - eii[2] * ejj[1];
68 eiXej[1] = eii[2] * ejj[0] - eii[0] * ejj[2];
69 eiXej[2] = eii[0] * ejj[1] - eii[1] * ejj[0];
70 eXrij[0] = eii[1] * dr[2] - eii[2] * dr[1];
71 eXrij[1] = eii[2] * dr[0] - eii[0] * dr[2];
72 eXrij[2] = eii[0] * dr[1] - eii[1] * dr[0];
73 for (unsigned short d = 0; d < 3; ++d)
74 m1[d] = -partialTiInvdr1 * eXrij[d] + (-partialGij + rffac) * eiXej[d];
75 eXrij[0] = ejj[1] * dr[2] - ejj[2] * dr[1];
76 eXrij[1] = ejj[2] * dr[0] - ejj[0] * dr[2];
77 eXrij[2] = ejj[0] * dr[1] - ejj[1] * dr[0];
78 for (unsigned short d = 0; d < 3; ++d)
79 m2[d] = -partialTjInvdr1 * eXrij[d] + (partialGij - rffac) * eiXej[d];
80}
81
82
86inline void PotForce2Quadrupole(const double dr[3], const double& dr2, const double* eii, const double* ejj,
87 const double& q2075,
88 double f[3], double m1[3], double m2[3], double& u)
89{
90 double invdr2 = 1. / dr2;
91 double invdr1 = sqrt(invdr2);
92 double qfac = q2075 * invdr2 * invdr2 * invdr1;
93 double costi = 0.;
94 double costj = 0.;
95 double cosgij = 0.;
96 for (unsigned short d = 0; d < 3; ++d) {
97 const double& drd = dr[d];
98 const double& eiid = eii[d];
99 const double& ejjd = ejj[d];
100 costi += eiid * drd;
101 costj += ejjd * drd;
102 cosgij += eiid * ejjd;
103 }
104 costi *= invdr1;
105 costj *= invdr1;
106 double cos2ti = costi * costi;
107 double cos2tj = costj * costj;
108 double term = (cosgij - 5. * costi * costj);
109 u = qfac * (1. - 5. * (cos2ti + cos2tj) - 15. * cos2ti * cos2tj + 2. * term * term);
110 const double partialRijInvdr1 = -5. * u * invdr2;
111 const double partialTiInvdr1 = -qfac * 10. * (costi + 3. * costi * cos2tj + 2. * costj * term) * invdr1;
112 const double partialTjInvdr1 = -qfac * 10. * (costj + 3. * cos2ti * costj + 2. * costi * term) * invdr1;
113 const double partialGij = qfac * 4. * term;
114 const double fac = -partialRijInvdr1 + (costi*partialTiInvdr1 + costj*partialTjInvdr1) * invdr1;
115 for (unsigned short d = 0; d < 3; ++d)
116 f[d] = fac * dr[d] - partialTiInvdr1*eii[d] - partialTjInvdr1*ejj[d];
117
118 double eiXej[3], eXrij[3];
119 eiXej[0] = eii[1] * ejj[2] - eii[2] * ejj[1];
120 eiXej[1] = eii[2] * ejj[0] - eii[0] * ejj[2];
121 eiXej[2] = eii[0] * ejj[1] - eii[1] * ejj[0];
122 eXrij[0] = eii[1] * dr[2] - eii[2] * dr[1];
123 eXrij[1] = eii[2] * dr[0] - eii[0] * dr[2];
124 eXrij[2] = eii[0] * dr[1] - eii[1] * dr[0];
125 for (unsigned short d = 0; d < 3; ++d)
126 m1[d] = -partialTiInvdr1 * eXrij[d] - partialGij * eiXej[d];
127
128 eXrij[0] = ejj[1] * dr[2] - ejj[2] * dr[1];
129 eXrij[1] = ejj[2] * dr[0] - ejj[0] * dr[2];
130 eXrij[2] = ejj[0] * dr[1] - ejj[1] * dr[0];
131 for (unsigned short d = 0; d < 3; ++d)
132 m2[d] = -partialTjInvdr1 * eXrij[d] + partialGij * eiXej[d];
133}
134
135
139inline void PotForceDiQuadrupole(const double dr[3], const double& dr2, const double* eii, const double* ejj,
140 const double& myq15,
141 double f[3], double m1[3], double m2[3], double& u)
142{
143 double invdr2 = 1. / dr2;
144 double invdr1 = sqrt(invdr2);
145 double myqfac = myq15 * invdr2 * invdr2;
146 double costi = 0.;
147 double costj = 0.;
148 double cosgij = 0.;
149 for (unsigned short d = 0; d < 3; ++d) {
150 const double& drd = dr[d];
151 const double& eiid = eii[d];
152 const double& ejjd = ejj[d];
153 costi += eiid * drd;
154 costj += ejjd * drd;
155 cosgij += eiid * ejjd;
156 }
157 costi *= invdr1;
158 costj *= invdr1;
159 double cos2tj = costj * costj;
160 u = myqfac * (-costi * (5. * cos2tj - 1.) + 2. * cosgij * costj);
161 const double partialRijInvdr1 = -4. * u * invdr2;
162 const double partialTiInvdr1 = myqfac * (-5. * cos2tj + 1.) * invdr1;
163 const double partialTjInvdr1 = myqfac * 2. * (-5. * costi * costj + cosgij) * invdr1;
164 const double partialGij = myqfac * 2. * costj;
165 const double fac = -partialRijInvdr1 + (costi*partialTiInvdr1 + costj*partialTjInvdr1) * invdr1;
166 for (unsigned short d = 0; d < 3; ++d)
167 f[d] = fac * dr[d] - partialTiInvdr1*eii[d] - partialTjInvdr1*ejj[d];
168
169 double eiXej[3], eXrij[3];
170 eiXej[0] = eii[1] * ejj[2] - eii[2] * ejj[1];
171 eiXej[1] = eii[2] * ejj[0] - eii[0] * ejj[2];
172 eiXej[2] = eii[0] * ejj[1] - eii[1] * ejj[0];
173 eXrij[0] = eii[1] * dr[2] - eii[2] * dr[1];
174 eXrij[1] = eii[2] * dr[0] - eii[0] * dr[2];
175 eXrij[2] = eii[0] * dr[1] - eii[1] * dr[0];
176 for (unsigned short d = 0; d < 3; ++d)
177 m1[d] = -partialTiInvdr1 * eXrij[d] - partialGij * eiXej[d];
178
179 eXrij[0] = ejj[1] * dr[2] - ejj[2] * dr[1];
180 eXrij[1] = ejj[2] * dr[0] - ejj[0] * dr[2];
181 eXrij[2] = ejj[0] * dr[1] - ejj[1] * dr[0];
182 for (unsigned short d = 0; d < 3; ++d)
183 m2[d] = -partialTjInvdr1 * eXrij[d] + partialGij * eiXej[d];
184}
185
186
190inline void PotForce2Charge(const double dr[3], const double& dr2,
191 const double& q1q2per4pie0, double f[3], double& u)
192{
193 double invdr2 = 1.0 / dr2;
194 double invdr = sqrt(invdr2);
195 u = q1q2per4pie0 * invdr;
196 const double fac = u * invdr2;
197 for (unsigned short d = 0; d < 3; d++)
198 f[d] = fac * dr[d];
199}
200
201
205inline void PotForceChargeQuadrupole(const double dr[3], const double& dr2,
206 const double* ejj, const double& qQ05per4pie0,
207 double f[3], double m2[3], double& u)
208{
209 double invdr2 = 1.0 / dr2;
210 double invdr = sqrt(invdr2);
211 double costj = 0;
212 for (unsigned short d = 0; d < 3; d++) {
213 costj += ejj[d] * dr[d];
214 }
215 costj *= invdr;
216 double qQinv4dr3 = qQ05per4pie0 * invdr * invdr2;
217 u = qQinv4dr3 * (3.0 * costj * costj - 1);
218
219 const double partialRijInvdr1 = -3.0 * u * invdr2;
220 const double partialTjInvdr1 = 6.0 * costj * qQinv4dr3 * invdr;
221 const double fac = costj * partialTjInvdr1 * invdr - partialRijInvdr1;
222 for (unsigned short d = 0; d < 3; ++d)
223 f[d] = fac*dr[d] - partialTjInvdr1 * ejj[d];
224
225 double minuseXrij[3];
226 minuseXrij[0] = ejj[2]*dr[1] - ejj[1]*dr[2];
227 minuseXrij[1] = ejj[0]*dr[2] - ejj[2]*dr[0];
228 minuseXrij[2] = ejj[1]*dr[0] - ejj[0]*dr[1];
229 for (unsigned short d = 0; d < 3; ++d)
230 m2[d] = partialTjInvdr1 * minuseXrij[d];
231}
232
233
237inline void PotForceChargeDipole(const double dr[3], const double& dr2,
238 const double* ejj, const double& minusqmyper4pie0,
239 double f[3], double m2[3], double& u)
240{
241 double invdr2 = 1.0 / dr2;
242 double invdr = sqrt(invdr2);
243 double costj = 0;
244 for (unsigned short d = 0; d < 3; d++) {
245 costj += ejj[d] * dr[d];
246 }
247 costj *= invdr;
248 double uInvcostj = minusqmyper4pie0 * invdr2;
249 u = uInvcostj * costj;
250
251 // const double partialRijInvdr1 = -2.0 * u * invdr2;
252 const double partialTjInvdr1 = uInvcostj*invdr;
253 const double fac = 3.0 * u * invdr2;
254 for (unsigned short d = 0; d < 3; ++d)
255 f[d] = fac*dr[d] - partialTjInvdr1 * ejj[d];
256
257 double minuseXrij[3];
258 minuseXrij[0] = ejj[2]*dr[1] - ejj[1]*dr[2];
259 minuseXrij[1] = ejj[0]*dr[2] - ejj[2]*dr[0];
260 minuseXrij[2] = ejj[1]*dr[0] - ejj[0]*dr[1];
261 for (unsigned short d = 0; d < 3; ++d)
262 m2[d] = partialTjInvdr1 * minuseXrij[d];
263}
264
265
282inline void PotForce(Molecule& mi, Molecule& mj, ParaStrm& params, double drm[3], double& Upot6LJ, double& UpotXpoles, double& MyRF, double Virial[3], bool calculateLJ)
283// ???better calc Virial, when molecule forces are calculated:
284// summing up molecule virials instead of site virials???
285{ // Force Calculation
286 double f[3];
287 double u;
288 double drs[3], dr2; // site distance vector & length^2
289 Virial[0]=0.;
290 Virial[1]=0.;
291 Virial[2]=0.;
292 // LJ centers
293 // no LJ interaction between solid atoms of the same component
294
295 const unsigned int nc1 = mi.numLJcenters();
296 const unsigned int nc2 = mj.numLJcenters();
297 for (unsigned int si = 0; si < nc1; ++si) {
298 const std::array<double,3> dii = mi.ljcenter_d_abs(si);
299 for (unsigned int sj = 0; sj < nc2; ++sj) {
300 const std::array<double,3> djj = mj.ljcenter_d_abs(sj);
301 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
302 double eps24;
303 params >> eps24;
304 double sig2;
305 params >> sig2;
306 double shift6;
307 params >> shift6; // must be 0.0 for full LJ
308 if (calculateLJ) {
309 PotForceLJ(drs, dr2, eps24, sig2, f, u);
310 u += shift6;
311
312 mi.Fljcenteradd(si, f);
313 mj.Fljcentersub(sj, f);
314 Upot6LJ += u;
315 for (unsigned short d = 0; d < 3; ++d)
316 Virial[d] += 0.5*drm[d] * f[d];
317 }
318 }
319 }
320
321
322 double m1[3], m2[3]; // angular momenta
323
324 const unsigned ne1 = mi.numCharges();
325 const unsigned ne2 = mj.numCharges();
326 const unsigned int nq1 = mi.numQuadrupoles();
327 const unsigned int nq2 = mj.numQuadrupoles();
328 const unsigned int nd1 = mi.numDipoles();
329 const unsigned int nd2 = mj.numDipoles();
330 for (unsigned si = 0; si < ne1; si++) {
331 const std::array<double,3> dii = mi.charge_d_abs(si);
332 // Charge-Charge
333 for (unsigned sj = 0; sj < ne2; sj++) {
334 const std::array<double,3> djj = mj.charge_d_abs(sj);
335 double q1q2per4pie0; // 4pie0 = 1 in reduced units
336 params >> q1q2per4pie0;
337 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
338 PotForce2Charge(drs, dr2, q1q2per4pie0, f, u);
339
340 mi.Fchargeadd(si, f);
341 mj.Fchargesub(sj, f);
342
343 UpotXpoles += u;
344 for (unsigned short d = 0; d < 3; d++)
345 Virial[d] += 0.5*drm[d] * f[d];
346 }
347 // Charge-Quadrupole
348 for (unsigned sj = 0; sj < nq2; sj++) {
349 const std::array<double,3> djj = mj.quadrupole_d_abs(sj);
350 double qQ05per4pie0; // 4pie0 = 1 in reduced units
351 params >> qQ05per4pie0;
352 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
353 const std::array<double,3> ejj = mj.quadrupole_e(sj);
354 PotForceChargeQuadrupole(drs, dr2, ejj.data(), qQ05per4pie0, f, m2, u);
355
356 mi.Fchargeadd(si, f);
357 mj.Fquadrupolesub(sj, f);
358 mj.Madd(m2);
359
360 UpotXpoles += u;
361 for (unsigned short d = 0; d < 3; d++)
362 Virial[d] += 0.5*drm[d] * f[d];
363 }
364 // Charge-Dipole
365 for (unsigned sj = 0; sj < nd2; sj++) {
366 const std::array<double,3> djj = mj.dipole_d_abs(sj);
367 double minusqmyper4pie0;
368 params >> minusqmyper4pie0;
369 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
370 const std::array<double,3> ejj = mj.dipole_e(sj);
371 PotForceChargeDipole(drs, dr2, ejj.data(), minusqmyper4pie0, f, m2, u);
372
373 mi.Fchargeadd(si, f);
374 mj.Fdipolesub(sj, f);
375 mj.Madd(m2);
376
377 UpotXpoles += u;
378 for (unsigned short d = 0; d < 3; d++)
379 Virial[d] += 0.5*drm[d] * f[d];
380 }
381 }
382 for (unsigned int si = 0; si < nq1; ++si) {
383 const std::array<double,3> dii = mi.quadrupole_d_abs(si);
384 const std::array<double,3> eii = mi.quadrupole_e(si);
385
386 // Quadrupole-Charge
387 for (unsigned sj = 0; sj < ne2; sj++) {
388 const std::array<double,3> djj = mj.charge_d_abs(sj);
389 double qQ05per4pie0; // 4pie0 = 1 in reduced units
390 params >> qQ05per4pie0;
391 minusSiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
392 PotForceChargeQuadrupole(drs, dr2, eii.data(), qQ05per4pie0, f, m1, u);
393
394 mi.Fquadrupolesub(si, f);
395 mj.Fchargeadd(sj, f);
396 mi.Madd(m1);
397
398 UpotXpoles += u;
399 for (unsigned short d = 0; d < 3; d++)
400 Virial[d] -= 0.5*drm[d] * f[d];
401 }
402 // Quadrupole-Quadrupole -------------------
403 for (unsigned int sj = 0; sj < nq2; ++sj) {
404 //double drs[3];
405 const std::array<double,3> djj = mj.quadrupole_d_abs(sj);
406 double q2075;
407 params >> q2075;
408 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
409 const std::array<double,3> ejj = mj.quadrupole_e(sj);
410 PotForce2Quadrupole(drs, dr2, eii.data(), ejj.data(), q2075, f, m1, m2, u);
411
412 mi.Fquadrupoleadd(si, f);
413 mj.Fquadrupolesub(sj, f);
414 mi.Madd(m1);
415 mj.Madd(m2);
416
417 UpotXpoles += u;
418 for (unsigned short d = 0; d < 3; ++d)
419 Virial[d] += 0.5*drm[d] * f[d];
420 }
421 // Quadrupole-Dipole -----------------------
422 for (unsigned int sj = 0; sj < nd2; ++sj) {
423 //double drs[3];
424 const std::array<double,3> djj = mj.dipole_d_abs(sj);
425 double qmy15;
426 params >> qmy15;
427 minusSiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
428 const std::array<double,3> ejj = mj.dipole_e(sj);
429 PotForceDiQuadrupole(drs, dr2, ejj.data(), eii.data(), qmy15, f, m2, m1, u);
430
431 mi.Fquadrupolesub(si, f);
432 mj.Fdipoleadd(sj, f);
433 mi.Madd(m1);
434 mj.Madd(m2);
435 UpotXpoles += u;
436 for (unsigned short d = 0; d < 3; d++)
437 Virial[d] -= 0.5*drm[d] * f[d];
438 }
439 }
440 for (unsigned int si = 0; si < nd1; ++si) {
441 const std::array<double,3> dii = mi.dipole_d_abs(si);
442 const std::array<double,3> eii = mi.dipole_e(si);
443 // Dipole-Charge
444 for (unsigned sj = 0; sj < ne2; sj++) {
445 const std::array<double,3> djj = mj.charge_d_abs(sj);
446 double minusqmyper4pie0;
447 params >> minusqmyper4pie0;
448 minusSiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
449 PotForceChargeDipole(drs, dr2, eii.data(), minusqmyper4pie0, f, m1, u);
450
451 mi.Fdipolesub(si, f);
452 mj.Fchargeadd(sj, f);
453 mi.Madd(m1);
454
455 UpotXpoles += u;
456 for (unsigned short d = 0; d < 3; d++)
457 Virial[d] -= 0.5*drm[d] * f[d];
458 }
459 // Dipole-Quadrupole -----------------------
460 for (unsigned int sj = 0; sj < nq2; ++sj) {
461 //double drs[3];
462 const std::array<double,3> djj = mj.quadrupole_d_abs(sj);
463 double myq15;
464 params >> myq15;
465 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
466 const std::array<double,3> ejj = mj.quadrupole_e(sj);
467 PotForceDiQuadrupole(drs, dr2, eii.data(), ejj.data(), myq15, f, m1, m2, u);
468
469 mi.Fdipoleadd(si, f);
470 mj.Fquadrupolesub(sj, f);
471 mi.Madd(m1);
472 mj.Madd(m2);
473 UpotXpoles += u;
474 for (unsigned short d = 0; d < 3; ++d)
475 Virial[d] += 0.5*drm[d] * f[d];
476 }
477 // Dipole-Dipole ---------------------------
478 for (unsigned int sj = 0; sj < nd2; ++sj) {
479 const std::array<double,3> djj = mj.dipole_d_abs(sj);
480 double my2;
481 params >> my2;
482 double rffac;
483 params >> rffac;
484 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
485 const std::array<double,3> ejj = mj.dipole_e(sj);
486 PotForce2Dipole(drs, dr2, eii.data(), ejj.data(), my2, rffac, f, m1, m2, u, MyRF);
487
488 mi.Fdipoleadd(si, f);
489 mj.Fdipolesub(sj, f);
490 mi.Madd(m1);
491 mj.Madd(m2);
492 UpotXpoles += u;
493 for (unsigned short d = 0; d < 3; ++d)
494 Virial[d] += 0.5*drm[d] * f[d];
495 }
496 }
497
498 mi.Viadd(Virial);
499 mj.Viadd(Virial);
500
501 // check whether all parameters were used
502 mardyn_assert(params.eos());
503}
504
506inline void FluidPot(Molecule& mi, Molecule& mj, ParaStrm& params, double /*drm*/[3], double& Upot6LJ, double& UpotXpoles, double& MyRF, bool calculateLJ)
507{
508 double f[3];
509 double u;
510 double drs[3], dr2; // site distance vector & length^2
511 // no LJ interaction between equal solid atoms
512
513 const unsigned int nc1 = mi.numLJcenters();
514 const unsigned int nc2 = mj.numLJcenters();
515 for (unsigned int si = 0; si < nc1; ++si) {
516 const std::array<double,3> dii = mi.ljcenter_d_abs(si);
517 for (unsigned int sj = 0; sj < nc2; ++sj) {
518 const std::array<double,3> djj = mj.ljcenter_d_abs(sj);
519 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
520 double eps24;
521 params >> eps24;
522 double sig2;
523 params >> sig2;
524 double shift6;
525 params >> shift6; // must be 0.0 for full LJ
526
527 if (calculateLJ) {
528 PotForceLJ(drs, dr2, eps24, sig2, f, u);
529 u += shift6;
530 Upot6LJ += u;
531 }
532 }
533 }
534
535
536 double m1[3], m2[3]; // angular momenta
537
538 const unsigned ne1 = mi.numCharges();
539 const unsigned ne2 = mj.numCharges();
540 const unsigned int nq1 = mi.numQuadrupoles();
541 const unsigned int nq2 = mj.numQuadrupoles();
542 const unsigned int nd1 = mi.numDipoles();
543 const unsigned int nd2 = mj.numDipoles();
544 for (unsigned si = 0; si < ne1; si++) {
545 const std::array<double,3> dii = mi.charge_d_abs(si);
546 // Charge-Charge
547 for (unsigned sj = 0; sj < ne2; sj++) {
548 const std::array<double,3> djj = mj.charge_d_abs(sj);
549 double q1q2per4pie0; // 4pie0 = 1 in reduced units
550 params >> q1q2per4pie0;
551 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
552 PotForce2Charge(drs, dr2, q1q2per4pie0, f, u);
553 UpotXpoles += u;
554 }
555 // Charge-Quadrupole
556 for (unsigned sj = 0; sj < nq2; sj++) {
557 const std::array<double,3> djj = mj.quadrupole_d_abs(sj);
558 double qQ05per4pie0; // 4pie0 = 1 in reduced units
559 params >> qQ05per4pie0;
560 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
561 const std::array<double,3> ejj = mj.quadrupole_e(sj);
562 PotForceChargeQuadrupole(drs, dr2, ejj.data(), qQ05per4pie0, f, m2, u);
563 UpotXpoles += u;
564 }
565 // Charge-Dipole
566 for (unsigned sj = 0; sj < nd2; sj++) {
567 const std::array<double,3> djj = mj.dipole_d_abs(sj);
568 double minusqmyper4pie0;
569 params >> minusqmyper4pie0;
570 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
571 const std::array<double,3> ejj = mj.dipole_e(sj);
572 PotForceChargeDipole(drs, dr2, ejj.data(), minusqmyper4pie0, f, m2, u);
573 UpotXpoles += u;
574 }
575 }
576 for (unsigned int si = 0; si < nq1; ++si) {
577 const std::array<double,3> dii = mi.quadrupole_d_abs(si);
578 const std::array<double,3> eii = mi.quadrupole_e(si);
579
580 // Quadrupole-Charge
581 for (unsigned sj = 0; sj < ne2; sj++) {
582 const std::array<double,3> djj = mj.charge_d_abs(sj);
583 double qQ05per4pie0; // 4pie0 = 1 in reduced units
584 params >> qQ05per4pie0;
585 minusSiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
586 PotForceChargeQuadrupole(drs, dr2, eii.data(), qQ05per4pie0, f, m1, u);
587 UpotXpoles += u;
588 }
589 // Quadrupole-Quadrupole -------------------
590 for (unsigned int sj = 0; sj < nq2; ++sj) {
591 const std::array<double,3> djj = mj.quadrupole_d_abs(sj);
592 double q2075;
593 params >> q2075;
594 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
595 const std::array<double,3> ejj = mj.quadrupole_e(sj);
596 PotForce2Quadrupole(drs, dr2, eii.data(), ejj.data(), q2075, f, m1, m2, u);
597 UpotXpoles += u;
598 }
599 // Quadrupole-Dipole -----------------------
600 for (unsigned int sj = 0; sj < nd2; ++sj) {
601 const std::array<double,3> djj = mj.dipole_d_abs(sj);
602 double qmy15;
603 params >> qmy15;
604 minusSiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
605 const std::array<double,3> ejj = mj.dipole_e(sj);
606 //for(unsigned short d=0;d<3;++d) drs[d]=-drs[d]; // avoid that and toggle add/sub below
607 PotForceDiQuadrupole(drs, dr2, ejj.data(), eii.data(), qmy15, f, m2, m1, u);
608 UpotXpoles += u;
609 }
610 }
611 for (unsigned int si = 0; si < nd1; ++si) {
612 const std::array<double,3> dii = mi.dipole_d_abs(si);
613 const std::array<double,3> eii = mi.dipole_e(si);
614 // Dipole-Charge
615 for (unsigned sj = 0; sj < ne2; sj++) {
616 const std::array<double,3> djj = mj.charge_d_abs(sj);
617 double minusqmyper4pie0;
618 params >> minusqmyper4pie0;
619 minusSiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
620 PotForceChargeDipole(drs, dr2, eii.data(), minusqmyper4pie0, f, m1, u);
621 UpotXpoles += u;
622 }
623 // Dipole-Quadrupole -----------------------
624 for (unsigned int sj = 0; sj < nq2; ++sj) {
625 //double drs[3];
626 const std::array<double,3> djj = mj.quadrupole_d_abs(sj);
627 double myq15;
628 params >> myq15;
629 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
630 const std::array<double,3> ejj = mj.quadrupole_e(sj);
631 PotForceDiQuadrupole(drs, dr2, eii.data(), ejj.data(), myq15, f, m1, m2, u);
632 UpotXpoles += u;
633 }
634 // Dipole-Dipole ---------------------------
635 for (unsigned int sj = 0; sj < nd2; ++sj) {
636 //double drs[3];
637 const std::array<double,3> djj = mj.dipole_d_abs(sj);
638 double my2;
639 params >> my2;
640 double rffac;
641 params >> rffac;
642 SiteSiteDistanceAbs(dii.data(), djj.data(), drs, dr2);
643 const std::array<double,3> ejj = mj.dipole_e(sj);
644 PotForce2Dipole(drs, dr2, eii.data(), ejj.data(), my2, rffac, f, m1, m2, u, MyRF);
645 UpotXpoles += u;
646 }
647 }
648
649 // check whether all parameters were used
650 mardyn_assert(params.eos());
651}
652
653#endif /* POTFORCE_H_ */
654
FullMolecule
FullMolecule modeled as LJ sphere with point polarities.
Definition: FullMolecule.h:18
ParaStrm
Definition: ParaStrm.h:21
ParaStrm::eos
bool eos() const
end of stream reached?
Definition: ParaStrm.h:43
Generated on Thu Oct 21 2021 02:54:36 for ls1-MarDyn by doxygen 1.9.2