00001
00007 #include "InfoStream.h"
00008 #include "Node.h"
00009 #include "PatchMap.h"
00010 #include "PatchMap.inl"
00011 #include "AtomMap.h"
00012 #include "ComputeExt.h"
00013 #include "ComputeExtMgr.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 <stdlib.h>
00026 #include <stdio.h>
00027 #include <errno.h>
00028
00029 #ifndef SQRT_PI
00030 #define SQRT_PI 1.7724538509055160273
00031 #endif
00032
00033 class ExtCoordMsg : public CMessage_ExtCoordMsg {
00034 public:
00035 int sourceNode;
00036 int numAtoms;
00037 Lattice lattice;
00038 CompAtom *coord;
00039 };
00040
00041 class ExtForceMsg : public CMessage_ExtForceMsg {
00042 public:
00043 BigReal energy;
00044 BigReal virial[3][3];
00045 ExtForce *force;
00046 };
00047
00048 class ComputeExtMgr : public BOCclass {
00049 public:
00050 ComputeExtMgr();
00051 ~ComputeExtMgr();
00052
00053 void setCompute(ComputeExt *c) { extCompute = c; }
00054
00055 void recvCoord(ExtCoordMsg *);
00056 void recvForce(ExtForceMsg *);
00057
00058 private:
00059 CProxy_ComputeExtMgr extProxy;
00060 ComputeExt *extCompute;
00061
00062 int numSources;
00063 int numArrived;
00064 ExtCoordMsg **coordMsgs;
00065 int numAtoms;
00066 CompAtom *coord;
00067 ExtForce *force;
00068 ExtForceMsg *oldmsg;
00069 };
00070
00071 ComputeExtMgr::ComputeExtMgr() :
00072 extProxy(thisgroup), extCompute(0), numSources(0), numArrived(0),
00073 coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0) {
00074 CkpvAccess(BOCclass_group).computeExtMgr = thisgroup;
00075 }
00076
00077 ComputeExtMgr::~ComputeExtMgr() {
00078 for ( int i=0; i<numSources; ++i ) { delete coordMsgs[i]; }
00079 delete [] coordMsgs;
00080 delete [] coord;
00081 delete [] force;
00082 delete oldmsg;
00083 }
00084
00085 ComputeExt::ComputeExt(ComputeID c) :
00086 ComputeHomePatches(c)
00087 {
00088 CProxy_ComputeExtMgr::ckLocalBranch(
00089 CkpvAccess(BOCclass_group).computeExtMgr)->setCompute(this);
00090
00091 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00092
00093 }
00094
00095 ComputeExt::~ComputeExt()
00096 {
00097 }
00098
00099 void ComputeExt::doWork()
00100 {
00101 ResizeArrayIter<PatchElem> ap(patchList);
00102
00103 #if 0
00104
00105 if ( ! patchList[0].p->flags.doFullElectrostatics )
00106 {
00107 for (ap = ap.begin(); ap != ap.end(); ap++) {
00108 CompAtom *x = (*ap).positionBox->open();
00109 Results *r = (*ap).forceBox->open();
00110 (*ap).positionBox->close(&x);
00111 (*ap).forceBox->close(&r);
00112 }
00113 reduction->submit();
00114 return;
00115 }
00116 #endif
00117
00118
00119 int numLocalAtoms = 0;
00120 for (ap = ap.begin(); ap != ap.end(); ap++) {
00121 numLocalAtoms += (*ap).p->getNumAtoms();
00122 }
00123
00124 ExtCoordMsg *msg = new (numLocalAtoms, 0) ExtCoordMsg;
00125 msg->sourceNode = CkMyPe();
00126 msg->numAtoms = numLocalAtoms;
00127 msg->lattice = patchList[0].p->flags.lattice;
00128 CompAtom *data_ptr = msg->coord;
00129
00130
00131 for (ap = ap.begin(); ap != ap.end(); ap++) {
00132 CompAtom *x = (*ap).positionBox->open();
00133 #if 0
00134 if ( patchList[0].p->flags.doMolly ) {
00135 (*ap).positionBox->close(&x);
00136 x = (*ap).avgPositionBox->open();
00137 }
00138 #endif
00139 int numAtoms = (*ap).p->getNumAtoms();
00140
00141 for(int i=0; i<numAtoms; ++i)
00142 {
00143 *data_ptr = x[i];
00144 ++data_ptr;
00145 }
00146
00147 #if 0
00148 if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00149 else { (*ap).positionBox->close(&x); }
00150 #endif
00151 (*ap).positionBox->close(&x);
00152 }
00153
00154 CProxy_ComputeExtMgr extProxy(CkpvAccess(BOCclass_group).computeExtMgr);
00155 #if CHARM_VERSION > 050402
00156 extProxy[0].recvCoord(msg);
00157 #else
00158 extProxy.recvCoord(msg,0);
00159 #endif
00160
00161 }
00162
00163 void ComputeExtMgr::recvCoord(ExtCoordMsg *msg) {
00164 if ( ! numSources ) {
00165 numSources = (PatchMap::Object())->numNodesWithPatches();
00166 coordMsgs = new ExtCoordMsg*[numSources];
00167 for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
00168 numArrived = 0;
00169 numAtoms = Node::Object()->molecule->numAtoms;
00170 coord = new CompAtom[numAtoms];
00171 force = new ExtForce[numAtoms];
00172 }
00173
00174 int i;
00175 for ( i=0; i < msg->numAtoms; ++i ) {
00176 coord[msg->coord[i].id] = msg->coord[i];
00177 }
00178
00179 coordMsgs[numArrived] = msg;
00180 ++numArrived;
00181
00182 if ( numArrived < numSources ) return;
00183 numArrived = 0;
00184
00185
00186 Lattice lattice = msg->lattice;
00187 SimParameters *simParams = Node::Object()->simParameters;
00188 FILE *file;
00189 int iret;
00190
00191
00192
00193 file = fopen(simParams->extCoordFilename,"w");
00194 if ( ! file ) { NAMD_die(strerror(errno)); }
00195 for ( i=0; i<numAtoms; ++i ) {
00196 int id = coord[i].id + 1;
00197 double charge = coord[i].charge;
00198 double x = coord[i].position.x;
00199 double y = coord[i].position.y;
00200 double z = coord[i].position.z;
00201 iret = fprintf(file,"%d %f %f %f %f\n",id,charge,x,y,z);
00202 if ( iret < 0 ) { NAMD_die(strerror(errno)); }
00203 }
00204
00205 Vector a = lattice.a(); if ( ! lattice.a_p() ) a = Vector(0,0,0);
00206 Vector b = lattice.b(); if ( ! lattice.b_p() ) b = Vector(0,0,0);
00207 Vector c = lattice.c(); if ( ! lattice.c_p() ) c = Vector(0,0,0);
00208 iret = fprintf(file,"%f %f %f\n%f %f %f\n%f %f %f\n",
00209 a.x, a.y, a.z, b.x, b.y, b.z, c.x, c.y, c.z);
00210 if ( iret < 0 ) { NAMD_die(strerror(errno)); }
00211 fclose(file);
00212
00213
00214
00215 iret = system(simParams->extForcesCommand);
00216 if ( iret == -1 ) { NAMD_die(strerror(errno)); }
00217 if ( iret ) { NAMD_die("Error running command for external forces."); }
00218
00219
00220 iret = remove(simParams->extCoordFilename);
00221 if ( iret ) { NAMD_die(strerror(errno)); }
00222
00223
00224
00225 file = fopen(simParams->extForceFilename,"r");
00226 if ( ! file ) { NAMD_die(strerror(errno)); }
00227 for ( i=0; i<numAtoms; ++i ) {
00228 int id, replace;
00229 double x, y, z;
00230 iret = fscanf(file,"%d %d %lf %lf %lf\n", &id, &replace, &x, &y, &z);
00231 if ( iret != 5 ) { NAMD_die("Error reading external forces file."); }
00232 if ( id != i + 1 ) { NAMD_die("Atom ID error in external forces file."); }
00233 force[i].force.x = x; force[i].force.y = y; force[i].force.z = z;
00234 force[i].replace = replace;
00235 }
00236
00237
00238
00239 double energy;
00240 double virial[3][3];
00241 iret = fscanf(file,"%lf\n", &energy);
00242 if ( iret != 1 ) {
00243 energy = 0;
00244 for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] = 0;
00245 } else {
00246 iret = fscanf(file,"%lf %lf %lf\n%lf %lf %lf\n%lf %lf %lf\n",
00247 &virial[0][0], &virial[0][1], &virial[0][2],
00248 &virial[1][0], &virial[1][1], &virial[1][2],
00249 &virial[2][0], &virial[2][1], &virial[2][2]);
00250 if ( iret != 9 ) {
00251 for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] = 0;
00252 } else {
00253
00254 for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] *= -1.0;
00255 }
00256 }
00257 fclose(file);
00258
00259
00260 iret = remove(simParams->extForceFilename);
00261 if ( iret ) { NAMD_die(strerror(errno)); }
00262
00263
00264
00265 for ( int j=0; j < numSources; ++j ) {
00266 ExtCoordMsg *cmsg = coordMsgs[j];
00267 coordMsgs[j] = 0;
00268 ExtForceMsg *fmsg = new (cmsg->numAtoms, 0) ExtForceMsg;
00269 for ( int i=0; i < cmsg->numAtoms; ++i ) {
00270 fmsg->force[i] = force[cmsg->coord[i].id];
00271 }
00272 if ( ! j ) {
00273 fmsg->energy = energy;
00274 for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l )
00275 fmsg->virial[k][l] = virial[k][l];
00276 } else {
00277 fmsg->energy = 0;
00278 for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l )
00279 fmsg->virial[k][l] = 0;
00280 }
00281 #if CHARM_VERSION > 050402
00282 extProxy[cmsg->sourceNode].recvForce(fmsg);
00283 #else
00284 extProxy.recvForce(fmsg,cmsg->sourceNode);
00285 #endif
00286 delete cmsg;
00287 }
00288
00289 }
00290
00291 void ComputeExtMgr::recvForce(ExtForceMsg *msg) {
00292 extCompute->saveResults(msg);
00293 delete oldmsg;
00294 oldmsg = msg;
00295 }
00296
00297 void ComputeExt::saveResults(ExtForceMsg *msg)
00298 {
00299 ResizeArrayIter<PatchElem> ap(patchList);
00300
00301 ExtForce *results_ptr = msg->force;
00302
00303
00304 for (ap = ap.begin(); ap != ap.end(); ap++) {
00305 Results *r = (*ap).forceBox->open();
00306 Force *f = r->f[Results::normal];
00307 int numAtoms = (*ap).p->getNumAtoms();
00308
00309 int replace = 0;
00310 ExtForce *replacementForces = results_ptr;
00311 for(int i=0; i<numAtoms; ++i) {
00312 if ( results_ptr->replace ) replace = 1;
00313 else f[i] += results_ptr->force;
00314 ++results_ptr;
00315 }
00316 if ( replace ) (*ap).p->replaceForces(replacementForces);
00317
00318 (*ap).forceBox->close(&r);
00319 }
00320
00321 reduction->item(REDUCTION_MISC_ENERGY) += msg->energy;
00322 reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += msg->virial[0][0];
00323 reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += msg->virial[0][1];
00324 reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += msg->virial[0][2];
00325 reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += msg->virial[1][0];
00326 reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += msg->virial[1][1];
00327 reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += msg->virial[1][2];
00328 reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += msg->virial[2][0];
00329 reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += msg->virial[2][1];
00330 reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += msg->virial[2][2];
00331 reduction->submit();
00332 }
00333
00334 #include "ComputeExtMgr.def.h"
00335