ComputeGlobal.C

Go to the documentation of this file.
00001 
00007 /*
00008    Forwards atoms to master node for force evaluation.
00009 */
00010 
00011 #include "InfoStream.h"
00012 #include "Node.h"
00013 #include "PatchMap.h"
00014 #include "PatchMap.inl"
00015 #include "AtomMap.h"
00016 #include "ComputeGlobal.h"
00017 #include "ComputeGlobalMsgs.h"
00018 #include "PatchMgr.h"
00019 #include "Molecule.h"
00020 #include "ReductionMgr.h"
00021 #include "ComputeMgr.h"
00022 #include "ComputeMgr.decl.h"
00023 #include "SimParameters.h"
00024 #include <stdio.h>
00025 #include <algorithm>
00026 
00027 //#define DEBUGM
00028 #define MIN_DEBUG_LEVEL 1
00029 #include "Debug.h"
00030 
00031 // CLIENTS
00032 
00033 ComputeGlobal::ComputeGlobal(ComputeID c, ComputeMgr *m)
00034         : ComputeHomePatches(c)
00035 {
00036   DebugM(3,"Constructing client\n");
00037   aid.resize(0);
00038   gdef.resize(0);
00039   comm = m;
00040   firsttime = 1;
00041   isRequested = 0;
00042   isRequestedAllocSize = 0;
00043   endRequested = 0;
00044   numGroupsRequested = 0;
00045   SimParameters *sp = Node::Object()->simParameters;
00046   dofull = (sp->GBISserOn || sp->GBISOn || sp->fullDirectOn || sp->FMAOn || sp->PMEOn);
00047   forceSendEnabled = 0;
00048   if ( sp->tclForcesOn ) forceSendEnabled = 1;
00049   if ( sp->colvarsOn ) forceSendEnabled = 1;
00050   forceSendActive = 0;
00051   fid.resize(0);
00052   totalForce.resize(0);
00053   gfcount = 0;
00054   groupTotalForce.resize(0);
00055   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00056   int numPatches = PatchMap::Object()->numPatches();
00057   forcePtrs = new Force*[numPatches];
00058   atomPtrs = new FullAtom*[numPatches];
00059   for ( int i = 0; i < numPatches; ++i ) { forcePtrs[i] = 0; atomPtrs[i] = 0; }
00060 }
00061 
00062 ComputeGlobal::~ComputeGlobal()
00063 {
00064   delete[] isRequested;
00065   delete[] forcePtrs;
00066   delete[] atomPtrs;
00067   delete reduction;
00068 }
00069 
00070 void ComputeGlobal::configure(AtomIDList &newaid, AtomIDList &newgdef) {
00071   DebugM(4,"Receiving configuration (" << newaid.size() <<
00072         " atoms and " << newgdef.size() << " atoms/groups) on client\n");
00073 
00074   AtomIDList::iterator a, a_e;
00075   
00076  if ( forceSendEnabled ) {
00077   // clear previous data
00078   int max = -1;
00079   for (a=newaid.begin(),a_e=newaid.end(); a!=a_e; ++a) {
00080     if ( *a > max ) max = *a;
00081   }
00082   for (a=newgdef.begin(),a_e=newgdef.end(); a!=a_e; ++a) {
00083     if ( *a > max ) max = *a;
00084   }
00085   endRequested = max+1;
00086   if ( endRequested > isRequestedAllocSize ) {
00087     delete [] isRequested;
00088     isRequestedAllocSize = endRequested+10;
00089     isRequested = new char[isRequestedAllocSize];
00090     memset(isRequested, 0, isRequestedAllocSize);
00091   } else {
00092     for (a=aid.begin(),a_e=aid.end(); a!=a_e; ++a) {
00093       isRequested[*a] = 0;
00094     }
00095     for (a=gdef.begin(),a_e=gdef.end(); a!=a_e; ++a) {
00096       if ( *a != -1 ) isRequested[*a] = 0;
00097     }
00098   }
00099   // reserve space
00100   gpair.resize(0);
00101   gpair.resize(newgdef.size());
00102   gpair.resize(0);
00103  }
00104 
00105   // store data
00106   aid.swap(newaid);
00107   gdef.swap(newgdef);
00108   
00109  if ( forceSendEnabled ) {
00110   int newgcount = 0;
00111   for (a=aid.begin(),a_e=aid.end(); a!=a_e; ++a) {
00112     isRequested[*a] = 1;
00113   }
00114   for (a=gdef.begin(),a_e=gdef.end(); a!=a_e; ++a) {
00115     if ( *a == -1 ) ++newgcount;
00116     else {
00117       isRequested[*a] |= 2;
00118       gpair.add(intpair(*a,newgcount));
00119     }
00120   }
00121   std::sort(gpair.begin(),gpair.end());
00122   numGroupsRequested = newgcount;
00123  }
00124 }
00125 
00126 #if 0
00127 void ComputeGlobal::recvConfig(ComputeGlobalConfigMsg *msg) {
00128   DebugM(3,"Receiving configure on client\n");
00129   configure(msg->aid,msg->gdef);
00130   delete msg;
00131   sendData();
00132 }
00133 #endif
00134 
00135 void ComputeGlobal::recvResults(ComputeGlobalResultsMsg *msg) {
00136   DebugM(3,"Receiving results (" << msg->aid.size() << " forces, "
00137          << msg->newgdef.size() << " new group atoms) on client\n");
00138 
00139   forceSendActive = msg->totalforces;
00140   if ( forceSendActive && ! forceSendEnabled ) NAMD_bug("ComputeGlobal::recvResults forceSendActive without forceSendEnabled");
00141 
00142   // set the forces only if we aren't going to resend the data
00143   int setForces = !msg->resendCoordinates;
00144 
00145   if(setForces) { // we are requested to 
00146     // Store forces to patches
00147     AtomMap *atomMap = AtomMap::Object();
00148     const Lattice & lattice = patchList[0].p->lattice;
00149     ResizeArrayIter<PatchElem> ap(patchList);
00150     Force **f = forcePtrs;
00151     FullAtom **t = atomPtrs;
00152     Force extForce = 0.;
00153     Tensor extVirial;
00154 
00155     for (ap = ap.begin(); ap != ap.end(); ap++) {
00156       (*ap).r = (*ap).forceBox->open();
00157       f[(*ap).patchID] = (*ap).r->f[Results::normal];
00158       t[(*ap).patchID] = (*ap).p->getAtomList().begin();
00159     }
00160 
00161     AtomIDList::iterator a = msg->aid.begin();
00162     AtomIDList::iterator a_e = msg->aid.end();
00163     ForceList::iterator f2 = msg->f.begin();
00164     for ( ; a != a_e; ++a, ++f2 ) {
00165       DebugM(1,"processing atom "<<(*a)<<", F="<<(*f2)<<"...\n");
00166       /* XXX if (*a) is out of bounds here we get a segfault */
00167       LocalID localID = atomMap->localID(*a);
00168       if ( localID.pid == notUsed || ! f[localID.pid] ) continue;
00169       Force f_atom = (*f2);
00170       f[localID.pid][localID.index] += f_atom;
00171       FullAtom &atom = t[localID.pid][localID.index];
00172       Position x_orig = atom.position;
00173       Transform trans = atom.transform;
00174       Position x_atom = lattice.reverse_transform(x_orig,trans);
00175       extForce += f_atom;
00176       extVirial += outer(f_atom,x_atom);
00177     }
00178     DebugM(1,"done with the loop\n");
00179 
00180   // calculate forces for atoms in groups
00181     AtomIDList::iterator g_i, g_e;
00182     g_i = gdef.begin(); g_e = gdef.end();
00183     ForceList::iterator gf_i = msg->gforce.begin();
00184     //iout << iDEBUG << "recvResults\n" << endi;
00185     for ( ; g_i != g_e; ++g_i, ++gf_i ) {
00186       //iout << iDEBUG << *gf_i << '\n' << endi;
00187       Vector accel = (*gf_i);
00188       for ( ; *g_i != -1; ++g_i ) {
00189         //iout << iDEBUG << *g_i << '\n' << endi;
00190         LocalID localID = atomMap->localID(*g_i);
00191         if ( localID.pid == notUsed || ! f[localID.pid] ) continue;
00192         FullAtom &atom = t[localID.pid][localID.index];
00193         Force f_atom = accel * atom.mass;
00194         f[localID.pid][localID.index] += f_atom;
00195         Position x_orig = atom.position;
00196         Transform trans = atom.transform;
00197         Position x_atom = lattice.reverse_transform(x_orig,trans);
00198         extForce += f_atom;
00199         extVirial += outer(f_atom,x_atom);
00200       }
00201     }
00202     DebugM(1,"done with the groups\n");
00203 
00204     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
00205     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
00206     reduction->submit();
00207   }
00208   // done setting the forces, close boxes below
00209 
00210   // Get reconfiguration if present
00211   if ( msg->reconfig ) configure(msg->newaid, msg->newgdef);
00212 
00213   // send another round of data if requested
00214 
00215   if(msg->resendCoordinates) {
00216     DebugM(3,"Sending requested data right away\n");
00217     sendData();
00218   }
00219 
00220   groupTotalForce.resize(numGroupsRequested);
00221   for ( int i=0; i<numGroupsRequested; ++i ) groupTotalForce[i] = 0;
00222 
00223   if(setForces) {
00224     ResizeArrayIter<PatchElem> ap(patchList);
00225     Force **f = forcePtrs;
00226     FullAtom **t = atomPtrs;
00227     for (ap = ap.begin(); ap != ap.end(); ap++) {
00228       CompAtom *x;
00229       (*ap).positionBox->close(&x);
00230       (*ap).forceBox->close(&((*ap).r));
00231       f[(*ap).patchID] = 0;
00232       t[(*ap).patchID] = 0;
00233     }
00234   }
00235 
00236   delete msg;
00237   DebugM(3,"Done processing results\n");
00238 }
00239 
00240 void ComputeGlobal::doWork()
00241 {
00242   DebugM(2,"doWork\n");
00243 
00244   ResizeArrayIter<PatchElem> ap(patchList);
00245   FullAtom **t = atomPtrs;
00246 
00247   for (ap = ap.begin(); ap != ap.end(); ap++) {
00248     CompAtom *x = (*ap).positionBox->open();
00249     t[(*ap).patchID] = (*ap).p->getAtomList().begin();
00250   }
00251 
00252   if(!firsttime) sendData();
00253   else {
00254     if ( hasPatchZero ) {
00255       ComputeGlobalDataMsg *msg = new ComputeGlobalDataMsg;
00256       msg->lat.add(patchList[0].p->lattice);
00257       msg->step = -1;
00258       msg->count = 1;
00259       comm->sendComputeGlobalData(msg);
00260     }
00261     firsttime = 0;
00262     comm->enableComputeGlobalResults();
00263   }
00264   DebugM(2,"done with doWork\n");
00265 }
00266 
00267 void ComputeGlobal::sendData()
00268 {
00269   DebugM(2,"sendData\n");
00270   // Get positions from patches
00271   AtomMap *atomMap = AtomMap::Object();
00272   const Lattice & lattice = patchList[0].p->lattice;
00273   ResizeArrayIter<PatchElem> ap(patchList);
00274   FullAtom **t = atomPtrs;
00275 
00276   ComputeGlobalDataMsg *msg = new  ComputeGlobalDataMsg;
00277 
00278   msg->count = 0;
00279   msg->step = patchList[0].p->flags.step;
00280 
00281   AtomIDList::iterator a = aid.begin();
00282   AtomIDList::iterator a_e = aid.end();
00283   for ( ; a != a_e; ++a ) {
00284     LocalID localID = atomMap->localID(*a);
00285     if ( localID.pid == notUsed || ! t[localID.pid] ) continue;
00286     msg->aid.add(*a);
00287     msg->count++;
00288     FullAtom &atom = t[localID.pid][localID.index];
00289     Position x_orig = atom.position;
00290     Transform trans = atom.transform;
00291     msg->p.add(lattice.reverse_transform(x_orig,trans));
00292   }
00293 
00294   // calculate group centers of mass
00295   AtomIDList::iterator g_i, g_e;
00296   g_i = gdef.begin(); g_e = gdef.end();
00297   for ( ; g_i != g_e; ++g_i ) {
00298     Vector com(0,0,0);
00299     BigReal mass = 0.;
00300     for ( ; *g_i != -1; ++g_i ) {
00301       LocalID localID = atomMap->localID(*g_i);
00302       if ( localID.pid == notUsed || ! t[localID.pid] ) continue;
00303       msg->count++;
00304       FullAtom &atom = t[localID.pid][localID.index];
00305       Position x_orig = atom.position;
00306       Transform trans = atom.transform;
00307       com += lattice.reverse_transform(x_orig,trans) * atom.mass;
00308       mass += atom.mass;
00309     }
00310     DebugM(1,"Adding center of mass "<<com<<"\n");
00311     msg->gcom.add(com);
00312     msg->gmass.add(mass);
00313   }
00314 
00315   msg->fid.swap(fid);
00316   msg->tf.swap(totalForce);
00317   fid.resize(0);
00318   totalForce.resize(0);
00319 
00320   if ( gfcount ) msg->gtf.swap(groupTotalForce);
00321   msg->count += ( msg->fid.size() + gfcount );
00322   gfcount = 0;
00323 
00324   DebugM(3,"Sending data (" << msg->aid.size() << " positions) on client\n");
00325   if ( hasPatchZero ) { msg->count++;  msg->lat.add(lattice); }
00326   if ( msg->count ) comm->sendComputeGlobalData(msg);
00327   else delete msg;
00328   comm->enableComputeGlobalResults();
00329 }
00330 
00331 
00332 // This function is called by each HomePatch after force
00333 // evaluation. It stores the indices and forces of the requested
00334 // atoms here, to be sent to GlobalMasterServer during the next
00335 // time step. The total force is the sum of three components:
00336 // "normal", "nbond" and "slow", the latter two may be calculated
00337 // less frequently, so their most recent values are stored in
00338 // "f_saved" and used here. If we don't do full electrostatics,
00339 // there's no "slow" part.
00340 void ComputeGlobal::saveTotalForces(HomePatch *homePatch)
00341 {
00342   if ( ! forceSendEnabled ) NAMD_bug("ComputeGlobal::saveTotalForces called unexpectedly");
00343   if ( ! forceSendActive ) return;
00344 
00345   if ( Node::Object()->simParameters->accelMDOn && Node::Object()->simParameters->accelMDDebugOn && Node::Object()->simParameters->accelMDdihe ) {
00346     int num=homePatch->numAtoms;
00347     FullAtomList &atoms = homePatch->atom;
00348     ForceList &af=homePatch->f[Results::amdf];
00349 
00350     for (int i=0; i<num; ++i) {
00351       int index = atoms[i].id;
00352       if (index < endRequested && isRequested[index] & 1) {
00353         fid.add(index);
00354         totalForce.add(af[i]);
00355       }
00356     }
00357     return;
00358   }
00359 
00360   int fixedAtomsOn = Node::Object()->simParameters->fixedAtomsOn;
00361   int num=homePatch->numAtoms;
00362   FullAtomList &atoms = homePatch->atom;
00363   ForceList &f1=homePatch->f[Results::normal], &f2=homePatch->f_saved[Results::nbond],
00364             &f3=homePatch->f_saved[Results::slow];
00365   Force f_sum;
00366   
00367   for (int i=0; i<num; ++i) {
00368     int index = atoms[i].id;
00369     char reqflag;
00370     if (index < endRequested && (reqflag = isRequested[index])) {
00371      f_sum = f1[i]+f2[i];
00372      if (dofull)
00373        f_sum += f3[i];
00374      if ( fixedAtomsOn && atoms[i].atomFixed ) f_sum = 0.;
00375      if ( reqflag  & 1 ) {  // individual atom
00376       fid.add(index);
00377       totalForce.add(f_sum);
00378      }
00379      if ( reqflag  & 2 ) {  // part of group
00380        intpair *gpend = gpair.end();
00381        intpair *gpi = std::lower_bound(gpair.begin(),gpend,intpair(index,0));
00382        if ( gpi == gpend || gpi->first != index )
00383          NAMD_bug("ComputeGlobal::saveTotalForces gpair corrupted.");
00384        do {
00385          ++gfcount;
00386          groupTotalForce[gpi->second] += f_sum;
00387        } while ( ++gpi != gpend && gpi->first == index );
00388      }
00389     }
00390   }
00391 }

Generated on Thu Sep 21 01:17:10 2017 for NAMD by  doxygen 1.4.7