ComputeQM Class Reference

#include <ComputeQM.h>

Inheritance diagram for ComputeQM:

ComputeHomePatches Compute List of all members.

Public Member Functions

 ComputeQM (ComputeID c)
virtual ~ComputeQM ()
void initialize ()
void doWork ()
void saveResults (QMForceMsg *)
void processFullQM (QMCoordMsg *)

Detailed Description

Definition at line 114 of file ComputeQM.h.


Constructor & Destructor Documentation

ComputeQM::ComputeQM ( ComputeID  c  ) 

Definition at line 598 of file ComputeQM.C.

References ReductionMgr::Object(), REDUCTIONS_BASIC, and ReductionMgr::willSubmit().

00598                                 :
00599   ComputeHomePatches(c), oldForces(0)
00600 {
00601   CProxy_ComputeQMMgr::ckLocalBranch(
00602         CkpvAccess(BOCclass_group).computeQMMgr)->setCompute(this);
00603 
00604   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00605   
00606 }

ComputeQM::~ComputeQM (  )  [virtual]

Definition at line 608 of file ComputeQM.C.

00609 {
00610     if (oldForces != 0)
00611         delete [] oldForces;
00612 }


Member Function Documentation

void ComputeQM::doWork (  )  [virtual]

Reimplemented from Compute.

Definition at line 674 of file ComputeQM.C.

References ResizeArrayIter< T >::begin(), ComputeQMAtom::charge, QMCoordMsg::coord, DebugM, ResizeArrayIter< T >::end(), SortedArray< Elem >::find(), ComputeQMAtom::homeIndx, ComputeQMAtom::id, QMCoordMsg::numAtoms, ComputeHomePatches::patchList, ComputeQMAtom::position, ComputeQMAtom::qmGrpID, ResizeArray< Elem >::size(), QMCoordMsg::sourceNode, QMCoordMsg::timestep, CompAtom::vdwType, ComputeQMAtom::vdwType, and x.

00675 {
00676     CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
00677     
00678     ResizeArrayIter<PatchElem> ap(patchList);
00679     
00680         int timeStep;
00681         
00682     #ifdef DEBUG_QM
00683     DebugM(4,"----> Initiating QM work on rank " << CkMyPe() <<
00684     " with " << patchList.size() << " patches." << std::endl );
00685     #endif
00686     
00687     patchData.clear();
00688     
00689     // Determines hou many QM atoms are in this node.
00690     int numLocalQMAtoms = 0, localMM1atoms = 0;
00691     for (ap = ap.begin(); ap != ap.end(); ap++) {
00692         CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
00693         int localNumAtoms = (*ap).p->getNumAtoms();
00694         
00695         patchData.push_back(patchDataStrc((*ap).positionBox,
00696                     (*ap).positionBox->open(),
00697                     (*ap).p) );
00698         
00699         for(int i=0; i<localNumAtoms; ++i) {
00700             if ( qmAtomGroup[xExt[i].id] > 0 ) {
00701                 numLocalQMAtoms++;
00702             }
00703             else if (meNumMMIndx > 0) {
00704                 // If we have a mechanical embedding simulation with no 
00705                 // point charges, but with QM-MM bonds, we look for the MM1 atoms 
00706                 // here and pass them directly. No need for an extra round of 
00707                 // comms just to get atoms we already know we need.
00708                 
00709                 auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
00710                 
00711                 if (retIt != NULL) {
00712                     localMM1atoms++;
00713                 }
00714             }
00715         }
00716         
00717         timeStep = (*ap).p->flags.step ;
00718     }
00719     
00720     DebugM(4, "Node " << CkMyPe() << " has " << numLocalQMAtoms
00721     << " QM atoms." << std::endl) ;
00722     #ifdef DEBUG_QM
00723     if (localMM1atoms > 0)
00724         DebugM(4, "Node " << CkMyPe() << " has " << localMM1atoms
00725             << " MM1 atoms." << std::endl) ;
00726     #endif
00727     // Send QM atoms
00728     
00729     // This pointer will be freed in rank zero.
00730     QMCoordMsg *msg = new (numLocalQMAtoms + localMM1atoms,0, 0) QMCoordMsg;
00731     msg->sourceNode = CkMyPe();
00732     msg->numAtoms = numLocalQMAtoms + localMM1atoms;
00733     msg->timestep = timeStep;
00734     ComputeQMAtom *data_ptr = msg->coord;
00735     
00736     // Build a message with coordinates, charge, atom ID and QM group
00737     // for all QM atoms in the node, then send it to node zero.
00738     int homeCounter = 0;
00739     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
00740         
00741         CompAtom *x = (*pdIt).compAtomP;
00742         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
00743         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
00744         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
00745         
00746         for(int i=0; i<localNumAtoms; ++i) {
00747             
00748             if ( qmAtomGroup[xExt[i].id] > 0 ) {
00749                 
00750                 Real charge = 0;
00751                 
00752                 for (int qi=0; qi<numQMAtms; qi++) {
00753                     if (qmAtmIndx[qi] == xExt[i].id) {
00754                         charge = qmAtmChrg[qi];
00755                         break;
00756                     }
00757                 }
00758                 
00759                 data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position, 
00760                                                           fullAtms[i].transform) ;
00761 //                 data_ptr->charge = x[i].charge;
00762                 data_ptr->charge = charge;
00763                 data_ptr->id = xExt[i].id;
00764                 data_ptr->qmGrpID = qmAtomGroup[xExt[i].id] ;
00765                 data_ptr->homeIndx = homeCounter;
00766                 data_ptr->vdwType = fullAtms[i].vdwType;
00767                 
00768                 ++data_ptr;
00769             }
00770             else if (meNumMMIndx > 0) {
00771                 // If we have a mechanical embedding simulation with no 
00772                 // point charges, but with QM-MM bonds, we look for the MM1 atoms 
00773                 // connected here and pass them directly.
00774                 
00775                 auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
00776                 
00777                 if (retIt != NULL) {
00778                     
00779                     DebugM(3,"Found atom " << retIt->mmIndx << "," << retIt->qmGrp << "\n" );
00780                     
00781                     data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position, 
00782                                                           fullAtms[i].transform) ;
00783                     // We force the charge to be zero since we would only get
00784                     // here if mechanical embeding was selected.
00785                     data_ptr->charge = 0;
00786                     data_ptr->id = xExt[i].id;
00787                     data_ptr->qmGrpID = retIt->qmGrp ;
00788                     data_ptr->homeIndx = homeCounter;
00789                     
00790                     // We re-use this varaible to indicate this is an MM1 atom.
00791                     data_ptr->vdwType = -1;
00792                     
00793                     ++data_ptr;
00794                     
00795                 }
00796                 
00797             }
00798                 
00799             homeCounter++;
00800             
00801         }
00802         
00803     }
00804     
00805     if (noPC) {
00806         for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
00807             CompAtom *x = (*pdIt).compAtomP;
00808             (*pdIt).posBoxP->close(&x);
00809         }
00810     }
00811     
00812     QMProxy[0].recvPartQM(msg);
00813     
00814 }

void ComputeQM::initialize (  )  [virtual]

Reimplemented from ComputeHomePatches.

Definition at line 614 of file ComputeQM.C.

References SortedArray< Elem >::add(), SimParameters::cutoff, Molecule::get_noPC(), Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), Molecule::get_qmCustomPCIdxs(), Molecule::get_qmCustPCSizes(), Molecule::get_qmGrpID(), Molecule::get_qmMeMMindx(), Molecule::get_qmMeNumBonds(), Molecule::get_qmMeQMGrp(), Molecule::get_qmNumGrps(), Molecule::get_qmTotCustPCs(), ComputeHomePatches::initialize(), Node::molecule, Node::Object(), SimParameters::qmCustomPCSel, ResizeArray< Elem >::resize(), and Node::simParameters.

00615 {
00616     ComputeHomePatches::initialize();
00617     
00618     // Get some pointers that will be used in the future,
00619     simParams = Node::Object()->simParameters ;
00620     molPtr = Node::Object()->molecule;
00621     // Future comes fast...
00622     qmAtomGroup = molPtr->get_qmAtomGroup() ;
00623     
00624     noPC = molPtr->get_noPC();
00625     if (noPC) {
00626         meNumMMIndx = molPtr->get_qmMeNumBonds();
00627         meMMindx = molPtr->get_qmMeMMindx();
00628         meQMGrp =molPtr->get_qmMeQMGrp();
00629         
00630         for (int i=0; i<meNumMMIndx; i++)
00631             meQMBonds.add(meMMQMGrp(meMMindx[i], meQMGrp[i]));
00632     }
00633     
00634     numQMAtms = molPtr->get_numQMAtoms();
00635     qmAtmChrg = molPtr->get_qmAtmChrg() ;
00636     qmAtmIndx = molPtr->get_qmAtmIndx() ;
00637     
00638     numQMGrps = molPtr->get_qmNumGrps();
00639     qmGrpIDArray = molPtr->get_qmGrpID() ;
00640     
00641     cutoff = simParams->cutoff;
00642     
00643     customPC = simParams->qmCustomPCSel;
00644     if (customPC) {
00645         
00646         customPCLists.resize(numQMGrps);
00647         
00648         int totCustPCs = molPtr->get_qmTotCustPCs();
00649         const int *custPCSizes = molPtr->get_qmCustPCSizes();
00650         const int *customPCIdxs = molPtr->get_qmCustomPCIdxs();
00651         
00652         int minI = 0, maxI = 0, grpIter = 0;
00653         
00654         while (grpIter < numQMGrps) {
00655             
00656             maxI += custPCSizes[grpIter];
00657             
00658             for (size_t i=minI; i < maxI; i++) {
00659                 
00660                 customPCLists[grpIter].add( customPCIdxs[i] ) ;
00661             }
00662             
00663             // Minimum index gets the size of the current group
00664             minI += custPCSizes[grpIter];
00665             // Group iterator gets incremented
00666             grpIter++;
00667             
00668         }
00669     }
00670     
00671 }

void ComputeQM::processFullQM ( QMCoordMsg  ) 

Definition at line 1290 of file ComputeQM.C.

References ResizeArray< Elem >::add(), ResizeArray< Elem >::clear(), QMCoordMsg::coord, DebugM, SortedArray< Elem >::find(), if(), j, Vector::length(), SortedArray< Elem >::load(), FullAtom::mass, QMCoordMsg::numAtoms, QMCoordMsg::numPCIndxs, ComputeHomePatches::patchList, QMCoordMsg::pcIndxs, PCMODECUSTOMSEL, PCMODEUPDATEPOS, PCMODEUPDATESEL, QMCoordMsg::pcSelMode, ResizeArray< Elem >::size(), SortedArray< Elem >::sort(), QMCoordMsg::sourceNode, FullAtom::transform, CompAtom::vdwType, and x.

Referenced by ComputeQMMgr::recvFullQM().

01290                                                    {
01291     
01292     ResizeArrayIter<PatchElem> ap(patchList); 
01293     
01294     // Dynamically accumulates point charges as they are found.
01295     // The same MM atom may be added to this vector more than once if it is 
01296     // within the cutoff region of two or more QM regions.
01297     ResizeArray<ComputeQMPntChrg> compPCVec ;
01298     
01299     // This will keep the set of QM groups with which each 
01300     // point charge should be used. It is re-used for each atom in this
01301     // patch so we can controll the creation of the compPCVec vector.
01302     // The first item in the pair is the QM Group ID, the second is the shortest
01303     // distance between the point charge and an atom of this QM group.
01304     GrpDistMap chrgGrpMap ;
01305     
01306     DebugM(4,"Rank " << CkMyPe() << " receiving from rank " << qmFullMsg->sourceNode
01307     << " a total of " << qmFullMsg->numAtoms << " QM atoms and " 
01308     << qmFullMsg->numPCIndxs << " PC IDs." << std::endl);
01309     
01310     // If this is the firts step after a step of PC re-selection,
01311     // we store all PC IDs that all PEs found, in case some atoms end up
01312     // here untill the next re-selection step.
01313     if (qmFullMsg->numPCIndxs) {
01314         
01315         pcIDSortList.clear();
01316         
01317         int *pcIndx = qmFullMsg->pcIndxs;
01318         for (int i=0; i< qmFullMsg->numPCIndxs;i++) {
01319             pcIDSortList.load(pcIndx[i]);
01320         }
01321         
01322         pcIDSortList.sort();
01323     }
01324     
01325     int totalNodeAtoms = 0;
01326     int atomIter = 0;
01327     int uniqueCounter = 0;
01328     
01329     switch ( qmFullMsg->pcSelMode ) {
01330     
01331     case PCMODEUPDATESEL:
01332     {
01333         
01334     DebugM(4,"Updating PC selection.\n")
01335     
01336     // Loops over all atoms in this node and checks if any MM atom is within 
01337     // the cutof radius from a QM atom
01338     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
01339         
01340         CompAtom *x = (*pdIt).compAtomP;
01341         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
01342         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
01343         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
01344         
01345         totalNodeAtoms += localNumAtoms;
01346         
01347         // Iterates over the local atoms in a PE, all patches.
01348         for(int i=0; i<localNumAtoms; ++i) {
01349             
01350             const ComputeQMAtom *qmDataPtn = qmFullMsg->coord;
01351             
01352             // A new subset for each atom in this node.
01353             chrgGrpMap.clear();
01354             
01355             Position posMM = patchList[0].p->lattice.reverse_transform(x[i].position, 
01356                                                           fullAtms[i].transform) ;
01357             
01358             Real pcGrpID = qmAtomGroup[xExt[i].id];
01359             
01360             Real charge = x[i].charge;
01361             
01362             // If the atom is a QM atom, there will be no charge info in the 
01363             // atom box. Therefore, we take the charge in the previous step
01364             // in case this atom is a point charge for a neighboring QM region.
01365             if (pcGrpID > 0) {
01366                 for (int qi=0; qi<numQMAtms; qi++) {
01367                     if (qmAtmIndx[qi] == xExt[i].id) {
01368                         charge = qmAtmChrg[qi];
01369                         break;
01370                     }
01371                     
01372                 }
01373             }
01374             
01375             for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
01376                 
01377                 Real qmGrpID = qmDataPtn->qmGrpID;
01378                 
01379                 // We check If a QM atom and the node atom "i"
01380                 // belong to the same QM group. If not, atom "i" is either an MM
01381                 // atom or a QM atom from another group, which will be seen
01382                 // as an MM point charge.
01383                 // If they belong to the same group, just skip the distance 
01384                 // calculation and move on to the next QM atom.
01385                 // This loop needs to go over all QM atoms since a point charge
01386                 // may me be close to two different QM groups, in which case it 
01387                 // will be sent to both.
01388                 if (qmGrpID == pcGrpID) {
01389                     continue;
01390                 }
01391                 
01392                 Position posQM = qmDataPtn->position;
01393                 
01394                 Vector distV = posMM - posQM;
01395                 BigReal dist = distV.length();
01396                 
01397                 if ( dist < cutoff ){
01398                     
01399                     auto ret = chrgGrpMap.find(qmGrpID) ;
01400                     
01401                     // If 'ret' points to end, it means the item was not found
01402                     // and will be added, which means a new QM region has this 
01403                     // atom within its cutoff region,
01404                     // which means a new point charge will be 
01405                     // created in the QM system in question.
01406                     // 'ret' means a lot to me!
01407                     if ( ret ==  chrgGrpMap.end()) {
01408                         chrgGrpMap.insert(std::pair<Real,BigReal>(qmGrpID, dist));
01409                     }
01410                     else {
01411                         // If the current distance is smaller than the one in 
01412                         // the map, we update it so that we have the smallest
01413                         // distance between the point charge and any QM atom in 
01414                         // this group.
01415                         if (dist < ret->second) {
01416                             ret->second = dist;
01417                         }
01418                     }
01419                     
01420                 }
01421             }
01422             
01423             for(auto mapIt = chrgGrpMap.begin(); 
01424                 mapIt != chrgGrpMap.end(); mapIt++) {
01425                 
01426                 // We now add the final info about this point charge 
01427                 // to the vector, repeating it for each QM group that has it
01428                 // within its cuttoff radius.
01429                 compPCVec.add(
01430                     ComputeQMPntChrg(posMM, charge, xExt[i].id,
01431                                   mapIt->first, atomIter, mapIt->second, 
01432                                   fullAtms[i].mass, fullAtms[i].vdwType)
01433                                );
01434             }
01435             
01436             // Counts how many atoms are seens as point charges, by one or more
01437             // QM groups.
01438             if (chrgGrpMap.size() > 0)
01439                 uniqueCounter++;
01440             
01441             atomIter++;
01442         }
01443         
01444         (*pdIt).posBoxP->close(&x);
01445     }
01446     
01447     }break; // End case PCMODEUPDATESEL
01448     
01449     case PCMODEUPDATEPOS: 
01450     {
01451     
01452     DebugM(4,"Updating PC positions ONLY.\n")
01453     
01454     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
01455         
01456         CompAtom *x = (*pdIt).compAtomP;
01457         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
01458         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
01459         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
01460         
01461         totalNodeAtoms += localNumAtoms;
01462         
01463         // Iterates over the local atoms in a PE, all patches.
01464         for(int i=0; i<localNumAtoms; ++i) {
01465             
01466             if (pcIDSortList.find(xExt[i].id) != NULL ) {
01467                 Position posMM = patchList[0].p->lattice.reverse_transform(
01468                                                         x[i].position, 
01469                                                         fullAtms[i].transform) ;
01470                 
01471                 Real pcGrpID = qmAtomGroup[xExt[i].id];
01472                 Real charge = x[i].charge;
01473                 
01474                 // If the atom is a QM atom, there will be no charge info in the 
01475                 // atom box. Therefore, we take the charge in the previous step
01476                 // in case this atom is a point charge for a neighboring QM region.
01477                 if (pcGrpID > 0) {
01478                     for (int qi=0; qi<numQMAtms; qi++) {
01479                         if (qmAtmIndx[qi] == xExt[i].id) {
01480                             charge = qmAtmChrg[qi];
01481                             break;
01482                         }
01483                         
01484                     }
01485                 }
01486                 
01487                 compPCVec.add(
01488                     ComputeQMPntChrg(posMM, charge, xExt[i].id,
01489                                   0, atomIter, 0, fullAtms[i].mass, 
01490                                   fullAtms[i].vdwType));
01491             }
01492             
01493             atomIter++;
01494         }
01495         
01496         (*pdIt).posBoxP->close(&x);
01497     }
01498     }break ; // End case PCMODEUPDATEPOS
01499     
01500     case PCMODECUSTOMSEL:
01501     {
01502     
01503     DebugM(4,"Updating PC positions for custom PC selection.\n")
01504     
01505     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
01506         
01507         CompAtom *x = (*pdIt).compAtomP;
01508         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
01509         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
01510         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
01511         
01512         totalNodeAtoms += localNumAtoms;
01513         
01514         // Iterates over the local atoms in a PE, all patches.
01515         for(int i=0; i<localNumAtoms; ++i) {
01516             
01517             // Checks if the atom ID belongs to a point charge list,
01518             // which indicates it is to be sent to rank zero for QM calculations.
01519             // With one list per QM group, the same point charge can be added to 
01520             // different QM groups, as if it was being dynamically selected.
01521             for (int grpIndx=0; grpIndx<numQMGrps; grpIndx++) {
01522                 
01523                 if (customPCLists[grpIndx].find(xExt[i].id) != NULL){
01524                     
01525                     Position posMM = patchList[0].p->lattice.reverse_transform(
01526                                                         x[i].position, 
01527                                                         fullAtms[i].transform) ;
01528                     
01529                     Real pcGrpID = qmAtomGroup[xExt[i].id];
01530                     Real charge = x[i].charge;
01531                     
01532                     // If the atom is a QM atom, there will be no charge info in the 
01533                     // atom box. Therefore, we take the charge in the previous step
01534                     // in case this atom is a point charge for a neighboring QM region.
01535                     if (pcGrpID > 0) {
01536                         for (int qi=0; qi<numQMAtms; qi++) {
01537                             if (qmAtmIndx[qi] == xExt[i].id) {
01538                                 charge = qmAtmChrg[qi];
01539                                 break;
01540                             }
01541                             
01542                         }
01543                     }
01544                     
01545                     compPCVec.add(
01546                         ComputeQMPntChrg(posMM, charge, xExt[i].id,
01547                                       qmGrpIDArray[grpIndx], atomIter, 
01548                                       0, fullAtms[i].mass, fullAtms[i].vdwType));
01549                     
01550                 }
01551                 
01552             }
01553             
01554             atomIter++;
01555         }
01556         
01557         (*pdIt).posBoxP->close(&x);
01558     }
01559     
01560     } break;
01561     }
01562     
01563     DebugM(4,"Rank " << CkMyPe() << " found a total of " << compPCVec.size()
01564     << " point charges, out of " << totalNodeAtoms
01565     << " atoms in this node. " << uniqueCounter << " are unique." << std::endl);
01566     
01567     // Send only the MM atoms within radius
01568     
01569     // this pointer should be freed in rank zero, after receiving it.
01570     QMPntChrgMsg *pntChrgMsg = new (compPCVec.size(), 0) QMPntChrgMsg;
01571     pntChrgMsg->sourceNode = CkMyPe();
01572     pntChrgMsg->numAtoms = compPCVec.size();
01573     
01574     for (int i=0; i<compPCVec.size(); i++ ) {
01575         
01576         // Only sends the positions and charges of atoms within
01577         // cutoff of a (any) QM atom. Repeats for the number of QM groups
01578         // this charge is near to, and indicates the QM group it should 
01579         // be used with.
01580         pntChrgMsg->coord[i] = compPCVec[i] ;
01581         
01582     }
01583     
01584     DebugM(4,"Rank " << pntChrgMsg->sourceNode << " sending a total of " 
01585     << compPCVec.size() << " elements to rank zero." << std::endl);
01586     
01587     CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
01588     QMProxy[0].recvPntChrg(pntChrgMsg);
01589     
01590     delete qmFullMsg;
01591 }

void ComputeQM::saveResults ( QMForceMsg  ) 

Definition at line 2577 of file ComputeQM.C.

References ResizeArrayIter< T >::begin(), QMForce::charge, DebugM, ResizeArrayIter< T >::end(), QMForce::force, QMForceMsg::force, Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), QMForce::homeIndx, QMForce::id, if(), Node::molecule, QMForceMsg::numForces, Node::Object(), ComputeHomePatches::patchList, QMForceMsg::PMEOn, and QMForce::replace.

Referenced by ComputeQMMgr::recvForce().

02577                                             {
02578     
02579     ResizeArrayIter<PatchElem> ap(patchList);
02580     
02581     bool callReplaceForces = false;
02582     
02583     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
02584     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
02585     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
02586     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
02587     
02588     int totalNumbAtoms = 0;
02589     for (ap = ap.begin(); ap != ap.end(); ap++) {
02590         totalNumbAtoms += (*ap).p->getNumAtoms();
02591     }
02592     
02593     // This is kept to be deleted in the next step so that the message can be 
02594     // used in "replaceForces" routine later on, in case there are forces to 
02595     // be replaced. The "replaceForces" function uses but DOES NOT free the pointer,
02596     // so we free the data from the previous iteration and allocate a new one for
02597     // the current iteration (remember that the number of atoms can change in a 
02598     // home patch between iterations).
02599     if (oldForces != 0)
02600         delete [] oldForces;
02601     oldForces = new ExtForce[totalNumbAtoms] ;
02602     
02603     for (int i=0; i < totalNumbAtoms; ++i) {
02604         oldForces[i].force = Force(0) ;
02605     }
02606     
02607     DebugM(1,"Force array has been created and zeroed in rank " 
02608     << CkMyPe() << std::endl);
02609     
02610     DebugM(1,"Preparing " << fmsg->numForces << " forces in rank " 
02611     << CkMyPe() << std::endl);
02612     
02613     QMForce *results_ptr = fmsg->force;
02614     // Iterates over the received forces and place them on the full array.
02615     for (int i=0; i < fmsg->numForces; ++i, ++results_ptr) {
02616         // For now we may have more than one item in the message acting on the same 
02617         // atom, such as an MM atom which is a point charge for two or more QM regions.
02618         
02619         oldForces[results_ptr->homeIndx].force += results_ptr->force;
02620         oldForces[results_ptr->homeIndx].replace = results_ptr->replace;
02621         
02622         if (results_ptr->replace == 1)
02623             callReplaceForces = true;
02624         
02625         // If the atom is in a QM group, update its charge to the local (this homePatch)
02626         // copy of the qmAtmChrg array.
02627         if (qmAtomGroup[results_ptr->id] > 0 && fmsg->PMEOn) {
02628             
02629             // Loops over all QM atoms (in all QM groups) comparing their global indices
02630             for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
02631                 
02632                 if (qmAtmIndx[qmIter] == results_ptr->id) {
02633                     qmAtmChrg[qmIter] = results_ptr->charge;
02634                     break;
02635                 }
02636                 
02637             }
02638             
02639         }
02640         
02641     }
02642     
02643     DebugM(1,"Placing forces on force boxes in rank " 
02644     << CkMyPe() << std::endl);
02645     
02646     // Places the received forces on the force array for each patch.
02647     int homeIndxIter = 0;
02648     for (ap = ap.begin(); ap != ap.end(); ap++) {
02649         Results *r = (*ap).forceBox->open();
02650         Force *f = r->f[Results::normal];
02651         const CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
02652         int localNumAtoms = (*ap).p->getNumAtoms();
02653         
02654         for(int i=0; i<localNumAtoms; ++i) {
02655             
02656             f[i] += oldForces[homeIndxIter].force;
02657             
02658             ++homeIndxIter;
02659         }
02660         
02661         if ( callReplaceForces )
02662             (*ap).p->replaceForces(oldForces);
02663         
02664         (*ap).forceBox->close(&r);
02665         
02666     }
02667     
02668     DebugM(1,"Forces placed on force boxes in rank " 
02669     << CkMyPe() << std::endl);
02670     
02671     if (fmsg->PMEOn) {
02672         
02673         DebugM(1,"PME ON! Accessing PMEmgr in rank " << CkMyPe() << std::endl);
02674         
02675         ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch(
02676             CkpvAccess(BOCclass_group).computePmeMgr) ;
02677         
02678         DebugM(4, "Initiating ComputePme::doQMWork on rank " << CkMyPe() << " over " 
02679             << getComputes(mgrP).size() << " pmeComputes." << std::endl) ;
02680         
02681         for ( int i=0; i< getComputes(mgrP).size(); ++i ) {
02682     //         WorkDistrib::messageEnqueueWork(pmeComputes[i]);
02683             getComputes(mgrP)[i]->doQMWork();
02684         }
02685     }
02686     
02687     DebugM(1,"Submitting reduction in rank " << CkMyPe() << std::endl);
02688     
02689     reduction->item(REDUCTION_MISC_ENERGY) += fmsg->energy;
02690     reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += fmsg->virial[0][0];
02691     reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += fmsg->virial[0][1];
02692     reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += fmsg->virial[0][2];
02693     reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += fmsg->virial[1][0];
02694     reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += fmsg->virial[1][1];
02695     reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += fmsg->virial[1][2];
02696     reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += fmsg->virial[2][0];
02697     reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += fmsg->virial[2][1];
02698     reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += fmsg->virial[2][2];
02699     reduction->submit();
02700     
02701     delete fmsg ;
02702 }


The documentation for this class was generated from the following files:
Generated on Tue Nov 21 01:17:18 2017 for NAMD by  doxygen 1.4.7