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 599 of file ComputeQM.C.

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

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

ComputeQM::~ComputeQM (  )  [virtual]

Definition at line 609 of file ComputeQM.C.

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


Member Function Documentation

void ComputeQM::doWork (  )  [virtual]

Reimplemented from Compute.

Definition at line 675 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.

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

void ComputeQM::initialize (  )  [virtual]

Reimplemented from ComputeHomePatches.

Definition at line 615 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.

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

void ComputeQM::processFullQM ( QMCoordMsg  ) 

Definition at line 1292 of file ComputeQM.C.

References ResizeArray< Elem >::add(), ResizeArray< Elem >::clear(), QMCoordMsg::coord, DebugM, Lattice::delta(), 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, ComputeQMAtom::position, ComputeQMAtom::qmGrpID, ResizeArray< Elem >::size(), SortedArray< Elem >::sort(), QMCoordMsg::sourceNode, CompAtom::vdwType, and x.

Referenced by ComputeQMMgr::recvFullQM().

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

void ComputeQM::saveResults ( QMForceMsg  ) 

Definition at line 2656 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().

02656                                             {
02657     
02658     ResizeArrayIter<PatchElem> ap(patchList);
02659     
02660     bool callReplaceForces = false;
02661     
02662     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
02663     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
02664     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
02665     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
02666     
02667     int totalNumbAtoms = 0;
02668     for (ap = ap.begin(); ap != ap.end(); ap++) {
02669         totalNumbAtoms += (*ap).p->getNumAtoms();
02670     }
02671     
02672     // This is kept to be deleted in the next step so that the message can be 
02673     // used in "replaceForces" routine later on, in case there are forces to 
02674     // be replaced. The "replaceForces" function uses but DOES NOT free the pointer,
02675     // so we free the data from the previous iteration and allocate a new one for
02676     // the current iteration (remember that the number of atoms can change in a 
02677     // home patch between iterations).
02678     if (oldForces != 0)
02679         delete [] oldForces;
02680     oldForces = new ExtForce[totalNumbAtoms] ;
02681     
02682     for (int i=0; i < totalNumbAtoms; ++i) {
02683         oldForces[i].force = Force(0) ;
02684     }
02685     
02686     DebugM(1,"Force array has been created and zeroed in rank " 
02687     << CkMyPe() << std::endl);
02688     
02689     DebugM(1,"Preparing " << fmsg->numForces << " forces in rank " 
02690     << CkMyPe() << std::endl);
02691     
02692     QMForce *results_ptr = fmsg->force;
02693     // Iterates over the received forces and place them on the full array.
02694     for (int i=0; i < fmsg->numForces; ++i, ++results_ptr) {
02695         // For now we may have more than one item in the message acting on the same 
02696         // atom, such as an MM atom which is a point charge for two or more QM regions.
02697         
02698         oldForces[results_ptr->homeIndx].force += results_ptr->force;
02699         oldForces[results_ptr->homeIndx].replace = results_ptr->replace;
02700         
02701         if (results_ptr->replace == 1)
02702             callReplaceForces = true;
02703         
02704         // If the atom is in a QM group, update its charge to the local (this homePatch)
02705         // copy of the qmAtmChrg array.
02706         if (qmAtomGroup[results_ptr->id] > 0 && fmsg->PMEOn) {
02707             
02708             // Loops over all QM atoms (in all QM groups) comparing their global indices
02709             for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
02710                 
02711                 if (qmAtmIndx[qmIter] == results_ptr->id) {
02712                     qmAtmChrg[qmIter] = results_ptr->charge;
02713                     break;
02714                 }
02715                 
02716             }
02717             
02718         }
02719         
02720     }
02721     
02722     DebugM(1,"Placing forces on force boxes in rank " 
02723     << CkMyPe() << std::endl);
02724     
02725     // Places the received forces on the force array for each patch.
02726     int homeIndxIter = 0;
02727     for (ap = ap.begin(); ap != ap.end(); ap++) {
02728         Results *r = (*ap).forceBox->open();
02729         Force *f = r->f[Results::normal];
02730         const CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
02731         int localNumAtoms = (*ap).p->getNumAtoms();
02732         
02733         for(int i=0; i<localNumAtoms; ++i) {
02734             
02735             f[i] += oldForces[homeIndxIter].force;
02736             
02737             ++homeIndxIter;
02738         }
02739         
02740         if ( callReplaceForces )
02741             (*ap).p->replaceForces(oldForces);
02742         
02743         (*ap).forceBox->close(&r);
02744         
02745     }
02746     
02747     DebugM(1,"Forces placed on force boxes in rank " 
02748     << CkMyPe() << std::endl);
02749     
02750     if (fmsg->PMEOn) {
02751         
02752         DebugM(1,"PME ON! Accessing PMEmgr in rank " << CkMyPe() << std::endl);
02753         
02754         ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch(
02755             CkpvAccess(BOCclass_group).computePmeMgr) ;
02756         
02757         DebugM(4, "Initiating ComputePme::doQMWork on rank " << CkMyPe() << " over " 
02758             << getComputes(mgrP).size() << " pmeComputes." << std::endl) ;
02759         
02760         for ( int i=0; i< getComputes(mgrP).size(); ++i ) {
02761     //         WorkDistrib::messageEnqueueWork(pmeComputes[i]);
02762             getComputes(mgrP)[i]->doQMWork();
02763         }
02764     }
02765     
02766     DebugM(1,"Submitting reduction in rank " << CkMyPe() << std::endl);
02767     
02768     reduction->item(REDUCTION_MISC_ENERGY) += fmsg->energy;
02769     reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += fmsg->virial[0][0];
02770     reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += fmsg->virial[0][1];
02771     reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += fmsg->virial[0][2];
02772     reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += fmsg->virial[1][0];
02773     reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += fmsg->virial[1][1];
02774     reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += fmsg->virial[1][2];
02775     reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += fmsg->virial[2][0];
02776     reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += fmsg->virial[2][1];
02777     reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += fmsg->virial[2][2];
02778     reduction->submit();
02779     
02780     delete fmsg ;
02781 }


The documentation for this class was generated from the following files:
Generated on Sat Feb 24 01:17:19 2018 for NAMD by  doxygen 1.4.7