ComputeQM Class Reference

#include <ComputeQM.h>

Inheritance diagram for ComputeQM:

Public Member Functions
	ComputeQM (ComputeID c)
virtual	~ComputeQM ()
void	initialize ()
void	doWork ()
void	saveResults (QMForceMsg *)
void	processFullQM (QMCoordMsg *)
Public Member Functions inherited from ComputeHomePatches
	ComputeHomePatches (ComputeID c)
virtual	~ComputeHomePatches ()
virtual void	atomUpdate ()
Flags *	getFlags (void)
Public Member Functions inherited from Compute
	Compute (ComputeID)
int	type ()
virtual	~Compute ()
void	setNumPatches (int n)
int	getNumPatches ()
virtual void	patchReady (PatchID, int doneMigration, int seq)
virtual int	noWork ()
virtual void	finishPatch (int)
int	sequence (void)
int	priority (void)
int	getGBISPhase (void)
virtual void	gbisP2PatchReady (PatchID, int seq)
virtual void	gbisP3PatchReady (PatchID, int seq)

Additional Inherited Members
Public Attributes inherited from Compute
const ComputeID	cid
LDObjHandle	ldObjHandle
LocalWorkMsg *const	localWorkMsg
Protected Member Functions inherited from Compute
void	enqueueWork ()
Protected Attributes inherited from ComputeHomePatches
int	useAvgPositions
int	hasPatchZero
ComputeHomePatchList	patchList
PatchMap *	patchMap
Protected Attributes inherited from Compute
int	computeType
int	basePriority
int	gbisPhase
int	gbisPhasePriority [3]

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().

                                :
  ComputeHomePatches(c), oldForces(0)
{
  CProxy_ComputeQMMgr::ckLocalBranch(
        CkpvAccess(BOCclass_group).computeQMMgr)->setCompute(this);

  reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
  
}

ComputeQM::~ComputeQM ( )

virtual

Definition at line 609 of file ComputeQM.C.

{
    if (oldForces != 0)
        delete [] oldForces;
}

Member Function Documentation

void ComputeQM::doWork ( void  )

virtual

Reimplemented from Compute.

Definition at line 675 of file ComputeQM.C.

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

{
    CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
    
    ResizeArrayIter<PatchElem> ap(patchList);
    
        int timeStep;
        
    #ifdef DEBUG_QM
    DebugM(4,"----> Initiating QM work on rank " << CkMyPe() <<
    " with " << patchList.size() << " patches." << std::endl );
    #endif
    
    patchData.clear();
    
    // Determines hou many QM atoms are in this node.
    int numLocalQMAtoms = 0, localMM1atoms = 0;
    for (ap = ap.begin(); ap != ap.end(); ap++) {
        CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
        int localNumAtoms = (*ap).p->getNumAtoms();
        
        patchData.push_back(patchDataStrc((*ap).positionBox,
                    (*ap).positionBox->open(),
                    (*ap).p) );
        
        for(int i=0; i<localNumAtoms; ++i) {
            if ( qmAtomGroup[xExt[i].id] > 0 ) {
                numLocalQMAtoms++;
            }
            else if (meNumMMIndx > 0) {
                // If we have a mechanical embedding simulation with no 
                // point charges, but with QM-MM bonds, we look for the MM1 atoms 
                // here and pass them directly. No need for an extra round of 
                // comms just to get atoms we already know we need.
                
                auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
                
                if (retIt != NULL) {
                    localMM1atoms++;
                }
            }
        }
        
        timeStep = (*ap).p->flags.step ;
    }
    
    DebugM(4, "Node " << CkMyPe() << " has " << numLocalQMAtoms
    << " QM atoms." << std::endl) ;
    #ifdef DEBUG_QM
    if (localMM1atoms > 0)
        DebugM(4, "Node " << CkMyPe() << " has " << localMM1atoms
            << " MM1 atoms." << std::endl) ;
    #endif
    // Send QM atoms
    
    // This pointer will be freed in rank zero.
    QMCoordMsg *msg = new (numLocalQMAtoms + localMM1atoms,0, 0) QMCoordMsg;
    msg->sourceNode = CkMyPe();
    msg->numAtoms = numLocalQMAtoms + localMM1atoms;
    msg->timestep = timeStep;
    ComputeQMAtom *data_ptr = msg->coord;
    
    // Build a message with coordinates, charge, atom ID and QM group
    // for all QM atoms in the node, then send it to node zero.
    int homeCounter = 0;
    for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
        
        CompAtom *x = (*pdIt).compAtomP;
        const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
        int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
        const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
        
        for(int i=0; i<localNumAtoms; ++i) {
            
            if ( qmAtomGroup[xExt[i].id] > 0 ) {
                
                Real charge = 0;
                
                for (int qi=0; qi<numQMAtms; qi++) {
                    if (qmAtmIndx[qi] == xExt[i].id) {
                        charge = qmAtmChrg[qi];
                        break;
                    }
                }
                
                data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position, 
                                                          fullAtms[i].transform) ;
//                 data_ptr->charge = x[i].charge;
                data_ptr->charge = charge;
                data_ptr->id = xExt[i].id;
                data_ptr->qmGrpID = qmAtomGroup[xExt[i].id] ;
                data_ptr->homeIndx = homeCounter;
                data_ptr->vdwType = fullAtms[i].vdwType;
                
                ++data_ptr;
            }
            else if (meNumMMIndx > 0) {
                // If we have a mechanical embedding simulation with no 
                // point charges, but with QM-MM bonds, we look for the MM1 atoms 
                // connected here and pass them directly.
                
                auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
                
                if (retIt != NULL) {
                    
                    DebugM(3,"Found atom " << retIt->mmIndx << "," << retIt->qmGrp << "\n" );
                    
                    data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position, 
                                                          fullAtms[i].transform) ;
                    // We force the charge to be zero since we would only get
                    // here if mechanical embeding was selected.
                    data_ptr->charge = 0;
                    data_ptr->id = xExt[i].id;
                    data_ptr->qmGrpID = retIt->qmGrp ;
                    data_ptr->homeIndx = homeCounter;
                    
                    // We re-use this varaible to indicate this is an MM1 atom.
                    data_ptr->vdwType = -1;
                    
                    ++data_ptr;
                    
                }
                
            }
                
            homeCounter++;
            
        }
        
    }
    
    if (noPC) {
        for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
            CompAtom *x = (*pdIt).compAtomP;
            (*pdIt).posBoxP->close(&x);
        }
    }
    
    QMProxy[0].recvPartQM(msg);
    
}

void ComputeQM::initialize ( void  )

virtual

Reimplemented from ComputeHomePatches.

Definition at line 615 of file ComputeQM.C.

References SortedArray< Type >::add(), ResizeArray< T >::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< T >::resize(), and Node::simParameters.

{
    ComputeHomePatches::initialize();
    
    // Get some pointers that will be used in the future,
    simParams = Node::Object()->simParameters ;
    molPtr = Node::Object()->molecule;
    // Future comes fast...
    qmAtomGroup = molPtr->get_qmAtomGroup() ;
    
    noPC = molPtr->get_noPC();
    if (noPC) {
        meNumMMIndx = molPtr->get_qmMeNumBonds();
        meMMindx = molPtr->get_qmMeMMindx();
        meQMGrp =molPtr->get_qmMeQMGrp();
        
        for (int i=0; i<meNumMMIndx; i++)
            meQMBonds.add(meMMQMGrp(meMMindx[i], meQMGrp[i]));
    }
    
    numQMAtms = molPtr->get_numQMAtoms();
    qmAtmChrg = molPtr->get_qmAtmChrg() ;
    qmAtmIndx = molPtr->get_qmAtmIndx() ;
    
    numQMGrps = molPtr->get_qmNumGrps();
    qmGrpIDArray = molPtr->get_qmGrpID() ;
    
    cutoff = simParams->cutoff;
    
    customPC = simParams->qmCustomPCSel;
    if (customPC) {
        
        customPCLists.resize(numQMGrps);
        
        int totCustPCs = molPtr->get_qmTotCustPCs();
        const int *custPCSizes = molPtr->get_qmCustPCSizes();
        const int *customPCIdxs = molPtr->get_qmCustomPCIdxs();
        
        int minI = 0, maxI = 0, grpIter = 0;
        
        while (grpIter < numQMGrps) {
            
            maxI += custPCSizes[grpIter];
            
            for (size_t i=minI; i < maxI; i++) {
                
                customPCLists[grpIter].add( customPCIdxs[i] ) ;
            }
            
            // Minimum index gets the size of the current group
            minI += custPCSizes[grpIter];
            // Group iterator gets incremented
            grpIter++;
            
        }
    }
    
}

void ComputeQM::processFullQM

(

QMCoordMsg *

qmFullMsg

)

Definition at line 1292 of file ComputeQM.C.

References ResizeArray< T >::add(), CompAtom::charge, charge, ResizeArray< T >::clear(), QMCoordMsg::coord, DebugM, Lattice::delta(), SortedArray< Type >::find(), CompAtomExt::id, Vector::length(), SortedArray< Type >::load(), FullAtom::mass, QMCoordMsg::numAtoms, QMCoordMsg::numPCIndxs, ComputeHomePatches::patchList, QMCoordMsg::pcIndxs, PCMODECUSTOMSEL, PCMODEUPDATEPOS, PCMODEUPDATESEL, QMCoordMsg::pcSelMode, ComputeQMAtom::position, ComputeQMAtom::qmGrpID, ResizeArray< T >::size(), SortedArray< Type >::sort(), QMCoordMsg::sourceNode, QMPntChrgMsg::sourceNode, CompAtom::vdwType, and x.

Referenced by ComputeQMMgr::recvFullQM().

                                                   {
    
    ResizeArrayIter<PatchElem> ap(patchList); 
    
    const Lattice & lattice = patchList[0].p->lattice;
    
    // Dynamically accumulates point charges as they are found.
    // The same MM atom may be added to this vector more than once if it is 
    // within the cutoff region of two or more QM regions.
    ResizeArray<ComputeQMPntChrg> compPCVec ;
    
    // This will keep the set of QM groups with which each 
    // point charge should be used. It is re-used for each atom in this
    // patch so we can controll the creation of the compPCVec vector.
    // The first item in the pair is the QM Group ID, the second is a pair 
    // with the point charge position (after PBC wraping) and the shortest
    // distance between the point charge and an atom of this QM group.
    GrpDistMap chrgGrpMap ;
    
    DebugM(4,"Rank " << CkMyPe() << " receiving from rank " << qmFullMsg->sourceNode
    << " a total of " << qmFullMsg->numAtoms << " QM atoms and " 
    << qmFullMsg->numPCIndxs << " PC IDs." << std::endl);
    
    // If this is the firts step after a step of PC re-selection,
    // we store all PC IDs that all PEs found, in case some atoms end up
    // here untill the next re-selection step.
    if (qmFullMsg->numPCIndxs) {
        
        pcIDSortList.clear();
        
        int *pcIndx = qmFullMsg->pcIndxs;
        for (int i=0; i< qmFullMsg->numPCIndxs;i++) {
            pcIDSortList.load(pcIndx[i]);
        }
        
        pcIDSortList.sort();
    }
    
    int totalNodeAtoms = 0;
    int atomIter = 0;
    int uniqueCounter = 0;
    
    const ComputeQMAtom *qmDataPtn = qmFullMsg->coord;
    
    switch ( qmFullMsg->pcSelMode ) {
    
    case PCMODEUPDATESEL:
    {
        
    DebugM(4,"Updating PC selection.\n")
    
    // Loops over all atoms in this node and checks if any MM atom is within 
    // the cutof radius from a QM atom
    for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
        
        CompAtom *x = (*pdIt).compAtomP;
        const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
        int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
        const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
        
        totalNodeAtoms += localNumAtoms;
        
        // Iterates over the local atoms in a PE, all patches.
        for(int i=0; i<localNumAtoms; ++i) {
            
            // A new subset for each atom in this node.
            chrgGrpMap.clear();
            
            Real pcGrpID = qmAtomGroup[xExt[i].id];
            
            Real charge = x[i].charge;
            
            // If the atom is a QM atom, there will be no charge info in the 
            // atom box. Therefore, we take the charge in the previous step
            // in case this atom is a point charge for a neighboring QM region.
            if (pcGrpID > 0) {
                for (int qi=0; qi<numQMAtms; qi++) {
                    if (qmAtmIndx[qi] == xExt[i].id) {
                        charge = qmAtmChrg[qi];
                        break;
                    }
                }
            }
            
            qmDataPtn = qmFullMsg->coord;
            for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
                
                Real qmGrpID = qmDataPtn->qmGrpID;
                
                // We check If a QM atom and the node atom "i"
                // belong to the same QM group. If not, atom "i" is either an MM
                // atom or a QM atom from another group, which will be seen
                // as an MM point charge.
                // If they belong to the same group, just skip the distance 
                // calculation and move on to the next QM atom.
                // This loop needs to go over all QM atoms since a point charge
                // may me be close to two different QM groups, in which case it 
                // will be sent to both.
                if (qmGrpID == pcGrpID) {
                    continue;
                }
                
                // calculate distance between QM atom and Poitn Charge
                Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
                
                BigReal dist = r12.length();  // Distance between atoms
                
                if ( dist < cutoff ){
                    
                    // Since the QM region may have been wrapped around the periodic box, we
                    // reposition point charges with respect to unwrapped coordinates of QM atoms,
                    // finding the shortest distance between the point charge and the QM region.
                    
                    Position posMM = qmDataPtn->position + r12;
                    auto ret = chrgGrpMap.find(qmGrpID) ;
                    
                    // If 'ret' points to end, it means the item was not found
                    // and will be added, which means a new QM region has this 
                    // atom within its cutoff region,
                    // which means a new point charge will be 
                    // created in the QM system in question.
                    // 'ret' means a lot to me!
                    if ( ret ==  chrgGrpMap.end()) {
                        chrgGrpMap.insert(std::pair<Real,PosDistPair>(qmGrpID, 
                                                         PosDistPair(posMM,dist) ) );
                    }
                    else {
                        // If the current distance is smaller than the one in 
                        // the map, we update it so that we have the smallest
                        // distance between the point charge and any QM atom in 
                        // this group.
                        if (dist < ret->second.second) {
                            ret->second.first = posMM ;
                            ret->second.second = dist ;
                        }
                    }
                }
            }
            
            for(auto mapIt = chrgGrpMap.begin(); 
                mapIt != chrgGrpMap.end(); mapIt++) {
                
                // We now add the final info about this point charge 
                // to the vector, repeating it for each QM group that has it
                // within its cuttoff radius.
                compPCVec.add(
                    ComputeQMPntChrg(mapIt->second.first, charge, xExt[i].id,
                                  mapIt->first, atomIter, mapIt->second.second, 
                                  fullAtms[i].mass, fullAtms[i].vdwType)
                               );
            }
            
            // Counts how many atoms are seens as point charges, by one or more
            // QM groups.
            if (chrgGrpMap.size() > 0)
                uniqueCounter++;
            
            atomIter++;
        }
        
        (*pdIt).posBoxP->close(&x);
    }
    
    }break; // End case PCMODEUPDATESEL
    
    case PCMODEUPDATEPOS: 
    {
    
    DebugM(4,"Updating PC positions ONLY.\n")
    
    for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
        
        CompAtom *x = (*pdIt).compAtomP;
        const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
        int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
        const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
        
        totalNodeAtoms += localNumAtoms;
        
        // Iterates over the local atoms in a PE, all patches.
        for(int i=0; i<localNumAtoms; ++i) {
            
            if (pcIDSortList.find(xExt[i].id) != NULL ) {
                
                BigReal dist =  0;
                Position posMM = 0; // Temporary point charge position for QM calculation.
                
                bool firstIngroupQM = true;
                qmDataPtn = qmFullMsg->coord;
                for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
                    
                    // Assuming we only have one QM region for now
                    // This optio will be deprecated
                    
                    // Only verify distances to atoms in the QM group for which we are
                    // gathering point charges.
//                     if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
                    
                    // calculate distance between QM atom and Poitn Charge
                    Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
                    
                    // We wait to find the first QM atom in the target QM group before we
                    // set the first guess at a minimal distance and point charge position.
                    if ( firstIngroupQM ) {
                        dist = r12.length();
                        
                        // Reposition classical point charge with respect to unwrapped QM position.
                        posMM = qmDataPtn->position + r12;
                        
                        firstIngroupQM = false;
                    }
                    
                    // If we find a QM atom with a shorter dstance to the point charge,
                    // set that as the PC distance and reposition the PC
                    if ( r12.length() < dist ) {
                        dist = r12.length();
                        posMM = qmDataPtn->position + r12 ;
                    }
                    
                }
                
                Real pcGrpID = qmAtomGroup[xExt[i].id];
                Real charge = x[i].charge;
                
                // If the atom is a QM atom, there will be no charge info in the 
                // atom box. Therefore, we take the charge in the previous step
                // in case this atom is a point charge for a neighboring QM region.
                if (pcGrpID > 0) {
                    for (int qi=0; qi<numQMAtms; qi++) {
                        if (qmAtmIndx[qi] == xExt[i].id) {
                            charge = qmAtmChrg[qi];
                            break;
                        }
                        
                    }
                }
                
                compPCVec.add(
                    ComputeQMPntChrg(posMM, charge, xExt[i].id,
                                  0, atomIter, 0, 
                                  fullAtms[i].mass, fullAtms[i].vdwType));
            }
            
            atomIter++;
        }
        
        (*pdIt).posBoxP->close(&x);
    }
    }break ; // End case PCMODEUPDATEPOS
    
    case PCMODECUSTOMSEL:
    {
    
    DebugM(4,"Updating PC positions for custom PC selection.\n")
    
    for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
        
        CompAtom *x = (*pdIt).compAtomP;
        const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
        int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
        const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
        
        totalNodeAtoms += localNumAtoms;
        
        // Iterates over the local atoms in a PE, all patches.
        for(int i=0; i<localNumAtoms; ++i) {
            
            // Checks if the atom ID belongs to a point charge list,
            // which indicates it is to be sent to rank zero for QM calculations.
            // With one list per QM group, the same point charge can be added to 
            // different QM groups, as if it was being dynamically selected.
            for (int grpIndx=0; grpIndx<numQMGrps; grpIndx++) {
                
                if (customPCLists[grpIndx].find(xExt[i].id) != NULL){
                    
                    // Initialize the search for the smallest distance between
                    // point charge and QM atoms. This will determine the closest position
                    // of the point charge accounting for periodic boundary conditions.
                    
                    // Since the QM region may have been wrapped around the periodic box, we
                    // reposition point charges with respect to unwrapped coordinates of QM atoms.
                    
                    BigReal dist =  0;
                    Position posMM = 0; // Temporary point charge position for QM calculation.
                    
                    bool firstIngroupQM = true;
                    qmDataPtn = qmFullMsg->coord;
                    for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
                        
                        // Only verify distances to atoms in the QM group for which we are
                        // gathering point charges.
                        if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
                        
                        // calculate distance between QM atom and Poitn Charge
                        Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
                        
                        // We wait to find the first QM atom in the target QM group before we
                        // set the first guess at a minimal distance and point charge position.
                        if ( firstIngroupQM ) {
                            dist = r12.length();
                            
                            // Reposition classical point charge with respect to unwrapped QM position.
                            posMM = qmDataPtn->position + r12;
                            
                            firstIngroupQM = false;
                        }
                        
                        // If we find a QM atom with a shorter dstance to the point charge,
                        // set that as the PC distance and reposition the PC
                        if ( r12.length() < dist ) {
                            dist = r12.length();
                            posMM = qmDataPtn->position + r12 ;
                        }
                        
                    }
                    
                    // After determining the position o fthe point charge, get its partial charge.
                    Real pcGrpID = qmAtomGroup[xExt[i].id];
                    Real charge = x[i].charge;
                    
                    // If the atom is a QM atom, there will be no charge info in the 
                    // atom box. Therefore, we take the charge in the previous step
                    // in case this atom is a point charge for a neighboring QM region.
                    if (pcGrpID > 0) {
                        for (int qi=0; qi<numQMAtms; qi++) {
                            if (qmAtmIndx[qi] == xExt[i].id) {
                                charge = qmAtmChrg[qi];
                                break;
                            }
                        }
                    }
                    
                    compPCVec.add(
                        ComputeQMPntChrg(posMM, charge, xExt[i].id,
                                      qmGrpIDArray[grpIndx], atomIter, dist, 
                                      fullAtms[i].mass, fullAtms[i].vdwType));
                }
                
            }
            
            atomIter++;
        }
        
        (*pdIt).posBoxP->close(&x);
    }
    
    } break;
    }
    
    DebugM(4,"Rank " << CkMyPe() << " found a total of " << compPCVec.size()
    << " point charges, out of " << totalNodeAtoms
    << " atoms in this node. " << uniqueCounter << " are unique." << std::endl);
    
    // Send only the MM atoms within radius
    
    // this pointer should be freed in rank zero, after receiving it.
    QMPntChrgMsg *pntChrgMsg = new (compPCVec.size(), 0) QMPntChrgMsg;
    pntChrgMsg->sourceNode = CkMyPe();
    pntChrgMsg->numAtoms = compPCVec.size();
    
    for (int i=0; i<compPCVec.size(); i++ ) {
        
        // Only sends the positions and charges of atoms within
        // cutoff of a (any) QM atom. Repeats for the number of QM groups
        // this charge is near to, and indicates the QM group it should 
        // be used with.
        pntChrgMsg->coord[i] = compPCVec[i] ;
        
    }
    
    DebugM(4,"Rank " << pntChrgMsg->sourceNode << " sending a total of " 
    << compPCVec.size() << " elements to rank zero." << std::endl);
    
    CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
    QMProxy[0].recvPntChrg(pntChrgMsg);
    
    delete qmFullMsg;
}

void ComputeQM::saveResults

(

QMForceMsg *

fmsg

)

Definition at line 2673 of file ComputeQM.C.

References ResizeArrayIter< Type >::begin(), QMForce::charge, DebugM, ResizeArrayIter< Type >::end(), QMForceMsg::energy, Results::f, QMForce::force, QMForceMsg::force, ExtForce::force, Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), Molecule::get_qmNumGrps(), getComputes(), QMForce::homeIndx, QMForce::id, SubmitReduction::item(), Node::molecule, Results::normal, QMForceMsg::numForces, Node::Object(), ComputeHomePatches::patchList, QMForceMsg::PMEOn, REDUCTION_MISC_ENERGY, QMForce::replace, ExtForce::replace, ResizeArray< T >::size(), SubmitReduction::submit(), and QMForceMsg::virial.

Referenced by ComputeQMMgr::recvForce().

                                            {
    
    ResizeArrayIter<PatchElem> ap(patchList);
    
    bool callReplaceForces = false;
    
    int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
    int numQMGrps = Node::Object()->molecule->get_qmNumGrps();
    const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
    const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
    Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
    
    int totalNumbAtoms = 0;
    for (ap = ap.begin(); ap != ap.end(); ap++) {
        totalNumbAtoms += (*ap).p->getNumAtoms();
    }
    
    // This is kept to be deleted in the next step so that the message can be 
    // used in "replaceForces" routine later on, in case there are forces to 
    // be replaced. The "replaceForces" function uses but DOES NOT free the pointer,
    // so we free the data from the previous iteration and allocate a new one for
    // the current iteration (remember that the number of atoms can change in a 
    // home patch between iterations).
    if (oldForces != 0)
        delete [] oldForces;
    oldForces = new ExtForce[totalNumbAtoms] ;
    
    for (int i=0; i < totalNumbAtoms; ++i) {
        oldForces[i].force = Force(0) ;
    }
    
    DebugM(1,"Force array has been created and zeroed in rank " 
    << CkMyPe() << std::endl);
    
    DebugM(1,"Preparing " << fmsg->numForces << " forces in rank " 
    << CkMyPe() << std::endl);
    
    QMForce *results_ptr = fmsg->force;
    // Iterates over the received forces and place them on the full array.
    for (int i=0; i < fmsg->numForces; ++i, ++results_ptr) {
        // For now we may have more than one item in the message acting on the same 
        // atom, such as an MM atom which is a point charge for two or more QM regions.
        
        oldForces[results_ptr->homeIndx].force += results_ptr->force;
        oldForces[results_ptr->homeIndx].replace = results_ptr->replace;
        
        if (results_ptr->replace == 1)
            callReplaceForces = true;
        
        // If the atom is in a QM group, update its charge to the local (this homePatch)
        // copy of the qmAtmChrg array.
        if (qmAtomGroup[results_ptr->id] > 0 && (fmsg->PMEOn || (numQMGrps > 1) ) ) {
            
            // Loops over all QM atoms (in all QM groups) comparing their global indices
            for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
                
                if (qmAtmIndx[qmIter] == results_ptr->id) {
                    qmAtmChrg[qmIter] = results_ptr->charge;
                    break;
                }
                
            }
            
        }
        
    }
    
    DebugM(1,"Placing forces on force boxes in rank " 
    << CkMyPe() << std::endl);
    
    // Places the received forces on the force array for each patch.
    int homeIndxIter = 0;
    for (ap = ap.begin(); ap != ap.end(); ap++) {
        Results *r = (*ap).forceBox->open();
        Force *f = r->f[Results::normal];
        const CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
        int localNumAtoms = (*ap).p->getNumAtoms();
        
        for(int i=0; i<localNumAtoms; ++i) {
            
            f[i] += oldForces[homeIndxIter].force;
            
            ++homeIndxIter;
        }
        
        if ( callReplaceForces )
            (*ap).p->replaceForces(oldForces);
        
        (*ap).forceBox->close(&r);
        
    }
    
    DebugM(1,"Forces placed on force boxes in rank " 
    << CkMyPe() << std::endl);
    
    if (fmsg->PMEOn) {
        
        DebugM(1,"PME ON! Accessing PMEmgr in rank " << CkMyPe() << std::endl);
        
        ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch(
            CkpvAccess(BOCclass_group).computePmeMgr) ;
        
        DebugM(4, "Initiating ComputePme::doQMWork on rank " << CkMyPe() << " over " 
            << getComputes(mgrP).size() << " pmeComputes." << std::endl) ;
        
        for ( int i=0; i< getComputes(mgrP).size(); ++i ) {
    //         WorkDistrib::messageEnqueueWork(pmeComputes[i]);
            getComputes(mgrP)[i]->doQMWork();
        }
    }
    
    DebugM(1,"Submitting reduction in rank " << CkMyPe() << std::endl);
    
    reduction->item(REDUCTION_MISC_ENERGY) += fmsg->energy;
    reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += fmsg->virial[0][0];
    reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += fmsg->virial[0][1];
    reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += fmsg->virial[0][2];
    reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += fmsg->virial[1][0];
    reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += fmsg->virial[1][1];
    reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += fmsg->virial[1][2];
    reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += fmsg->virial[2][0];
    reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += fmsg->virial[2][1];
    reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += fmsg->virial[2][2];
    reduction->submit();
    
    delete fmsg ;
}

---

The documentation for this class was generated from the following files:

• ComputeQM.h
• ComputeQM.C