NAMD
Public Member Functions | List of all members
ComputeQM Class Reference

#include <ComputeQM.h>

Inheritance diagram for ComputeQM:
ComputeHomePatches Compute

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

599  :
600  ComputeHomePatches(c), oldForces(0)
601 {
602  CProxy_ComputeQMMgr::ckLocalBranch(
603  CkpvAccess(BOCclass_group).computeQMMgr)->setCompute(this);
604 
606 
607 }
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
ComputeHomePatches(ComputeID c)
ComputeQM::~ComputeQM ( )
virtual

Definition at line 609 of file ComputeQM.C.

610 {
611  if (oldForces != 0)
612  delete [] oldForces;
613 }

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.

676 {
677  CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
678 
680 
681  int timeStep;
682 
683  #ifdef DEBUG_QM
684  DebugM(4,"----> Initiating QM work on rank " << CkMyPe() <<
685  " with " << patchList.size() << " patches." << std::endl );
686  #endif
687 
688  patchData.clear();
689 
690  // Determines hou many QM atoms are in this node.
691  int numLocalQMAtoms = 0, localMM1atoms = 0;
692  for (ap = ap.begin(); ap != ap.end(); ap++) {
693  CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
694  int localNumAtoms = (*ap).p->getNumAtoms();
695 
696  patchData.push_back(patchDataStrc((*ap).positionBox,
697  (*ap).positionBox->open(),
698  (*ap).p) );
699 
700  for(int i=0; i<localNumAtoms; ++i) {
701  if ( qmAtomGroup[xExt[i].id] > 0 ) {
702  numLocalQMAtoms++;
703  }
704  else if (meNumMMIndx > 0) {
705  // If we have a mechanical embedding simulation with no
706  // point charges, but with QM-MM bonds, we look for the MM1 atoms
707  // here and pass them directly. No need for an extra round of
708  // comms just to get atoms we already know we need.
709 
710  auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
711 
712  if (retIt != NULL) {
713  localMM1atoms++;
714  }
715  }
716  }
717 
718  timeStep = (*ap).p->flags.step ;
719  }
720 
721  DebugM(4, "Node " << CkMyPe() << " has " << numLocalQMAtoms
722  << " QM atoms." << std::endl) ;
723  #ifdef DEBUG_QM
724  if (localMM1atoms > 0)
725  DebugM(4, "Node " << CkMyPe() << " has " << localMM1atoms
726  << " MM1 atoms." << std::endl) ;
727  #endif
728  // Send QM atoms
729 
730  // This pointer will be freed in rank zero.
731  QMCoordMsg *msg = new (numLocalQMAtoms + localMM1atoms,0, 0) QMCoordMsg;
732  msg->sourceNode = CkMyPe();
733  msg->numAtoms = numLocalQMAtoms + localMM1atoms;
734  msg->timestep = timeStep;
735  ComputeQMAtom *data_ptr = msg->coord;
736 
737  // Build a message with coordinates, charge, atom ID and QM group
738  // for all QM atoms in the node, then send it to node zero.
739  int homeCounter = 0;
740  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
741 
742  CompAtom *x = (*pdIt).compAtomP;
743  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
744  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
745  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
746 
747  for(int i=0; i<localNumAtoms; ++i) {
748 
749  if ( qmAtomGroup[xExt[i].id] > 0 ) {
750 
751  Real charge = 0;
752 
753  for (int qi=0; qi<numQMAtms; qi++) {
754  if (qmAtmIndx[qi] == xExt[i].id) {
755  charge = qmAtmChrg[qi];
756  break;
757  }
758  }
759 
760  data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position,
761  fullAtms[i].transform) ;
762 // data_ptr->charge = x[i].charge;
763  data_ptr->charge = charge;
764  data_ptr->id = xExt[i].id;
765  data_ptr->qmGrpID = qmAtomGroup[xExt[i].id] ;
766  data_ptr->homeIndx = homeCounter;
767  data_ptr->vdwType = fullAtms[i].vdwType;
768 
769  ++data_ptr;
770  }
771  else if (meNumMMIndx > 0) {
772  // If we have a mechanical embedding simulation with no
773  // point charges, but with QM-MM bonds, we look for the MM1 atoms
774  // connected here and pass them directly.
775 
776  auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
777 
778  if (retIt != NULL) {
779 
780  DebugM(3,"Found atom " << retIt->mmIndx << "," << retIt->qmGrp << "\n" );
781 
782  data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position,
783  fullAtms[i].transform) ;
784  // We force the charge to be zero since we would only get
785  // here if mechanical embeding was selected.
786  data_ptr->charge = 0;
787  data_ptr->id = xExt[i].id;
788  data_ptr->qmGrpID = retIt->qmGrp ;
789  data_ptr->homeIndx = homeCounter;
790 
791  // We re-use this varaible to indicate this is an MM1 atom.
792  data_ptr->vdwType = -1;
793 
794  ++data_ptr;
795 
796  }
797 
798  }
799 
800  homeCounter++;
801 
802  }
803 
804  }
805 
806  if (noPC) {
807  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
808  CompAtom *x = (*pdIt).compAtomP;
809  (*pdIt).posBoxP->close(&x);
810  }
811  }
812 
813  QMProxy[0].recvPartQM(msg);
814 
815 }
int timestep
Definition: ComputeQM.C:107
ComputeHomePatchList patchList
float charge
Definition: ComputeQM.C:67
float Real
Definition: common.h:109
#define DebugM(x, y)
Definition: Debug.h:59
ComputeQMAtom * coord
Definition: ComputeQM.C:110
Position position
Definition: ComputeQM.C:66
int sourceNode
Definition: ComputeQM.C:105
Real qmGrpID
Definition: ComputeQM.C:69
short vdwType
Definition: NamdTypes.h:55
Elem * find(const Elem &elem)
int numAtoms
Definition: ComputeQM.C:106
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
int size(void) const
Definition: ResizeArray.h:127
gridSize x
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.

616 {
618 
619  // Get some pointers that will be used in the future,
620  simParams = Node::Object()->simParameters ;
621  molPtr = Node::Object()->molecule;
622  // Future comes fast...
623  qmAtomGroup = molPtr->get_qmAtomGroup() ;
624 
625  noPC = molPtr->get_noPC();
626  if (noPC) {
627  meNumMMIndx = molPtr->get_qmMeNumBonds();
628  meMMindx = molPtr->get_qmMeMMindx();
629  meQMGrp =molPtr->get_qmMeQMGrp();
630 
631  for (int i=0; i<meNumMMIndx; i++)
632  meQMBonds.add(meMMQMGrp(meMMindx[i], meQMGrp[i]));
633  }
634 
635  numQMAtms = molPtr->get_numQMAtoms();
636  qmAtmChrg = molPtr->get_qmAtmChrg() ;
637  qmAtmIndx = molPtr->get_qmAtmIndx() ;
638 
639  numQMGrps = molPtr->get_qmNumGrps();
640  qmGrpIDArray = molPtr->get_qmGrpID() ;
641 
642  cutoff = simParams->cutoff;
643 
644  customPC = simParams->qmCustomPCSel;
645  if (customPC) {
646 
647  customPCLists.resize(numQMGrps);
648 
649  int totCustPCs = molPtr->get_qmTotCustPCs();
650  const int *custPCSizes = molPtr->get_qmCustPCSizes();
651  const int *customPCIdxs = molPtr->get_qmCustomPCIdxs();
652 
653  int minI = 0, maxI = 0, grpIter = 0;
654 
655  while (grpIter < numQMGrps) {
656 
657  maxI += custPCSizes[grpIter];
658 
659  for (size_t i=minI; i < maxI; i++) {
660 
661  customPCLists[grpIter].add( customPCIdxs[i] ) ;
662  }
663 
664  // Minimum index gets the size of the current group
665  minI += custPCSizes[grpIter];
666  // Group iterator gets incremented
667  grpIter++;
668 
669  }
670  }
671 
672 }
static Node * Object()
Definition: Node.h:86
Bool get_noPC()
Definition: Molecule.h:856
int get_qmNumGrps() const
Definition: Molecule.h:827
const int * get_qmAtmIndx()
Definition: Molecule.h:823
int add(const Elem &elem)
Definition: SortedArray.h:55
SimParameters * simParameters
Definition: Node.h:178
int get_numQMAtoms()
Definition: Molecule.h:825
int get_qmTotCustPCs()
Definition: Molecule.h:863
Real * get_qmMeQMGrp()
Definition: Molecule.h:859
const Real * get_qmAtomGroup() const
Definition: Molecule.h:819
int get_qmMeNumBonds()
Definition: Molecule.h:857
Real * get_qmGrpID()
Definition: Molecule.h:829
int * get_qmMeMMindx()
Definition: Molecule.h:858
Real * get_qmAtmChrg()
Definition: Molecule.h:822
int * get_qmCustPCSizes()
Definition: Molecule.h:864
int add(const Elem &elem)
Definition: ResizeArray.h:97
void resize(int i)
Definition: ResizeArray.h:84
virtual void initialize()
int * get_qmCustomPCIdxs()
Definition: Molecule.h:865
Molecule * molecule
Definition: Node.h:176
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().

1292  {
1293 
1295 
1296  const Lattice & lattice = patchList[0].p->lattice;
1297 
1298  // Dynamically accumulates point charges as they are found.
1299  // The same MM atom may be added to this vector more than once if it is
1300  // within the cutoff region of two or more QM regions.
1301  ResizeArray<ComputeQMPntChrg> compPCVec ;
1302 
1303  // This will keep the set of QM groups with which each
1304  // point charge should be used. It is re-used for each atom in this
1305  // patch so we can controll the creation of the compPCVec vector.
1306  // The first item in the pair is the QM Group ID, the second is a pair
1307  // with the point charge position (after PBC wraping) and the shortest
1308  // distance between the point charge and an atom of this QM group.
1309  GrpDistMap chrgGrpMap ;
1310 
1311  DebugM(4,"Rank " << CkMyPe() << " receiving from rank " << qmFullMsg->sourceNode
1312  << " a total of " << qmFullMsg->numAtoms << " QM atoms and "
1313  << qmFullMsg->numPCIndxs << " PC IDs." << std::endl);
1314 
1315  // If this is the firts step after a step of PC re-selection,
1316  // we store all PC IDs that all PEs found, in case some atoms end up
1317  // here untill the next re-selection step.
1318  if (qmFullMsg->numPCIndxs) {
1319 
1320  pcIDSortList.clear();
1321 
1322  int *pcIndx = qmFullMsg->pcIndxs;
1323  for (int i=0; i< qmFullMsg->numPCIndxs;i++) {
1324  pcIDSortList.load(pcIndx[i]);
1325  }
1326 
1327  pcIDSortList.sort();
1328  }
1329 
1330  int totalNodeAtoms = 0;
1331  int atomIter = 0;
1332  int uniqueCounter = 0;
1333 
1334  const ComputeQMAtom *qmDataPtn = qmFullMsg->coord;
1335 
1336  switch ( qmFullMsg->pcSelMode ) {
1337 
1338  case PCMODEUPDATESEL:
1339  {
1340 
1341  DebugM(4,"Updating PC selection.\n")
1342 
1343  // Loops over all atoms in this node and checks if any MM atom is within
1344  // the cutof radius from a QM atom
1345  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
1346 
1347  CompAtom *x = (*pdIt).compAtomP;
1348  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
1349  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
1350  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
1351 
1352  totalNodeAtoms += localNumAtoms;
1353 
1354  // Iterates over the local atoms in a PE, all patches.
1355  for(int i=0; i<localNumAtoms; ++i) {
1356 
1357  // A new subset for each atom in this node.
1358  chrgGrpMap.clear();
1359 
1360  Real pcGrpID = qmAtomGroup[xExt[i].id];
1361 
1362  Real charge = x[i].charge;
1363 
1364  // If the atom is a QM atom, there will be no charge info in the
1365  // atom box. Therefore, we take the charge in the previous step
1366  // in case this atom is a point charge for a neighboring QM region.
1367  if (pcGrpID > 0) {
1368  for (int qi=0; qi<numQMAtms; qi++) {
1369  if (qmAtmIndx[qi] == xExt[i].id) {
1370  charge = qmAtmChrg[qi];
1371  break;
1372  }
1373  }
1374  }
1375 
1376  qmDataPtn = qmFullMsg->coord;
1377  for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
1378 
1379  Real qmGrpID = qmDataPtn->qmGrpID;
1380 
1381  // We check If a QM atom and the node atom "i"
1382  // belong to the same QM group. If not, atom "i" is either an MM
1383  // atom or a QM atom from another group, which will be seen
1384  // as an MM point charge.
1385  // If they belong to the same group, just skip the distance
1386  // calculation and move on to the next QM atom.
1387  // This loop needs to go over all QM atoms since a point charge
1388  // may me be close to two different QM groups, in which case it
1389  // will be sent to both.
1390  if (qmGrpID == pcGrpID) {
1391  continue;
1392  }
1393 
1394  // calculate distance between QM atom and Poitn Charge
1395  Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
1396 
1397  BigReal dist = r12.length(); // Distance between atoms
1398 
1399  if ( dist < cutoff ){
1400 
1401  // Since the QM region may have been wrapped around the periodic box, we
1402  // reposition point charges with respect to unwrapped coordinates of QM atoms,
1403  // finding the shortest distance between the point charge and the QM region.
1404 
1405  Position posMM = qmDataPtn->position + r12;
1406  auto ret = chrgGrpMap.find(qmGrpID) ;
1407 
1408  // If 'ret' points to end, it means the item was not found
1409  // and will be added, which means a new QM region has this
1410  // atom within its cutoff region,
1411  // which means a new point charge will be
1412  // created in the QM system in question.
1413  // 'ret' means a lot to me!
1414  if ( ret == chrgGrpMap.end()) {
1415  chrgGrpMap.insert(std::pair<Real,PosDistPair>(qmGrpID,
1416  PosDistPair(posMM,dist) ) );
1417  }
1418  else {
1419  // If the current distance is smaller than the one in
1420  // the map, we update it so that we have the smallest
1421  // distance between the point charge and any QM atom in
1422  // this group.
1423  if (dist < ret->second.second) {
1424  ret->second.first = posMM ;
1425  ret->second.second = dist ;
1426  }
1427  }
1428  }
1429  }
1430 
1431  for(auto mapIt = chrgGrpMap.begin();
1432  mapIt != chrgGrpMap.end(); mapIt++) {
1433 
1434  // We now add the final info about this point charge
1435  // to the vector, repeating it for each QM group that has it
1436  // within its cuttoff radius.
1437  compPCVec.add(
1438  ComputeQMPntChrg(mapIt->second.first, charge, xExt[i].id,
1439  mapIt->first, atomIter, mapIt->second.second,
1440  fullAtms[i].mass, fullAtms[i].vdwType)
1441  );
1442  }
1443 
1444  // Counts how many atoms are seens as point charges, by one or more
1445  // QM groups.
1446  if (chrgGrpMap.size() > 0)
1447  uniqueCounter++;
1448 
1449  atomIter++;
1450  }
1451 
1452  (*pdIt).posBoxP->close(&x);
1453  }
1454 
1455  }break; // End case PCMODEUPDATESEL
1456 
1457  case PCMODEUPDATEPOS:
1458  {
1459 
1460  DebugM(4,"Updating PC positions ONLY.\n")
1461 
1462  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
1463 
1464  CompAtom *x = (*pdIt).compAtomP;
1465  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
1466  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
1467  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
1468 
1469  totalNodeAtoms += localNumAtoms;
1470 
1471  // Iterates over the local atoms in a PE, all patches.
1472  for(int i=0; i<localNumAtoms; ++i) {
1473 
1474  if (pcIDSortList.find(xExt[i].id) != NULL ) {
1475 
1476  BigReal dist = 0;
1477  Position posMM = 0; // Temporary point charge position for QM calculation.
1478 
1479  bool firstIngroupQM = true;
1480  qmDataPtn = qmFullMsg->coord;
1481  for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
1482 
1483  // Assuming we only have one QM region for now
1484  // This optio will be deprecated
1485 
1486  // Only verify distances to atoms in the QM group for which we are
1487  // gathering point charges.
1488 // if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
1489 
1490  // calculate distance between QM atom and Poitn Charge
1491  Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
1492 
1493  // We wait to find the first QM atom in the target QM group before we
1494  // set the first guess at a minimal distance and point charge position.
1495  if ( firstIngroupQM ) {
1496  dist = r12.length();
1497 
1498  // Reposition classical point charge with respect to unwrapped QM position.
1499  posMM = qmDataPtn->position + r12;
1500 
1501  firstIngroupQM = false;
1502  }
1503 
1504  // If we find a QM atom with a shorter dstance to the point charge,
1505  // set that as the PC distance and reposition the PC
1506  if ( r12.length() < dist ) {
1507  dist = r12.length();
1508  posMM = qmDataPtn->position + r12 ;
1509  }
1510 
1511  }
1512 
1513  Real pcGrpID = qmAtomGroup[xExt[i].id];
1514  Real charge = x[i].charge;
1515 
1516  // If the atom is a QM atom, there will be no charge info in the
1517  // atom box. Therefore, we take the charge in the previous step
1518  // in case this atom is a point charge for a neighboring QM region.
1519  if (pcGrpID > 0) {
1520  for (int qi=0; qi<numQMAtms; qi++) {
1521  if (qmAtmIndx[qi] == xExt[i].id) {
1522  charge = qmAtmChrg[qi];
1523  break;
1524  }
1525 
1526  }
1527  }
1528 
1529  compPCVec.add(
1530  ComputeQMPntChrg(posMM, charge, xExt[i].id,
1531  0, atomIter, 0,
1532  fullAtms[i].mass, fullAtms[i].vdwType));
1533  }
1534 
1535  atomIter++;
1536  }
1537 
1538  (*pdIt).posBoxP->close(&x);
1539  }
1540  }break ; // End case PCMODEUPDATEPOS
1541 
1542  case PCMODECUSTOMSEL:
1543  {
1544 
1545  DebugM(4,"Updating PC positions for custom PC selection.\n")
1546 
1547  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
1548 
1549  CompAtom *x = (*pdIt).compAtomP;
1550  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
1551  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
1552  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
1553 
1554  totalNodeAtoms += localNumAtoms;
1555 
1556  // Iterates over the local atoms in a PE, all patches.
1557  for(int i=0; i<localNumAtoms; ++i) {
1558 
1559  // Checks if the atom ID belongs to a point charge list,
1560  // which indicates it is to be sent to rank zero for QM calculations.
1561  // With one list per QM group, the same point charge can be added to
1562  // different QM groups, as if it was being dynamically selected.
1563  for (int grpIndx=0; grpIndx<numQMGrps; grpIndx++) {
1564 
1565  if (customPCLists[grpIndx].find(xExt[i].id) != NULL){
1566 
1567  // Initialize the search for the smallest distance between
1568  // point charge and QM atoms. This will determine the closest position
1569  // of the point charge accounting for periodic boundary conditions.
1570 
1571  // Since the QM region may have been wrapped around the periodic box, we
1572  // reposition point charges with respect to unwrapped coordinates of QM atoms.
1573 
1574  BigReal dist = 0;
1575  Position posMM = 0; // Temporary point charge position for QM calculation.
1576 
1577  bool firstIngroupQM = true;
1578  qmDataPtn = qmFullMsg->coord;
1579  for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
1580 
1581  // Only verify distances to atoms in the QM group for which we are
1582  // gathering point charges.
1583  if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
1584 
1585  // calculate distance between QM atom and Poitn Charge
1586  Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
1587 
1588  // We wait to find the first QM atom in the target QM group before we
1589  // set the first guess at a minimal distance and point charge position.
1590  if ( firstIngroupQM ) {
1591  dist = r12.length();
1592 
1593  // Reposition classical point charge with respect to unwrapped QM position.
1594  posMM = qmDataPtn->position + r12;
1595 
1596  firstIngroupQM = false;
1597  }
1598 
1599  // If we find a QM atom with a shorter dstance to the point charge,
1600  // set that as the PC distance and reposition the PC
1601  if ( r12.length() < dist ) {
1602  dist = r12.length();
1603  posMM = qmDataPtn->position + r12 ;
1604  }
1605 
1606  }
1607 
1608  // After determining the position o fthe point charge, get its partial charge.
1609  Real pcGrpID = qmAtomGroup[xExt[i].id];
1610  Real charge = x[i].charge;
1611 
1612  // If the atom is a QM atom, there will be no charge info in the
1613  // atom box. Therefore, we take the charge in the previous step
1614  // in case this atom is a point charge for a neighboring QM region.
1615  if (pcGrpID > 0) {
1616  for (int qi=0; qi<numQMAtms; qi++) {
1617  if (qmAtmIndx[qi] == xExt[i].id) {
1618  charge = qmAtmChrg[qi];
1619  break;
1620  }
1621  }
1622  }
1623 
1624  compPCVec.add(
1625  ComputeQMPntChrg(posMM, charge, xExt[i].id,
1626  qmGrpIDArray[grpIndx], atomIter, dist,
1627  fullAtms[i].mass, fullAtms[i].vdwType));
1628  }
1629 
1630  }
1631 
1632  atomIter++;
1633  }
1634 
1635  (*pdIt).posBoxP->close(&x);
1636  }
1637 
1638  } break;
1639  }
1640 
1641  DebugM(4,"Rank " << CkMyPe() << " found a total of " << compPCVec.size()
1642  << " point charges, out of " << totalNodeAtoms
1643  << " atoms in this node. " << uniqueCounter << " are unique." << std::endl);
1644 
1645  // Send only the MM atoms within radius
1646 
1647  // this pointer should be freed in rank zero, after receiving it.
1648  QMPntChrgMsg *pntChrgMsg = new (compPCVec.size(), 0) QMPntChrgMsg;
1649  pntChrgMsg->sourceNode = CkMyPe();
1650  pntChrgMsg->numAtoms = compPCVec.size();
1651 
1652  for (int i=0; i<compPCVec.size(); i++ ) {
1653 
1654  // Only sends the positions and charges of atoms within
1655  // cutoff of a (any) QM atom. Repeats for the number of QM groups
1656  // this charge is near to, and indicates the QM group it should
1657  // be used with.
1658  pntChrgMsg->coord[i] = compPCVec[i] ;
1659 
1660  }
1661 
1662  DebugM(4,"Rank " << pntChrgMsg->sourceNode << " sending a total of "
1663  << compPCVec.size() << " elements to rank zero." << std::endl);
1664 
1665  CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
1666  QMProxy[0].recvPntChrg(pntChrgMsg);
1667 
1668  delete qmFullMsg;
1669 }
#define PCMODEUPDATEPOS
Definition: ComputeQM.C:100
std::pair< Position, BigReal > PosDistPair
Definition: ComputeQM.C:1290
int sourceNode
Definition: ComputeQM.C:169
Definition: Vector.h:64
ComputeHomePatchList patchList
int * pcIndxs
Definition: ComputeQM.C:111
float Real
Definition: common.h:109
#define DebugM(x, y)
Definition: Debug.h:59
BigReal length(void) const
Definition: Vector.h:169
ComputeQMAtom * coord
Definition: ComputeQM.C:110
Charge charge
Definition: NamdTypes.h:54
void sort(void)
Definition: SortedArray.h:66
#define PCMODEUPDATESEL
Definition: ComputeQM.C:99
Position position
Definition: ComputeQM.C:66
void clear()
Definition: ResizeArray.h:87
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
int sourceNode
Definition: ComputeQM.C:105
Real qmGrpID
Definition: ComputeQM.C:69
int add(const Elem &elem)
Definition: ResizeArray.h:97
int pcSelMode
Definition: ComputeQM.C:109
short vdwType
Definition: NamdTypes.h:55
Mass mass
Definition: NamdTypes.h:108
std::map< Real, std::pair< Position, BigReal > > GrpDistMap
Definition: ComputeQM.C:1289
Elem * find(const Elem &elem)
int numAtoms
Definition: ComputeQM.C:106
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
int size(void) const
Definition: ResizeArray.h:127
#define PCMODECUSTOMSEL
Definition: ComputeQM.C:101
gridSize x
int load(const Elem &elem)
Definition: SortedArray.h:50
int numPCIndxs
Definition: ComputeQM.C:108
double BigReal
Definition: common.h:114
for(int i=0;i< n1;++i)
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().

2673  {
2674 
2676 
2677  bool callReplaceForces = false;
2678 
2679  int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
2680  int numQMGrps = Node::Object()->molecule->get_qmNumGrps();
2681  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
2682  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
2683  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
2684 
2685  int totalNumbAtoms = 0;
2686  for (ap = ap.begin(); ap != ap.end(); ap++) {
2687  totalNumbAtoms += (*ap).p->getNumAtoms();
2688  }
2689 
2690  // This is kept to be deleted in the next step so that the message can be
2691  // used in "replaceForces" routine later on, in case there are forces to
2692  // be replaced. The "replaceForces" function uses but DOES NOT free the pointer,
2693  // so we free the data from the previous iteration and allocate a new one for
2694  // the current iteration (remember that the number of atoms can change in a
2695  // home patch between iterations).
2696  if (oldForces != 0)
2697  delete [] oldForces;
2698  oldForces = new ExtForce[totalNumbAtoms] ;
2699 
2700  for (int i=0; i < totalNumbAtoms; ++i) {
2701  oldForces[i].force = Force(0) ;
2702  }
2703 
2704  DebugM(1,"Force array has been created and zeroed in rank "
2705  << CkMyPe() << std::endl);
2706 
2707  DebugM(1,"Preparing " << fmsg->numForces << " forces in rank "
2708  << CkMyPe() << std::endl);
2709 
2710  QMForce *results_ptr = fmsg->force;
2711  // Iterates over the received forces and place them on the full array.
2712  for (int i=0; i < fmsg->numForces; ++i, ++results_ptr) {
2713  // For now we may have more than one item in the message acting on the same
2714  // atom, such as an MM atom which is a point charge for two or more QM regions.
2715 
2716  oldForces[results_ptr->homeIndx].force += results_ptr->force;
2717  oldForces[results_ptr->homeIndx].replace = results_ptr->replace;
2718 
2719  if (results_ptr->replace == 1)
2720  callReplaceForces = true;
2721 
2722  // If the atom is in a QM group, update its charge to the local (this homePatch)
2723  // copy of the qmAtmChrg array.
2724  if (qmAtomGroup[results_ptr->id] > 0 && (fmsg->PMEOn || (numQMGrps > 1) ) ) {
2725 
2726  // Loops over all QM atoms (in all QM groups) comparing their global indices
2727  for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
2728 
2729  if (qmAtmIndx[qmIter] == results_ptr->id) {
2730  qmAtmChrg[qmIter] = results_ptr->charge;
2731  break;
2732  }
2733 
2734  }
2735 
2736  }
2737 
2738  }
2739 
2740  DebugM(1,"Placing forces on force boxes in rank "
2741  << CkMyPe() << std::endl);
2742 
2743  // Places the received forces on the force array for each patch.
2744  int homeIndxIter = 0;
2745  for (ap = ap.begin(); ap != ap.end(); ap++) {
2746  Results *r = (*ap).forceBox->open();
2747  Force *f = r->f[Results::normal];
2748  const CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
2749  int localNumAtoms = (*ap).p->getNumAtoms();
2750 
2751  for(int i=0; i<localNumAtoms; ++i) {
2752 
2753  f[i] += oldForces[homeIndxIter].force;
2754 
2755  ++homeIndxIter;
2756  }
2757 
2758  if ( callReplaceForces )
2759  (*ap).p->replaceForces(oldForces);
2760 
2761  (*ap).forceBox->close(&r);
2762 
2763  }
2764 
2765  DebugM(1,"Forces placed on force boxes in rank "
2766  << CkMyPe() << std::endl);
2767 
2768  if (fmsg->PMEOn) {
2769 
2770  DebugM(1,"PME ON! Accessing PMEmgr in rank " << CkMyPe() << std::endl);
2771 
2772  ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch(
2773  CkpvAccess(BOCclass_group).computePmeMgr) ;
2774 
2775  DebugM(4, "Initiating ComputePme::doQMWork on rank " << CkMyPe() << " over "
2776  << getComputes(mgrP).size() << " pmeComputes." << std::endl) ;
2777 
2778  for ( int i=0; i< getComputes(mgrP).size(); ++i ) {
2779  // WorkDistrib::messageEnqueueWork(pmeComputes[i]);
2780  getComputes(mgrP)[i]->doQMWork();
2781  }
2782  }
2783 
2784  DebugM(1,"Submitting reduction in rank " << CkMyPe() << std::endl);
2785 
2786  reduction->item(REDUCTION_MISC_ENERGY) += fmsg->energy;
2787  reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += fmsg->virial[0][0];
2788  reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += fmsg->virial[0][1];
2789  reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += fmsg->virial[0][2];
2790  reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += fmsg->virial[1][0];
2791  reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += fmsg->virial[1][1];
2792  reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += fmsg->virial[1][2];
2793  reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += fmsg->virial[2][0];
2794  reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += fmsg->virial[2][1];
2795  reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += fmsg->virial[2][2];
2796  reduction->submit();
2797 
2798  delete fmsg ;
2799 }
static Node * Object()
Definition: Node.h:86
int replace
Definition: ComputeQM.h:102
int get_qmNumGrps() const
Definition: Molecule.h:827
BigReal energy
Definition: ComputeQM.C:187
const int * get_qmAtmIndx()
Definition: Molecule.h:823
Vector Force
Definition: NamdTypes.h:26
Definition: Vector.h:64
ComputeHomePatchList patchList
int get_numQMAtoms()
Definition: Molecule.h:825
int numForces
Definition: ComputeQM.C:189
float Real
Definition: common.h:109
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
int replace
Definition: NamdTypes.h:201
const Real * get_qmAtomGroup() const
Definition: Molecule.h:819
int id
Definition: ComputeQM.h:106
ResizeArray< ComputePme * > & getComputes(ComputePmeMgr *mgr)
Definition: ComputePme.C:587
Force * f[maxNumForces]
Definition: PatchTypes.h:67
Real * get_qmAtmChrg()
Definition: Molecule.h:822
float charge
Definition: ComputeQM.h:105
QMForce * force
Definition: ComputeQM.C:191
Force force
Definition: NamdTypes.h:202
Force force
Definition: ComputeQM.h:103
BigReal virial[3][3]
Definition: ComputeQM.C:188
int homeIndx
Definition: ComputeQM.h:104
void submit(void)
Definition: ReductionMgr.h:323
int size(void) const
Definition: ResizeArray.h:127
Molecule * molecule
Definition: Node.h:176
bool PMEOn
Definition: ComputeQM.C:186

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