Skip to content

Commit 442e653

Browse files
momokchungmomokchung
committed
first implementation of expolinduce.cpp, rename invpolscale to polinv
1 parent 3a24887 commit 442e653

File tree

6 files changed

+266
-3
lines changed

6 files changed

+266
-3
lines changed

‎include/ff/hippomod.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,5 +74,6 @@ TINKER_EXTERN real* prepep;
7474
TINKER_EXTERN real* dmppep;
7575
TINKER_EXTERN int* lpep;
7676
TINKER_EXTERN real (*polscale)[3][3];
77+
TINKER_EXTERN real (*polinv)[3][3];
7778
TINKER_EXTERN ExpolScr scrtyp;
7879
}

‎src/acc/cmakesrc.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ hippo/edisp.cpp
4040
hippo/empole.cpp
4141
hippo/epolar.cpp
4242
hippo/erepel.cpp
43+
hippo/expolinduce.cpp
4344
hippo/field.cpp
4445
hippo/induce.cpp
4546
lflpiston.cpp

‎src/acc/hippo/expolinduce.cpp‎

Lines changed: 251 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,251 @@
1+
#include "ff/amoeba/induce.h"
2+
#include "ff/amoebamod.h"
3+
#include "ff/atom.h"
4+
#include "ff/hippo/induce.h"
5+
#include "ff/hippomod.h"
6+
#include "ff/switch.h"
7+
#include "tool/darray.h"
8+
#include "tool/ioprint.h"
9+
#include "tool/error.h"
10+
#include <tinker/detail/inform.hh>
11+
#include <tinker/detail/polpcg.hh>
12+
#include <tinker/detail/polpot.hh>
13+
#include <tinker/detail/units.hh>
14+
15+
namespace tinker {
16+
void induceMutualPcg4_acc(real (*uind)[3])
17+
{
18+
auto* field = work01_;
19+
auto* rsd = work02_;
20+
auto* zrsd = work03_;
21+
auto* conj = work04_;
22+
auto* vec = work05_;
23+
24+
bool dirguess = polpcg::pcgguess;
25+
// use sparse matrix preconditioner
26+
// or just use diagonal matrix preconditioner
27+
const bool sparse_prec = polpcg::pcgprec and (switchOff(Switch::USOLVE) > 0);
28+
29+
// zero out the induced dipoles at each site
30+
31+
bool predict = polpred != UPred::NONE;
32+
if (predict and nualt < maxualt) {
33+
predict = false;
34+
dirguess = true;
35+
}
36+
// get the electrostatic field due to permanent multipoles
37+
dfieldChgpen(field);
38+
// direct induced dipoles
39+
40+
#pragma acc parallel loop independent async\
41+
deviceptr(polarity,udir,field)
42+
for (int i = 0; i < n; ++i) {
43+
real poli = polarity[i];
44+
#pragma acc loop seq
45+
for (int j = 0; j < 3; ++j)
46+
udir[i][j] = poli * field[i][j];
47+
}
48+
49+
// call alterpol
50+
51+
// initial induced dipole
52+
if (predict) {
53+
ulspredSum(uind, nullptr);
54+
} else if (dirguess) {
55+
#pragma acc parallel loop independent async\
56+
deviceptr(polarity,uind,field,polinv)
57+
for (int i = 0; i < n; ++i) {
58+
real poli = polarity[i];
59+
#pragma acc loop seq
60+
for (int j = 0; j < 3; ++j) {
61+
uind[i][j] = poli * (polinv[i][j][1]*field[i][1]
62+
+ polinv[i][j][2]*field[i][2]
63+
+ polinv[i][j][3]*field[i][3]);
64+
}
65+
}
66+
} else {
67+
darray::zero(g::q0, n, uind);
68+
}
69+
// initial residual r(0)
70+
71+
// if do not use pcgguess, r(0) = E - T Zero = E
72+
// if use pcgguess, r(0) = E - (inv_alpha + Tu) alpha E
73+
// = E - E -Tu udir
74+
// = -Tu udir
75+
76+
if (predict) {
77+
ufieldChgpen(uind, field);
78+
#pragma acc parallel loop independent async\
79+
deviceptr(polarity_inv,rsd,udir,uind,field,polscale)
80+
for (int i = 0; i < n; ++i) {
81+
real pol = polarity_inv[i];
82+
#pragma acc loop seq
83+
for (int j = 0; j < 3; ++j)
84+
rsd[i][j] = (udir[i][j] - uind[i][1]*polscale[i][j][1]
85+
- uind[i][2]*polscale[i][j][2]
86+
- uind[i][3]*polscale[i][j][3])
87+
* pol + field[i][j];
88+
}
89+
} else if (dirguess) {
90+
ufieldChgpen(udir, rsd);
91+
} else {
92+
darray::copy(g::q0, n, rsd, field);
93+
}
94+
#pragma acc parallel loop independent async deviceptr(polarity,rsd)
95+
for (int i = 0; i < n; ++i) {
96+
if (polarity[i] == 0) {
97+
rsd[i][0] = 0;
98+
rsd[i][1] = 0;
99+
rsd[i][2] = 0;
100+
}
101+
}
102+
103+
// initial M r(0) and p(0)
104+
105+
if (sparse_prec) {
106+
sparsePrecondBuild2();
107+
sparsePrecondApply2(rsd, zrsd);
108+
} else {
109+
diagPrecond2(rsd, zrsd);
110+
}
111+
darray::copy(g::q0, n, conj, zrsd);
112+
113+
// initial r(0) M r(0)
114+
115+
real sum;
116+
sum = darray::dotThenReturn(g::q0, n, rsd, zrsd);
117+
118+
// conjugate gradient iteration of the mutual induced dipoles
119+
120+
const bool debug = inform::debug;
121+
const int politer = polpot::politer;
122+
const real poleps = polpot::poleps;
123+
const real debye = units::debye;
124+
const real pcgpeek = polpcg::pcgpeek;
125+
const int maxiter = 100; // see also subroutine induce0a in induce.f
126+
127+
bool done = false;
128+
int iter = 0;
129+
real eps = 100;
130+
real epsold;
131+
132+
while (not done) {
133+
++iter;
134+
135+
// vec = (inv_alpha + Tu) conj, field = -Tu conj
136+
// vec = inv_alpha * conj - field
137+
ufieldChgpen(conj, field);
138+
#pragma acc parallel loop independent async\
139+
deviceptr(polarity_inv,vec,conj,field,polscale)
140+
for (int i = 0; i < n; ++i) {
141+
real poli_inv = polarity_inv[i];
142+
#pragma acc loop seq
143+
for (int j = 0; j < 3; ++j)
144+
vec[i][j] = poli_inv * (conj[i][1]*polscale[i][j][1]
145+
+ conj[i][2]*polscale[i][j][2]
146+
+ conj[i][3]*polscale[i][j][3])
147+
- field[i][j];
148+
}
149+
150+
// a <- p T p
151+
real a;
152+
a = darray::dotThenReturn(g::q0, n, conj, vec);
153+
// a <- r M r / p T p
154+
if (a != 0)
155+
a = sum / a;
156+
157+
// u <- u + a p
158+
// r <- r - a T p
159+
#pragma acc parallel loop independent async\
160+
deviceptr(polarity,uind,conj,rsd,vec)
161+
for (int i = 0; i < n; ++i) {
162+
#pragma acc loop seq
163+
for (int j = 0; j < 3; ++j) {
164+
uind[i][j] += a * conj[i][j];
165+
rsd[i][j] -= a * vec[i][j];
166+
}
167+
if (polarity[i] == 0) {
168+
rsd[i][0] = 0;
169+
rsd[i][1] = 0;
170+
rsd[i][2] = 0;
171+
}
172+
}
173+
174+
// calculate/update M r
175+
if (sparse_prec)
176+
sparsePrecondApply2(rsd, zrsd);
177+
else
178+
diagPrecond2(rsd, zrsd);
179+
180+
real b;
181+
real sum1;
182+
sum1 = darray::dotThenReturn(g::q0, n, rsd, zrsd);
183+
b = sum1 / sum;
184+
if (sum == 0)
185+
b = 0;
186+
187+
// calculate/update p
188+
#pragma acc parallel loop independent async\
189+
deviceptr(conj,zrsd)
190+
for (int i = 0; i < n; ++i) {
191+
#pragma acc loop seq
192+
for (int j = 0; j < 3; ++j)
193+
conj[i][j] = zrsd[i][j] + b * conj[i][j];
194+
}
195+
196+
sum = sum1;
197+
198+
real epsd;
199+
epsd = darray::dotThenReturn(g::q0, n, rsd, rsd);
200+
epsold = eps;
201+
eps = epsd;
202+
eps = debye * REAL_SQRT(eps / n);
203+
204+
if (debug) {
205+
if (iter == 1) {
206+
print(stdout,
207+
"\n Determination of SCF Induced Dipole Moments\n\n"
208+
" Iter RMS Residual (Debye)\n\n");
209+
}
210+
print(stdout, " %8d %-16.10f\n", iter, eps);
211+
}
212+
213+
if (eps < poleps)
214+
done = true;
215+
// if (eps > epsold)
216+
// done = true;
217+
if (iter >= politer)
218+
done = true;
219+
220+
// apply a "peek" iteration to the mutual induced dipoles
221+
222+
if (done) {
223+
#pragma acc parallel loop independent async\
224+
deviceptr(polarity,uind,rsd)
225+
for (int i = 0; i < n; ++i) {
226+
real term = pcgpeek * polarity[i];
227+
#pragma acc loop seq
228+
for (int j = 0; j < 3; ++j)
229+
uind[i][j] += term * rsd[i][j];
230+
}
231+
}
232+
}
233+
234+
// print the results from the conjugate gradient iteration
235+
236+
if (debug) {
237+
print(stdout,
238+
" Induced Dipoles : Iterations %4d RMS "
239+
"Residual %14.10f\n",
240+
iter, eps);
241+
}
242+
243+
// terminate the calculation if dipoles failed to converge
244+
245+
// if (iter >= maxiter || eps > epsold) {
246+
if (iter >= maxiter) {
247+
printError();
248+
TINKER_THROW("INDUCE -- Warning, Induced Dipoles are not Converged");
249+
}
250+
}
251+
}

‎src/hippo/alterpol.cpp‎

Whitespace-only changes.

‎src/hippo/expol.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,14 +14,14 @@ void expolData(RcOp op)
1414

1515
if (op & RcOp::DEALLOC) {
1616
darray::deallocate(kpep, prepep, dmppep, lpep);
17-
darray::deallocate(polscale);
17+
darray::deallocate(polscale, polinv);
1818

1919
scrtyp = ExpolScr::NONE;
2020
}
2121

2222
if (op & RcOp::ALLOC) {
2323
darray::allocate(n, &kpep, &prepep, &dmppep, &lpep);
24-
darray::allocate(n, &polscale);
24+
darray::allocate(n, &polscale, &polinv);
2525
}
2626

2727
if (op & RcOp::INIT) {

‎src/hippo/induce.cpp‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,13 +39,23 @@ static void induceMutualPcg3(real (*uind)[3])
3939
TINKER_FCALL2(acc1, cu1, induceMutualPcg3, uind);
4040
}
4141

42+
// TODO add cuda implementation later
43+
void induceMutualPcg4_acc(real (*)[3]);
44+
static void induceMutualPcg4(real (*uind)[3])
45+
{
46+
induceMutualPcg4_acc(uind);
47+
}
48+
4249
void induce2(real (*ud)[3])
4350
{
4451
if (polpot::use_dirdamp) {
4552
induceMutualPcg3(ud);
4653
ulspredSave(ud, nullptr);
4754
} else {
48-
induceMutualPcg2(ud);
55+
if (polpot::use_expol)
56+
induceMutualPcg4(ud);
57+
else
58+
induceMutualPcg2(ud);
4959
ulspredSave(ud, nullptr);
5060
}
5161
inducePrint(ud);

0 commit comments

Comments
 (0)