Skip to content

Commit 3784a1c

Browse files
zhi-wangzhi-wang
committed
expol: add screening S2 and G
1 parent eb0c62c commit 3784a1c

File tree

6 files changed

+204
-23
lines changed

6 files changed

+204
-23
lines changed

‎include/seq/pair_alterpol.h‎

Lines changed: 46 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#pragma once
22
#include "ff/hippo/expolscr.h"
3+
#include "math/sinhc.h"
34
#include "math/switch.h"
45
#include "seq/seq.h"
56

@@ -13,9 +14,8 @@ inline void damp_expl(ExpolScr scrtyp, real& restrict s2, real& restrict ds2, re
1314
constexpr real inv2 = 1. / 2, inv3 = 1. / 3;
1415
constexpr real one = 1.;
1516

16-
real alphaik, dmpik2, dampik, dampik2, expik, s;
17-
1817
if (scrtyp == ExpolScr::S2U) {
18+
real alphaik, dmpik2, dampik, dampik2, expik, s;
1919
alphaik = REAL_SQRT(alphai * alphak);
2020
dmpik2 = inv2 * alphaik;
2121
dampik = dmpik2 * r;
@@ -24,6 +24,42 @@ inline void damp_expl(ExpolScr scrtyp, real& restrict s2, real& restrict ds2, re
2424
s = (one + dampik + dampik2 * inv3) * expik;
2525
s2 = s * s;
2626
if (DO_G) ds2 = s * (-alphaik * inv3) * (dampik + dampik2) * expik;
27+
} else if (scrtyp == ExpolScr::S2) {
28+
real pfac = 2 / (alphai + alphak);
29+
real r2 = r * r;
30+
pfac = pfac * pfac;
31+
pfac = pfac * alphai * alphak;
32+
pfac = pfac * pfac * pfac;
33+
pfac *= r2;
34+
35+
real a = alphai * r / 2, b = alphak * r / 2;
36+
real c = (a + b) / 2, d = (b - a) / 2;
37+
real expmc = REAL_EXP(-c);
38+
39+
real c2 = c * c;
40+
real d2 = d * d;
41+
real c2d2 = (c * d) * (c * d);
42+
real f1d, f2d, f3d;
43+
fsinhc3(d, f1d, f2d, f3d);
44+
45+
real s;
46+
s = f1d * (c + 1) + f2d * c2;
47+
s /= r;
48+
s *= expmc;
49+
s2 = pfac * s * s;
50+
51+
if (DO_G) {
52+
real ds;
53+
ds = f1d * c2 + f2d * ((c - 2) * c2 - (c + 1) * d2) - f3d * c2d2;
54+
ds /= -r2;
55+
ds *= expmc;
56+
ds2 = pfac * 2 * s * ds;
57+
}
58+
59+
} else if (scrtyp == ExpolScr::G) {
60+
real alphaik = REAL_SQRT(alphai * alphak);
61+
s2 = REAL_EXP(-alphaik / (real)10 * r * r);
62+
if (DO_G) ds2 = (-alphaik / (real)5) * r * s2;
2763
}
2864

2965
s2 = sizik * s2;
@@ -92,7 +128,7 @@ inline void pair_dexpol(ExpolScr scrtyp, real r, real pscale, real cut, real off
92128
real ds2k = springk * ds2 * pscale;
93129

94130
// compute rotation matrix
95-
real ai[3][3], ak[3][3];
131+
real ai[3][3];
96132
ai[0][2] = xr / r;
97133
ai[1][2] = yr / r;
98134
ai[2][2] = zr / r;
@@ -116,24 +152,16 @@ inline void pair_dexpol(ExpolScr scrtyp, real r, real pscale, real cut, real off
116152
ai[0][1] = ai[2][0] * ai[1][2] - ai[1][0] * ai[2][2];
117153
ai[1][1] = ai[0][0] * ai[2][2] - ai[2][0] * ai[0][2];
118154
ai[2][1] = ai[1][0] * ai[0][2] - ai[0][0] * ai[1][2];
119-
ak[0][0] = ai[0][0];
120-
ak[1][0] = ai[1][0];
121-
ak[2][0] = ai[2][0];
122-
ak[0][1] = -ai[0][1];
123-
ak[1][1] = -ai[1][1];
124-
ak[2][1] = -ai[2][1];
125-
ak[0][2] = -ai[0][2];
126-
ak[1][2] = -ai[1][2];
127-
ak[2][2] = -ai[2][2];
155+
// ak[][0] = ai[][0], ak[][1] = -ai[][1], ak[][2] = -ai[][2]
128156

129157
// local frame force
130158
real frcil[3], frckl[3];
131159
real uixl = uix * ai[0][0] + uiy * ai[1][0] + uiz * ai[2][0];
132160
real uiyl = uix * ai[0][1] + uiy * ai[1][1] + uiz * ai[2][1];
133161
real uizl = uix * ai[0][2] + uiy * ai[1][2] + uiz * ai[2][2];
134-
real ukxl = -(ukx * ak[0][0] + uky * ak[1][0] + ukz * ak[2][0]);
135-
real ukyl = -(ukx * ak[0][1] + uky * ak[1][1] + ukz * ak[2][1]);
136-
real ukzl = -(ukx * ak[0][2] + uky * ak[1][2] + ukz * ak[2][2]);
162+
real ukxl = -(ukx * ai[0][0] + uky * ai[1][0] + ukz * ai[2][0]);
163+
real ukyl = ukx * ai[0][1] + uky * ai[1][1] + ukz * ai[2][1];
164+
real ukzl = ukx * ai[0][2] + uky * ai[1][2] + ukz * ai[2][2];
137165
frcil[2] = uizl * uizl * ds2i;
138166
frckl[2] = ukzl * ukzl * ds2k;
139167
// local frame torque
@@ -151,9 +179,9 @@ inline void pair_dexpol(ExpolScr scrtyp, real r, real pscale, real cut, real off
151179
real frcxi = ai[0][0] * frcil[0] + ai[0][1] * frcil[1] + ai[0][2] * frcil[2];
152180
real frcyi = ai[1][0] * frcil[0] + ai[1][1] * frcil[1] + ai[1][2] * frcil[2];
153181
real frczi = ai[2][0] * frcil[0] + ai[2][1] * frcil[1] + ai[2][2] * frcil[2];
154-
real frcxk = ak[0][0] * frckl[0] + ak[0][1] * frckl[1] + ak[0][2] * frckl[2];
155-
real frcyk = ak[1][0] * frckl[0] + ak[1][1] * frckl[1] + ak[1][2] * frckl[2];
156-
real frczk = ak[2][0] * frckl[0] + ak[2][1] * frckl[1] + ak[2][2] * frckl[2];
182+
real frcxk = ai[0][0] * frckl[0] - ai[0][1] * frckl[1] - ai[0][2] * frckl[2];
183+
real frcyk = ai[1][0] * frckl[0] - ai[1][1] * frckl[1] - ai[1][2] * frckl[2];
184+
real frczk = ai[2][0] * frckl[0] - ai[2][1] * frckl[1] - ai[2][2] * frckl[2];
157185
frc[0] = f * (frcxk - frcxi);
158186
frc[1] = f * (frcyk - frcyi);
159187
frc[2] = f * (frczk - frczi);

‎test/expol.cpp‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -228,3 +228,121 @@ TEST_CASE("ExchangePolarization-4", "[ff][hippo][expol]")
228228
finish();
229229
testEnd();
230230
}
231+
232+
TEST_CASE("ExchangePolarization-5", "[ff][hippo][expol]")
233+
{
234+
const char* xn = "crys.xyz";
235+
const char* kn = "expols2.key";
236+
const char* argv[] = {"dummy", xn, "-k", kn};
237+
int argc = 4;
238+
239+
TestFile fx1(TINKER9_DIRSTR "/test/file/expol/crys.xyz");
240+
TestFile fk1(TINKER9_DIRSTR "/test/file/expol/expol.key", kn,
241+
"\n"
242+
"exchange-polar s2"
243+
"\n");
244+
TestFile fp1(TINKER9_DIRSTR "/test/file/expol/expol.prm");
245+
246+
const double eps_e = testGetEps(0.0001, 0.0001);
247+
const double eps_g = testGetEps(0.0001, 0.0001);
248+
const double eps_v = testGetEps(0.001, 0.001);
249+
250+
TestReference r(TINKER9_DIRSTR "/test/ref/expol.5.txt");
251+
auto ref_c = r.getCount();
252+
auto ref_e = r.getEnergy();
253+
auto ref_v = r.getVirial();
254+
auto ref_g = r.getGradient();
255+
256+
rc_flag = calc::xyz | calc::vmask;
257+
testBeginWithArgs(argc, argv);
258+
initialize();
259+
260+
energy(calc::v0);
261+
COMPARE_REALS(esum, ref_e, eps_e);
262+
263+
energy(calc::v1);
264+
COMPARE_REALS(esum, ref_e, eps_e);
265+
COMPARE_GRADIENT(ref_g, eps_g);
266+
for (int i = 0; i < 3; ++i)
267+
for (int j = 0; j < 3; ++j)
268+
COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v);
269+
270+
energy(calc::v3);
271+
COMPARE_REALS(esum, ref_e, eps_e);
272+
COMPARE_INTS(countReduce(nep), ref_c);
273+
274+
energy(calc::v4);
275+
COMPARE_REALS(esum, ref_e, eps_e);
276+
COMPARE_GRADIENT(ref_g, eps_g);
277+
278+
energy(calc::v5);
279+
COMPARE_GRADIENT(ref_g, eps_g);
280+
281+
energy(calc::v6);
282+
COMPARE_GRADIENT(ref_g, eps_g);
283+
for (int i = 0; i < 3; ++i)
284+
for (int j = 0; j < 3; ++j)
285+
COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v);
286+
287+
finish();
288+
testEnd();
289+
}
290+
291+
TEST_CASE("ExchangePolarization-6", "[ff][hippo][expol]")
292+
{
293+
const char* xn = "crys.xyz";
294+
const char* kn = "expolg.key";
295+
const char* argv[] = {"dummy", xn, "-k", kn};
296+
int argc = 4;
297+
298+
TestFile fx1(TINKER9_DIRSTR "/test/file/expol/crys.xyz");
299+
TestFile fk1(TINKER9_DIRSTR "/test/file/expol/expol.key", kn,
300+
"\n"
301+
"exchange-polar g"
302+
"\n");
303+
TestFile fp1(TINKER9_DIRSTR "/test/file/expol/expol.prm");
304+
305+
const double eps_e = testGetEps(0.0001, 0.0001);
306+
const double eps_g = testGetEps(0.0001, 0.0001);
307+
const double eps_v = testGetEps(0.001, 0.001);
308+
309+
TestReference r(TINKER9_DIRSTR "/test/ref/expol.6.txt");
310+
auto ref_c = r.getCount();
311+
auto ref_e = r.getEnergy();
312+
auto ref_v = r.getVirial();
313+
auto ref_g = r.getGradient();
314+
315+
rc_flag = calc::xyz | calc::vmask;
316+
testBeginWithArgs(argc, argv);
317+
initialize();
318+
319+
energy(calc::v0);
320+
COMPARE_REALS(esum, ref_e, eps_e);
321+
322+
energy(calc::v1);
323+
COMPARE_REALS(esum, ref_e, eps_e);
324+
COMPARE_GRADIENT(ref_g, eps_g);
325+
for (int i = 0; i < 3; ++i)
326+
for (int j = 0; j < 3; ++j)
327+
COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v);
328+
329+
energy(calc::v3);
330+
COMPARE_REALS(esum, ref_e, eps_e);
331+
COMPARE_INTS(countReduce(nep), ref_c);
332+
333+
energy(calc::v4);
334+
COMPARE_REALS(esum, ref_e, eps_e);
335+
COMPARE_GRADIENT(ref_g, eps_g);
336+
337+
energy(calc::v5);
338+
COMPARE_GRADIENT(ref_g, eps_g);
339+
340+
energy(calc::v6);
341+
COMPARE_GRADIENT(ref_g, eps_g);
342+
for (int i = 0; i < 3; ++i)
343+
for (int j = 0; j < 3; ++j)
344+
COMPARE_REALS(vir[i * 3 + j], ref_v[i][j], eps_v);
345+
346+
finish();
347+
testEnd();
348+
}

‎test/file/expol/crys.xyz‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
4 Na+_Cl-_Na+_Cl-
2+
100.000000 100.000000 100.000000 90.000000 90.000000 90.000000
3+
1 Na+ 0.000000 0.000000 0.000000 7
4+
2 Cl- 0.000000 1.000000 2.270000 15
5+
3 Na+ 1.500000 -1.500000 1.500000 7
6+
4 Cl- 2.270000 0.000000 1.000000 15

‎test/file/expol/expol.key‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,4 @@ chgpen 21 12.760026532535434 5.683167591839673
1717
polarize 1 1.21889959 2
1818
polarize 2 0.14304110 1
1919

20-
screening s2u
21-
digits 10
22-
23-
polar-eps 1.0E-14
24-
2520
polarizeterm only

‎test/ref/expol.5.txt‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
expol.cpp -- Screening S2
2+
3+
Energy Component Breakdown : Kcal/mole Interactions
4+
Polarization -88.59818254 6
5+
6+
Internal Virial Tensor : 76.339 50.971 -14.081
7+
50.971 125.437 -61.029
8+
-14.081 -61.029 69.980
9+
10+
Cartesian Gradient Breakdown over Individual Atoms :
11+
12+
Type Atom dE/dX dE/dY dE/dZ Norm
13+
14+
Anlyt 1 -13.78553573 11.67733212 -25.81538711 31.50928899
15+
Anlyt 2 -6.78341068 -4.70227541 14.85762924 16.99632905
16+
Anlyt 3 -38.50314991 -86.75934147 50.59101447 107.55987463
17+
Anlyt 4 59.07209632 79.78428476 -39.63325661 106.89171945

‎test/ref/expol.6.txt‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
expol.cpp -- Screening G
2+
3+
Energy Component Breakdown : Kcal/mole Interactions
4+
Polarization -24.89981598 6
5+
6+
Internal Virial Tensor : 22.065 -1.907 -7.544
7+
-1.907 18.320 -5.911
8+
-7.544 -5.911 9.928
9+
10+
Cartesian Gradient Breakdown over Individual Atoms :
11+
12+
Type Atom dE/dX dE/dY dE/dZ Norm
13+
14+
Anlyt 1 -2.20768138 2.73909441 -2.32249154 4.21550261
15+
Anlyt 2 -6.47893686 1.83005418 3.51447391 7.59455384
16+
Anlyt 3 -3.04773596 -10.99303159 6.28394314 13.02395406
17+
Anlyt 4 11.73435421 6.42388300 -7.47592551 15.32484270

0 commit comments

Comments
 (0)