13 #include "ComputeMsmSerialMgr.decl.h"
18 #include "ComputeMgr.decl.h"
20 #define MIN_DEBUG_LEVEL 3
65 CProxy_ComputeMsmSerialMgr msmProxy;
82 msmProxy(thisgroup), msmCompute(0), numSources(0), numArrived(0),
83 coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0),
84 msmsolver(0), msmcoord(0), msmforce(0)
86 CkpvAccess(BOCclass_group).computeMsmSerialMgr = thisgroup;
91 for (
int i=0; i < numSources; i++)
delete coordMsgs[i];
97 if (msmcoord)
delete[] msmcoord;
98 if (msmforce)
delete[] msmforce;
104 CProxy_ComputeMsmSerialMgr::ckLocalBranch(
105 CkpvAccess(BOCclass_group).computeMsmSerialMgr)->setCompute(
this);
118 if ( !
patchList[0].p->flags.doFullElectrostatics )
120 for (ap = ap.
begin(); ap != ap.
end(); ap++) {
122 Results *r = (*ap).forceBox->open();
123 (*ap).positionBox->close(&x);
124 (*ap).forceBox->close(&r);
131 int numLocalAtoms = 0;
132 for (ap = ap.
begin(); ap != ap.
end(); ap++) {
133 numLocalAtoms += (*ap).p->getNumAtoms();
143 for (ap = ap.
begin(); ap != ap.
end(); ap++) {
147 (*ap).positionBox->close(&x);
148 x = (*ap).avgPositionBox->open();
150 int numAtoms = (*ap).p->getNumAtoms();
152 for(
int i=0; i < numAtoms; i++)
156 data_ptr->
id = xExt[i].
id;
160 if (
patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
161 else { (*ap).positionBox->close(&x); }
164 CProxy_ComputeMsmSerialMgr msmProxy(
165 CkpvAccess(BOCclass_group).computeMsmSerialMgr);
166 msmProxy[0].recvCoord(msg);
174 double padding,
double gridspacing) {
175 const double NL_ZERO_TOLERANCE = 1e-4;
176 double xlen = 1, inv_xlen = 1;
177 double ylen = 1, inv_ylen = 1;
178 double zlen = 1, inv_zlen = 1;
182 double rpadx = padding * ulen_1;
183 double rpady = padding * vlen_1;
184 double rpadz = padding * wlen_1;
185 double rhx = gridspacing * ulen_1;
186 double rhy = gridspacing * vlen_1;
187 double rhz = gridspacing * wlen_1;
190 double xmin, xmax, ymin, ymax, zmin, zmax;
211 if (s.
x < xmin) xmin = s.
x;
212 else if (s.
x > xmax) xmax = s.
x;
213 if (s.
y < ymin) ymin = s.
y;
214 else if (s.
y > ymax) ymax = s.
y;
215 if (s.
z < zmin) zmin = s.
z;
216 else if (s.
z > zmax) zmax = s.
z;
219 printf(
"*** xmin=%.4f xmax=%.4f\n", xmin, xmax);
220 printf(
"*** ymin=%.4f ymax=%.4f\n", ymin, ymax);
221 printf(
"*** zmin=%.4f zmax=%.4f\n", zmin, zmax);
224 if ( ! is_periodic_x) {
228 double mupper = ceil(xmax / (2*rhx));
229 double mlower = floor(xmin / (2*rhx));
233 rc.
x = 0.5*(xmin + xmax);
235 if (xlen < NL_ZERO_TOLERANCE) {
243 if ( ! is_periodic_y) {
247 double mupper = ceil(ymax / (2*rhy));
248 double mlower = floor(ymin / (2*rhy));
252 rc.
y = 0.5*(ymin + ymax);
254 if (ylen < NL_ZERO_TOLERANCE) {
262 if ( ! is_periodic_z) {
266 double mupper = ceil(zmax / (2*rhz));
267 double mlower = floor(zmin / (2*rhz));
271 rc.
z = 0.5*(zmin + zmax);
273 if (zlen < NL_ZERO_TOLERANCE) {
282 printf(
"xmin=%g xmax=%g\n", xmin, xmax);
283 printf(
"ymin=%g ymax=%g\n", ymin, ymax);
284 printf(
"zmin=%g zmax=%g\n", zmin, zmax);
285 printf(
"xlen=%g ylen=%g zlen=%g\n", xlen, ylen, zlen);
286 printf(
"rc_x=%g rc_y=%g rc_z=%g\n", rc.
x, rc.
y, rc.
z);
290 c.
x = (u.
x*rc.
x + v.
x*rc.
y + w.
x*rc.
z) + c.
x;
291 c.
y = (u.
y*rc.
x + v.
y*rc.
y + w.
y*rc.
z) + c.
y;
292 c.
z = (u.
z*rc.
x + v.
z*rc.
y + w.
z*rc.
z) + c.
z;
295 printf(
"c_x=%g c_y=%g c_z=%g\n", c.
x, c.
y, c.
z);
309 if ( ! numSources ) {
312 for (
int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
320 for ( i=0; i < msg->
numAtoms; ++i ) {
324 coordMsgs[numArrived] = msg;
327 if ( numArrived < numSources )
return;
344 if ( ! msmsolver )
NAMD_die(
"unable to create MSM solver");
346 double cutoff = simParams->
cutoff;
360 if (rc)
NAMD_die(
"unable to configure MSM solver");
361 Vector u=lattice.
a(), v=lattice.
b(), w=lattice.
c(), c=lattice.
origin();
368 isperiodic, numAtoms, coord, padding, gridspacing);
370 double vec1[3], vec2[3], vec3[3], center[3];
384 printf(
"dielectric = %g\n", dielectric);
385 printf(
"vec1 = %g %g %g\n", vec1[0], vec1[1], vec1[2]);
386 printf(
"vec2 = %g %g %g\n", vec2[0], vec2[1], vec2[2]);
387 printf(
"vec3 = %g %g %g\n", vec3[0], vec3[1], vec3[2]);
388 printf(
"center = %g %g %g\n", center[0], center[1], center[2]);
389 printf(
"cutoff = %g\n", cutoff);
390 printf(
"numatoms = %d\n", numAtoms);
392 rc =
NL_msm_setup(msmsolver, cutoff, vec1, vec2, vec3, center, msmflags);
393 if (rc)
NAMD_die(
"unable to set up MSM solver");
394 msmcoord =
new double[4*numAtoms];
395 msmforce =
new double[3*numAtoms];
396 if (msmcoord==0 || msmforce==0)
NAMD_die(
"can't allocate MSM atom buffers");
399 for (i = 0; i < numAtoms; i++) {
400 msmcoord[4*i+3] = celec * coord[i].
charge;
405 for (i = 0; i < numAtoms; i++) {
410 for (i = 0; i < numAtoms; i++) {
416 if (rc)
NAMD_die(
"error evaluating MSM forces");
417 for (i = 0; i < numAtoms; i++) {
418 force[i].
x = msmforce[3*i ];
419 force[i].
y = msmforce[3*i+1];
420 force[i].
z = msmforce[3*i+2];
423 for (
int k=0; k < 3; k++) {
424 for (
int l=0; l < 3; l++) {
430 for (
int j=0; j < numSources; j++) {
435 for (
int i=0; i < cmsg->
numAtoms; i++) {
441 for (
int k=0; k < 3; k++) {
442 for (
int l=0; l < 3; l++) {
443 fmsg->
virial[k][l] = virial[k][l];
449 for (
int k=0; k < 3; k++) {
450 for (
int l=0; l < 3; l++) {
474 for (ap = ap.
begin(); ap != ap.
end(); ap++) {
475 Results *r = (*ap).forceBox->open();
477 int numAtoms = (*ap).p->getNumAtoms();
479 for(
int i=0; i<numAtoms; ++i) {
480 f[i].
x += results_ptr->
x;
481 f[i].
y += results_ptr->
y;
482 f[i].
z += results_ptr->
z;
486 (*ap).forceBox->close(&r);
490 reduction->
item(REDUCTION_VIRIAL_SLOW_XX) += msg->
virial[0][0];
491 reduction->
item(REDUCTION_VIRIAL_SLOW_XY) += msg->
virial[0][1];
492 reduction->
item(REDUCTION_VIRIAL_SLOW_XZ) += msg->
virial[0][2];
493 reduction->
item(REDUCTION_VIRIAL_SLOW_YX) += msg->
virial[1][0];
494 reduction->
item(REDUCTION_VIRIAL_SLOW_YY) += msg->
virial[1][1];
495 reduction->
item(REDUCTION_VIRIAL_SLOW_YZ) += msg->
virial[1][2];
496 reduction->
item(REDUCTION_VIRIAL_SLOW_ZX) += msg->
virial[2][0];
497 reduction->
item(REDUCTION_VIRIAL_SLOW_ZY) += msg->
virial[2][1];
498 reduction->
item(REDUCTION_VIRIAL_SLOW_ZZ) += msg->
virial[2][2];
502 #include "ComputeMsmSerialMgr.def.h"
void setCompute(ComputeMsmSerial *c)
static PatchMap * Object()
SimParameters * simParameters
ComputeHomePatchList patchList
ComputeMsmSerialAtom * coord
NL_Msm * NL_msm_create(void)
int NL_msm_configure(NL_Msm *pm, double gridspacing, int approx, int split, int nlevels)
SubmitReduction * willSubmit(int setID, int size=-1)
void recvForce(MsmSerialForceMsg *)
static ReductionMgr * Object(void)
ResizeArrayIter< T > end(void) const
void saveResults(MsmSerialForceMsg *)
int NL_msm_compute_force(NL_Msm *pm, double *felec, double *uelec, const double *atom, int natoms)
virtual ~ComputeMsmSerial()
void NAMD_die(const char *err_msg)
std::vector< std::string > split(const std::string &text, std::string delimiter)
void recvCoord(MsmSerialCoordMsg *)
ComputeMsmSerial(ComputeID c)
static void rescale_nonperiodic_cell(Vector &u, Vector &v, Vector &w, Vector &c, Vector &ru, Vector &rv, Vector &rw, int isperiodic, int numatoms, const ComputeMsmSerialAtom *coord, double padding, double gridspacing)
void NL_msm_destroy(NL_Msm *pm)
int NL_msm_setup(NL_Msm *msm, double cutoff, double cellvec1[3], double cellvec2[3], double cellvec3[3], double cellcenter[3], int msmflags)
ResizeArrayIter< T > begin(void) const