00001
00007 #include "InfoStream.h"
00008 #include "Node.h"
00009 #include "PatchMap.h"
00010 #include "PatchMap.inl"
00011 #include "AtomMap.h"
00012 #include "ComputeMsmSerial.h"
00013 #include "ComputeMsmSerialMgr.decl.h"
00014 #include "PatchMgr.h"
00015 #include "Molecule.h"
00016 #include "ReductionMgr.h"
00017 #include "ComputeMgr.h"
00018 #include "ComputeMgr.decl.h"
00019
00020 #define MIN_DEBUG_LEVEL 3
00021 #include "Debug.h"
00022 #include "SimParameters.h"
00023 #include "WorkDistrib.h"
00024 #include "varsizemsg.h"
00025 #include "msm.h"
00026 #include <stdlib.h>
00027 #include <stdio.h>
00028 #include <errno.h>
00029
00030
00031 struct ComputeMsmSerialAtom {
00032 Position position;
00033 float charge;
00034 int id;
00035 };
00036
00037 typedef Force MsmSerialForce;
00038
00039 class MsmSerialCoordMsg : public CMessage_MsmSerialCoordMsg {
00040 public:
00041 int sourceNode;
00042 int numAtoms;
00043 Lattice lattice;
00044 ComputeMsmSerialAtom *coord;
00045 };
00046
00047 class MsmSerialForceMsg : public CMessage_MsmSerialForceMsg {
00048 public:
00049 BigReal energy;
00050 BigReal virial[3][3];
00051 MsmSerialForce *force;
00052 };
00053
00054 class ComputeMsmSerialMgr : public CBase_ComputeMsmSerialMgr {
00055 public:
00056 ComputeMsmSerialMgr();
00057 ~ComputeMsmSerialMgr();
00058
00059 void setCompute(ComputeMsmSerial *c) { msmCompute = c; }
00060
00061 void recvCoord(MsmSerialCoordMsg *);
00062 void recvForce(MsmSerialForceMsg *);
00063
00064 private:
00065 CProxy_ComputeMsmSerialMgr msmProxy;
00066 ComputeMsmSerial *msmCompute;
00067
00068 int numSources;
00069 int numArrived;
00070 MsmSerialCoordMsg **coordMsgs;
00071 int numAtoms;
00072 ComputeMsmSerialAtom *coord;
00073 MsmSerialForce *force;
00074 MsmSerialForceMsg *oldmsg;
00075
00076 NL_Msm *msmsolver;
00077 double *msmcoord;
00078 double *msmforce;
00079 };
00080
00081 ComputeMsmSerialMgr::ComputeMsmSerialMgr() :
00082 msmProxy(thisgroup), msmCompute(0), numSources(0), numArrived(0),
00083 coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0),
00084 msmsolver(0), msmcoord(0), msmforce(0)
00085 {
00086 CkpvAccess(BOCclass_group).computeMsmSerialMgr = thisgroup;
00087 }
00088
00089 ComputeMsmSerialMgr::~ComputeMsmSerialMgr()
00090 {
00091 for (int i=0; i < numSources; i++) delete coordMsgs[i];
00092 delete [] coordMsgs;
00093 delete [] coord;
00094 delete [] force;
00095 delete oldmsg;
00096 if (msmsolver) NL_msm_destroy(msmsolver);
00097 if (msmcoord) delete[] msmcoord;
00098 if (msmforce) delete[] msmforce;
00099 }
00100
00101 ComputeMsmSerial::ComputeMsmSerial(ComputeID c) :
00102 ComputeHomePatches(c)
00103 {
00104 CProxy_ComputeMsmSerialMgr::ckLocalBranch(
00105 CkpvAccess(BOCclass_group).computeMsmSerialMgr)->setCompute(this);
00106 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00107 }
00108
00109 ComputeMsmSerial::~ComputeMsmSerial()
00110 {
00111 }
00112
00113 void ComputeMsmSerial::doWork()
00114 {
00115 ResizeArrayIter<PatchElem> ap(patchList);
00116
00117
00118 if ( ! patchList[0].p->flags.doFullElectrostatics )
00119 {
00120 for (ap = ap.begin(); ap != ap.end(); ap++) {
00121 CompAtom *x = (*ap).positionBox->open();
00122 Results *r = (*ap).forceBox->open();
00123 (*ap).positionBox->close(&x);
00124 (*ap).forceBox->close(&r);
00125 }
00126 reduction->submit();
00127 return;
00128 }
00129
00130
00131 int numLocalAtoms = 0;
00132 for (ap = ap.begin(); ap != ap.end(); ap++) {
00133 numLocalAtoms += (*ap).p->getNumAtoms();
00134 }
00135
00136 MsmSerialCoordMsg *msg = new (numLocalAtoms, 0) MsmSerialCoordMsg;
00137 msg->sourceNode = CkMyPe();
00138 msg->numAtoms = numLocalAtoms;
00139 msg->lattice = patchList[0].p->flags.lattice;
00140 ComputeMsmSerialAtom *data_ptr = msg->coord;
00141
00142
00143 for (ap = ap.begin(); ap != ap.end(); ap++) {
00144 CompAtom *x = (*ap).positionBox->open();
00145 CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
00146 if ( patchList[0].p->flags.doMolly ) {
00147 (*ap).positionBox->close(&x);
00148 x = (*ap).avgPositionBox->open();
00149 }
00150 int numAtoms = (*ap).p->getNumAtoms();
00151
00152 for(int i=0; i < numAtoms; i++)
00153 {
00154 data_ptr->position = x[i].position;
00155 data_ptr->charge = x[i].charge;
00156 data_ptr->id = xExt[i].id;
00157 ++data_ptr;
00158 }
00159
00160 if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00161 else { (*ap).positionBox->close(&x); }
00162 }
00163
00164 CProxy_ComputeMsmSerialMgr msmProxy(
00165 CkpvAccess(BOCclass_group).computeMsmSerialMgr);
00166 msmProxy[0].recvCoord(msg);
00167 }
00168
00169
00170 static void rescale_nonperiodic_cell(
00171 Vector &u, Vector &v, Vector &w, Vector &c,
00172 Vector &ru, Vector &rv, Vector &rw,
00173 int isperiodic, int numatoms, const ComputeMsmSerialAtom *coord,
00174 double padding, double gridspacing) {
00175 const double NL_ZERO_TOLERANCE = 1e-4;
00176 double xlen = 1, inv_xlen = 1;
00177 double ylen = 1, inv_ylen = 1;
00178 double zlen = 1, inv_zlen = 1;
00179 double ulen_1 = u.rlength();
00180 double vlen_1 = v.rlength();
00181 double wlen_1 = w.rlength();
00182 double rpadx = padding * ulen_1;
00183 double rpady = padding * vlen_1;
00184 double rpadz = padding * wlen_1;
00185 double rhx = gridspacing * ulen_1;
00186 double rhy = gridspacing * vlen_1;
00187 double rhz = gridspacing * wlen_1;
00188 Vector s, d;
00189 Vector rc = 0;
00190 double xmin, xmax, ymin, ymax, zmin, zmax;
00191 int is_periodic_x = (isperiodic & NL_MSM_PERIODIC_VEC1);
00192 int is_periodic_y = (isperiodic & NL_MSM_PERIODIC_VEC2);
00193 int is_periodic_z = (isperiodic & NL_MSM_PERIODIC_VEC3);
00194 int i;
00195
00196
00197
00198
00199 d = coord[0].position - c;
00200 s.x = ru * d;
00201 s.y = rv * d;
00202 s.z = rw * d;
00203 xmin = xmax = s.x;
00204 ymin = ymax = s.y;
00205 zmin = zmax = s.z;
00206 for (i = 1; i < numatoms; i++) {
00207 d = coord[i].position - c;
00208 s.x = ru * d;
00209 s.y = rv * d;
00210 s.z = rw * d;
00211 if (s.x < xmin) xmin = s.x;
00212 else if (s.x > xmax) xmax = s.x;
00213 if (s.y < ymin) ymin = s.y;
00214 else if (s.y > ymax) ymax = s.y;
00215 if (s.z < zmin) zmin = s.z;
00216 else if (s.z > zmax) zmax = s.z;
00217 }
00218 #if 0
00219 printf("*** xmin=%.4f xmax=%.4f\n", xmin, xmax);
00220 printf("*** ymin=%.4f ymax=%.4f\n", ymin, ymax);
00221 printf("*** zmin=%.4f zmax=%.4f\n", zmin, zmax);
00222 #endif
00223
00224 if ( ! is_periodic_x) {
00225 xmax += rpadx;
00226 xmin -= rpadx;
00227 if (rhx > 0) {
00228 double mupper = ceil(xmax / (2*rhx));
00229 double mlower = floor(xmin / (2*rhx));
00230 xmax = 2*rhx*mupper;
00231 xmin = 2*rhx*mlower;
00232 }
00233 rc.x = 0.5*(xmin + xmax);
00234 xlen = xmax - xmin;
00235 if (xlen < NL_ZERO_TOLERANCE) {
00236 xlen = 1;
00237 }
00238 else {
00239 inv_xlen = 1/xlen;
00240 }
00241 }
00242
00243 if ( ! is_periodic_y) {
00244 ymax += rpady;
00245 ymin -= rpady;
00246 if (rhy > 0) {
00247 double mupper = ceil(ymax / (2*rhy));
00248 double mlower = floor(ymin / (2*rhy));
00249 ymax = 2*rhy*mupper;
00250 ymin = 2*rhy*mlower;
00251 }
00252 rc.y = 0.5*(ymin + ymax);
00253 ylen = ymax - ymin;
00254 if (ylen < NL_ZERO_TOLERANCE) {
00255 ylen = 1;
00256 }
00257 else {
00258 inv_ylen = 1/ylen;
00259 }
00260 }
00261
00262 if ( ! is_periodic_z) {
00263 zmax += rpadz;
00264 zmin -= rpadz;
00265 if (rhz > 0) {
00266 double mupper = ceil(zmax / (2*rhz));
00267 double mlower = floor(zmin / (2*rhz));
00268 zmax = 2*rhz*mupper;
00269 zmin = 2*rhz*mlower;
00270 }
00271 rc.z = 0.5*(zmin + zmax);
00272 zlen = zmax - zmin;
00273 if (zlen < NL_ZERO_TOLERANCE) {
00274 zlen = 1;
00275 }
00276 else {
00277 inv_zlen = 1/zlen;
00278 }
00279 }
00280
00281 #if 0
00282 printf("xmin=%g xmax=%g\n", xmin, xmax);
00283 printf("ymin=%g ymax=%g\n", ymin, ymax);
00284 printf("zmin=%g zmax=%g\n", zmin, zmax);
00285 printf("xlen=%g ylen=%g zlen=%g\n", xlen, ylen, zlen);
00286 printf("rc_x=%g rc_y=%g rc_z=%g\n", rc.x, rc.y, rc.z);
00287 #endif
00288
00289
00290 c.x = (u.x*rc.x + v.x*rc.y + w.x*rc.z) + c.x;
00291 c.y = (u.y*rc.x + v.y*rc.y + w.y*rc.z) + c.y;
00292 c.z = (u.z*rc.x + v.z*rc.y + w.z*rc.z) + c.z;
00293
00294 #if 0
00295 printf("c_x=%g c_y=%g c_z=%g\n", c.x, c.y, c.z);
00296 #endif
00297
00298 u *= xlen;
00299 v *= ylen;
00300 w *= zlen;
00301
00302 ru *= inv_xlen;
00303 rv *= inv_ylen;
00304 rw *= inv_zlen;
00305 }
00306
00307
00308 void ComputeMsmSerialMgr::recvCoord(MsmSerialCoordMsg *msg) {
00309 if ( ! numSources ) {
00310 numSources = (PatchMap::Object())->numNodesWithPatches();
00311 coordMsgs = new MsmSerialCoordMsg*[numSources];
00312 for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
00313 numArrived = 0;
00314 numAtoms = Node::Object()->molecule->numAtoms;
00315 coord = new ComputeMsmSerialAtom[numAtoms];
00316 force = new MsmSerialForce[numAtoms];
00317 }
00318
00319 int i;
00320 for ( i=0; i < msg->numAtoms; ++i ) {
00321 coord[msg->coord[i].id] = msg->coord[i];
00322 }
00323
00324 coordMsgs[numArrived] = msg;
00325 ++numArrived;
00326
00327 if ( numArrived < numSources ) return;
00328 numArrived = 0;
00329
00330
00331 Lattice lattice = msg->lattice;
00332 SimParameters *simParams = Node::Object()->simParameters;
00333
00334 double energy = 0;
00335 double virial[3][3];
00336
00337 int rc = 0;
00338
00339 if ( ! msmsolver ) {
00340
00341
00342
00343 msmsolver = NL_msm_create();
00344 if ( ! msmsolver ) NAMD_die("unable to create MSM solver");
00345 double dielectric = simParams->dielectric;
00346 double cutoff = simParams->cutoff;
00347 double gridspacing = simParams->MSMGridSpacing;
00348 double padding = simParams->MSMPadding;
00349 int approx = simParams->MSMApprox;
00350 int split = simParams->MSMSplit;
00351 int nlevels = simParams->MSMLevels;
00352 int msmflags = 0;
00353 msmflags |= (lattice.a_p() ? NL_MSM_PERIODIC_VEC1 : 0);
00354 msmflags |= (lattice.b_p() ? NL_MSM_PERIODIC_VEC2 : 0);
00355 msmflags |= (lattice.c_p() ? NL_MSM_PERIODIC_VEC3 : 0);
00356 msmflags |= NL_MSM_COMPUTE_LONG_RANGE;
00357
00358
00359 rc = NL_msm_configure(msmsolver, gridspacing, approx, split, nlevels);
00360 if (rc) NAMD_die("unable to configure MSM solver");
00361 Vector u=lattice.a(), v=lattice.b(), w=lattice.c(), c=lattice.origin();
00362 Vector ru=lattice.a_r(), rv=lattice.b_r(), rw=lattice.c_r();
00363 if ((msmflags & NL_MSM_PERIODIC_ALL) != NL_MSM_PERIODIC_ALL) {
00364
00365 int isperiodic = (msmflags & NL_MSM_PERIODIC_ALL);
00366
00367 rescale_nonperiodic_cell(u, v, w, c, ru, rv, rw,
00368 isperiodic, numAtoms, coord, padding, gridspacing);
00369 }
00370 double vec1[3], vec2[3], vec3[3], center[3];
00371 vec1[0] = u.x;
00372 vec1[1] = u.y;
00373 vec1[2] = u.z;
00374 vec2[0] = v.x;
00375 vec2[1] = v.y;
00376 vec2[2] = v.z;
00377 vec3[0] = w.x;
00378 vec3[1] = w.y;
00379 vec3[2] = w.z;
00380 center[0] = c.x;
00381 center[1] = c.y;
00382 center[2] = c.z;
00383 #if 0
00384 printf("dielectric = %g\n", dielectric);
00385 printf("vec1 = %g %g %g\n", vec1[0], vec1[1], vec1[2]);
00386 printf("vec2 = %g %g %g\n", vec2[0], vec2[1], vec2[2]);
00387 printf("vec3 = %g %g %g\n", vec3[0], vec3[1], vec3[2]);
00388 printf("center = %g %g %g\n", center[0], center[1], center[2]);
00389 printf("cutoff = %g\n", cutoff);
00390 printf("numatoms = %d\n", numAtoms);
00391 #endif
00392 rc = NL_msm_setup(msmsolver, cutoff, vec1, vec2, vec3, center, msmflags);
00393 if (rc) NAMD_die("unable to set up MSM solver");
00394 msmcoord = new double[4*numAtoms];
00395 msmforce = new double[3*numAtoms];
00396 if (msmcoord==0 || msmforce==0) NAMD_die("can't allocate MSM atom buffers");
00397
00398 double celec = sqrt(COULOMB / simParams->dielectric);
00399 for (i = 0; i < numAtoms; i++) {
00400 msmcoord[4*i+3] = celec * coord[i].charge;
00401 }
00402 }
00403
00404
00405 for (i = 0; i < numAtoms; i++) {
00406 msmcoord[4*i ] = coord[i].position.x;
00407 msmcoord[4*i+1] = coord[i].position.y;
00408 msmcoord[4*i+2] = coord[i].position.z;
00409 }
00410 for (i = 0; i < numAtoms; i++) {
00411 msmforce[3*i ] = 0;
00412 msmforce[3*i+1] = 0;
00413 msmforce[3*i+2] = 0;
00414 }
00415 rc = NL_msm_compute_force(msmsolver, msmforce, &energy, msmcoord, numAtoms);
00416 if (rc) NAMD_die("error evaluating MSM forces");
00417 for (i = 0; i < numAtoms; i++) {
00418 force[i].x = msmforce[3*i ];
00419 force[i].y = msmforce[3*i+1];
00420 force[i].z = msmforce[3*i+2];
00421 }
00422
00423 for (int k=0; k < 3; k++) {
00424 for (int l=0; l < 3; l++) {
00425 virial[k][l] = 0;
00426 }
00427 }
00428
00429
00430 for (int j=0; j < numSources; j++) {
00431 MsmSerialCoordMsg *cmsg = coordMsgs[j];
00432 coordMsgs[j] = 0;
00433 MsmSerialForceMsg *fmsg = new (cmsg->numAtoms, 0) MsmSerialForceMsg;
00434
00435 for (int i=0; i < cmsg->numAtoms; i++) {
00436 fmsg->force[i] = force[cmsg->coord[i].id];
00437 }
00438
00439 if ( ! j ) {
00440 fmsg->energy = energy;
00441 for (int k=0; k < 3; k++) {
00442 for (int l=0; l < 3; l++) {
00443 fmsg->virial[k][l] = virial[k][l];
00444 }
00445 }
00446 }
00447 else {
00448 fmsg->energy = 0;
00449 for (int k=0; k < 3; k++) {
00450 for (int l=0; l < 3; l++) {
00451 fmsg->virial[k][l] = 0;
00452 }
00453 }
00454 }
00455
00456 msmProxy[cmsg->sourceNode].recvForce(fmsg);
00457 delete cmsg;
00458 }
00459 }
00460
00461 void ComputeMsmSerialMgr::recvForce(MsmSerialForceMsg *msg) {
00462 msmCompute->saveResults(msg);
00463 delete oldmsg;
00464 oldmsg = msg;
00465 }
00466
00467 void ComputeMsmSerial::saveResults(MsmSerialForceMsg *msg)
00468 {
00469 ResizeArrayIter<PatchElem> ap(patchList);
00470
00471 MsmSerialForce *results_ptr = msg->force;
00472
00473
00474 for (ap = ap.begin(); ap != ap.end(); ap++) {
00475 Results *r = (*ap).forceBox->open();
00476 Force *f = r->f[Results::slow];
00477 int numAtoms = (*ap).p->getNumAtoms();
00478
00479 for(int i=0; i<numAtoms; ++i) {
00480 f[i].x += results_ptr->x;
00481 f[i].y += results_ptr->y;
00482 f[i].z += results_ptr->z;
00483 ++results_ptr;
00484 }
00485
00486 (*ap).forceBox->close(&r);
00487 }
00488
00489 reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += msg->energy;
00490 reduction->item(REDUCTION_VIRIAL_SLOW_XX) += msg->virial[0][0];
00491 reduction->item(REDUCTION_VIRIAL_SLOW_XY) += msg->virial[0][1];
00492 reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += msg->virial[0][2];
00493 reduction->item(REDUCTION_VIRIAL_SLOW_YX) += msg->virial[1][0];
00494 reduction->item(REDUCTION_VIRIAL_SLOW_YY) += msg->virial[1][1];
00495 reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += msg->virial[1][2];
00496 reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += msg->virial[2][0];
00497 reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += msg->virial[2][1];
00498 reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += msg->virial[2][2];
00499 reduction->submit();
00500 }
00501
00502 #include "ComputeMsmSerialMgr.def.h"
00503