ComputeExt.C

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

Generated on Sun Nov 19 01:17:11 2017 for NAMD by  doxygen 1.4.7