13 myGrid(grid), N(natoms) {
15 M =
new double[3*N*
order];
16 dM =
new double[3*N*
order];
26 void LjPmeRealSpace::fill_b_spline(
double *p) {
35 fr[0] = p[index]; fr[0] -= (int)(fr[0]);
36 fr[1] = p[index+1]; fr[1] -= (int)(fr[1]);
37 fr[2] = p[index+2]; fr[2] -= (int)(fr[2]);
38 compute_LjPme_b_spline<order>(fr, Mi, dMi);
46 switch (myGrid.
order) {
48 fill_charges_order<4>(q_arr, f_arr, fz_arr, p);
51 fill_charges_order<6>(q_arr, f_arr, fz_arr, p);
54 fill_charges_order<8>(q_arr, f_arr, fz_arr, p);
57 fill_charges_order<10>(q_arr, f_arr, fz_arr, p);
59 default:
NAMD_die(
"unsupported LJPMEInterpOrder");
64 void LjPmeRealSpace::fill_charges_order(
float **q_arr,
char *f_arr,
char *fz_arr,
68 int K1, K2, K3, dim2, dim3;
71 K1=myGrid.
K1; K2=myGrid.
K2; K3=myGrid.
K3; dim2=myGrid.
dim2; dim3=myGrid.
dim3;
74 this->fill_b_spline<order>(p);
81 u2i = (int)(p[index+1]);
82 u3i = (int)(p[index+2]);
88 for (j=0; j<
order; j++) {
93 ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2;
95 for (k=0; k<
order; k++) {
98 m1m2 = m1*Mi[
order+k];
100 ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));
101 float *qline = q_arr[ind2];
103 q_arr[ind2] = qline =
new float[dim3];
104 memset( (
void*) qline, 0, dim3 *
sizeof(
float) );
107 for (l=0; l<
order; l++) {
110 m3 = Mi[2*
order + l];
112 ind = u3 + (u3 < 0 ? K3 : 0);
113 qline[ind] += m1m2*m3;
118 for (l=0; l<
order; l++) {
120 int ind = u3 + (u3 < 0 ? K3 : 0);
127 double *force,
const Lattice &lattice) {
128 switch (myGrid.
order) {
130 compute_scaledForces_order<4>(q_arr, pos, force, lattice);
133 compute_scaledForces_order<6>(q_arr, pos, force, lattice);
136 compute_scaledForces_order<8>(q_arr, pos, force, lattice);
139 compute_scaledForces_order<10>(q_arr, pos, force, lattice);
141 default:
NAMD_die(
"unsupported LJPMEInterpOrder");
146 void LjPmeRealSpace::compute_scaledForces_order(
const float *
const *q_arr,
double *pos,
147 double *force,
const Lattice &lattice) {
148 int i, j, k, l, stride;
151 int K1, K2, K3, dim2;
156 double r1x = recip1.
x;
157 double r1y = recip1.
y;
158 double r1z = recip1.
z;
159 double r2x = recip2.
x;
160 double r2y = recip2.
y;
161 double r2z = recip2.
z;
162 double r3x = recip3.
x;
163 double r3y = recip3.
y;
164 double r3z = recip3.
z;
166 K1=myGrid.
K1; K2=myGrid.
K2; K3=myGrid.
K3; dim2=myGrid.
dim2;
170 for (i=0; i<N; i++) {
172 int u1, u2, u2i, u3i;
174 u1 = (int)(pos[4*i]);
175 u2i = (int)(pos[4*i+1]);
176 u3i = (int)(pos[4*i+2]);
182 for (j=0; j<
order; j++) {
188 ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2;
190 for (k=0; k<
order; k++) {
191 double m2, d2, m1m2, m1d2, d1m2;
199 ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));
200 const float *qline = q_arr[ind2];
201 for (l=0; l<
order; l++) {
205 d3=K3*dMi[2*
order+l];
207 ind = u3 + (u3 < 0 ? K3 : 0);
209 f1 -= d1m2 * m3 * term;
210 f2 -= m1d2 * m3 * term;
211 f3 -= m1m2 * d3 * term;
218 force[3*i] += f1*r1x + f2*r2x + f3*r3x;
219 force[3*i+1] += f1*r1y + f2*r2y + f3*r3y;
220 force[3*i+2] += f1*r1z + f2*r2z + f3*r3z;
void compute_scaledForces(const float *const *q_arr, double *pos, double *force, const Lattice &lattice)
NAMD_HOST_DEVICE Vector a_r() const
NAMD_HOST_DEVICE Vector b_r() const
void NAMD_die(const char *err_msg)
NAMD_HOST_DEVICE Vector c_r() const
LjPmeRealSpace(LjPmeGrid grid, int natoms)
void fill_charges(float **q_arr, char *f_arr, char *fz_arr, double *p)