#include <HomePatch.h>

Inheritance diagram for HomePatch:

Classes
struct	checkpoint_t

struct	RattleList

Public Member Functions
	~HomePatch ()

void	registerProxy (RegisterProxyMsg *)

void	unregisterProxy (UnregisterProxyMsg *)

void	receiveResults (ProxyResultVarsizeMsg *msg)

void	receiveResults (ProxyResultMsg *msg)

void	receiveResult (ProxyGBISP1ResultMsg *msg)

void	receiveResult (ProxyGBISP2ResultMsg *msg)

void	receiveResults (ProxyCombinedResultRawMsg *msg)

void	depositMigration (MigrateAtomsMsg *)

void	useSequencer (Sequencer *sequencerPtr)

void	runSequencer (void)

void	positionsReady (int doMigration=0)

void	positionsReady_SOA (int doMigration=0)

void	positionsReady_GPU (int doMigration=0, int startup=0)

void	updateAtomCount (const int n, const int reallocate)

void	updateAtomBuffers ()

void	saveForce (const int ftag=Results::normal)

void	addForceToMomentum (FullAtom __restrict atom_arr, const Force __restrict force_arr, const BigReal dt, int num_atoms) __attribute__((__noinline__))

void	addForceToMomentum3 (FullAtom __restrict atom_arr, const Force __restrict force_arr1, const Force __restrict force_arr2, const Force __restrict force_arr3, const BigReal dt1, const BigReal dt2, const BigReal dt3, int num_atoms) __attribute__((__noinline__))

void	addVelocityToPosition (FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms) __attribute__((__noinline__))

int	hardWallDrude (const BigReal, Tensor virial, SubmitReduction )

void	addRattleForce (const BigReal invdt, Tensor &wc)

void	buildRattleList ()

int	rattle1old (const BigReal, Tensor virial, SubmitReduction )

int	rattle1 (const BigReal, Tensor virial, SubmitReduction )

void	rattle2 (const BigReal, Tensor *virial)

void	minimize_rattle2 (const BigReal, Tensor *virial, bool forces=false)

void	buildRattleList_SOA ()

int	rattle1_SOA (const BigReal, Tensor virial, SubmitReduction )

void	mollyAverage ()

void	mollyMollify (Tensor *virial)

void	loweAndersenVelocities ()

void	loweAndersenFinish ()

void	setGBISIntrinsicRadii ()

void	gbisComputeAfterP1 ()

void	gbisComputeAfterP2 ()

void	gbisP2Ready ()

void	gbisP3Ready ()

void	setLcpoType ()

void	checkpoint (void)

void	revert (void)

void	exchangeCheckpoint (int scriptTask, int &bpc)

void	recvCheckpointReq (int task, const char *key, int replica, int pe)

void	recvCheckpointLoad (CheckpointAtomsMsg *msg)

void	recvCheckpointStore (CheckpointAtomsMsg *msg)

void	recvCheckpointAck ()

void	exchangeAtoms (int scriptTask)

void	recvExchangeReq (int req)

void	recvExchangeMsg (ExchangeAtomsMsg *msg)

void	replaceForces (ExtForce *f)

void	qmSwapAtoms ()

void	submitLoadStats (int timestep)

FullAtomList &	getAtomList ()

ScaledPosition	getMin ()

ScaledPosition	getMax ()

void	buildSpanningTree (void)

void	sendNodeAwareSpanningTree ()

void	recvNodeAwareSpanningTree (ProxyNodeAwareSpanningTreeMsg *msg)

void	sendSpanningTree ()

void	recvSpanningTree (int *t, int n)

void	sendProxies ()

Public Member Functions inherited from Patch
	Patch (PatchID pd)

int	hasNewAtoms ()

virtual	~Patch ()

Box< Patch, CompAtom > *	registerPositionPickup (Compute *cid)

void	unregisterPositionPickup (Compute cid, Box< Patch, CompAtom > *const box)

Box< Patch, CompAtom > *	registerAvgPositionPickup (Compute *cid)

void	unregisterAvgPositionPickup (Compute cid, Box< Patch, CompAtom > *const box)

Box< Patch, CompAtom > *	registerVelocityPickup (Compute *cid)

void	unregisterVelocityPickup (Compute cid, Box< Patch, CompAtom > *const box)

Box< Patch, Real > *	registerIntRadPickup (Compute *cid)

void	unregisterIntRadPickup (Compute cid, Box< Patch, Real > *const box)

Box< Patch, GBReal > *	registerPsiSumDeposit (Compute *cid)

void	unregisterPsiSumDeposit (Compute cid, Box< Patch, GBReal > *const box)

Box< Patch, Real > *	registerBornRadPickup (Compute *cid)

void	unregisterBornRadPickup (Compute cid, Box< Patch, Real > *const box)

Box< Patch, GBReal > *	registerDEdaSumDeposit (Compute *cid)

void	unregisterDEdaSumDeposit (Compute cid, Box< Patch, GBReal > *const box)

Box< Patch, Real > *	registerDHdrPrefixPickup (Compute *cid)

void	unregisterDHdrPrefixPickup (Compute cid, Box< Patch, Real > *const box)

Box< Patch, int > *	registerLcpoTypePickup (Compute *cid)

void	unregisterLcpoTypePickup (Compute cid, Box< Patch, int > *const box)

Box< Patch, Results > *	registerForceDeposit (Compute *cid)

void	unregisterForceDeposit (Compute cid, Box< Patch, Results > *const box)

void	positionsReady (int n=0, int startup=1)

void	positionBoxClosed (void)

void	forceBoxClosed (void)

void	avgPositionBoxClosed (void)

void	velocityBoxClosed (void)

void	intRadBoxClosed (void)

void	psiSumBoxClosed (void)

void	bornRadBoxClosed (void)

void	dEdaSumBoxClosed (void)

void	dHdrPrefixBoxClosed (void)

void	gbisP2Ready ()

void	gbisP3Ready ()

void	lcpoTypeBoxClosed (void)

int	getNumAtoms () const

int	getNumFixedAtoms () const

void	setNumFixedAtoms (int numFixed)

PatchID	getPatchID () const

int	getNumComputes () const

CompAtomExt *	getCompAtomExtInfo ()

CompAtom *	getCompAtomInfo ()

CudaAtom *	getCudaAtomList ()

Public Attributes
int	marginViolations

bool	gridForceIdxChecked

std::vector< int >	settleList

std::vector< RattleList >	rattleList

std::vector< RattleParam >	rattleParam

std::vector< int >	noconstList

bool	rattleListValid

bool	rattleListValid_SOA

std::vector< Vector >	velNew

std::vector< Vector >	posNew

int	checkpoint_task

std::map< std::string, checkpoint_t * >	checkpoints

int	exchange_dst

int	exchange_src

int	exchange_req

ExchangeAtomsMsg *	exchange_msg

LDObjHandle	ldObjHandle

Public Attributes inherited from Patch
Lattice &	lattice

Flags	flags

Protected Member Functions
virtual void	boxClosed (int)

void	doPairlistCheck ()

void	doGroupSizeCheck ()

void	doGroupSizeCheck_SOA ()

void	doMarginCheck ()

void	doMarginCheck_SOA ()

void	doAtomMigration ()

Protected Attributes
int	inMigration

int	numMlBuf

MigrateAtomsMsg *	msgbuf [PatchMap::MaxOneAway]

Protected Attributes inherited from Patch
const PatchID	patchID

int	numAtoms

int	numFixedAtoms

CompAtomList	p

CompAtomList	p_avg

CompAtomList	v

AtomMapper *	atomMapper

RealList	intRad

GBRealList	psiSum

GBRealList	psiFin

RealList	bornRad

RealList	dHdrPrefix

GBRealList	dEdaSum

IntList	lcpoType

CompAtomExtList	pExt

CompAtom *	avgPositionPtrBegin

CompAtom *	avgPositionPtrEnd

CompAtom *	velocityPtrBegin

CompAtom *	velocityPtrEnd

CudaAtom *	cudaAtomPtr

ForceList	f [Results::maxNumForces]

Results	results

int	computesSortedByPriority

int	firstHoldableCompute

OwnerBox< Patch, CompAtom >	positionBox

ComputePtrList	positionComputeList

OwnerBox< Patch, CompAtom >	avgPositionBox

ComputePtrList	avgPositionComputeList

OwnerBox< Patch, CompAtom >	velocityBox

ComputePtrList	velocityComputeList

OwnerBox< Patch, Real >	intRadBox

ComputePtrList	intRadComputeList

OwnerBox< Patch, GBReal >	psiSumBox

ComputePtrList	psiSumComputeList

OwnerBox< Patch, Real >	bornRadBox

ComputePtrList	bornRadComputeList

OwnerBox< Patch, GBReal >	dEdaSumBox

ComputePtrList	dEdaSumComputeList

OwnerBox< Patch, Real >	dHdrPrefixBox

ComputePtrList	dHdrPrefixComputeList

OwnerBox< Patch, int >	lcpoTypeBox

ComputePtrList	lcpoTypeComputeList

OwnerBox< Patch, Results >	forceBox

ComputePtrList	forceComputeList

int	boxesOpen

int	_hasNewAtoms

int *	child

int	nChild

Friends
class	PatchMgr

class	Sequencer

class	SequencerCUDA

class	ComputeGlobal

Detailed Description

Definition at line 329 of file HomePatch.h.

Constructor & Destructor Documentation

◆ ~HomePatch()

HomePatch::~HomePatch ( )

Definition at line 357 of file HomePatch.C.

References Patch::atomMapper, ResizeArray< Elem >::begin(), ResizeArray< Elem >::end(), and AtomMapper::unregisterIDsFullAtom().

 {
     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
 #ifdef NODEAWARE_PROXY_SPANNINGTREE
     ptnTree.resize(0);
     #ifdef USE_NODEPATCHMGR
     delete [] nodeChildren;    
     #endif
 #endif
   delete [] child;
 }

Member Function Documentation

◆ addForceToMomentum()

void HomePatch::addForceToMomentum	(	FullAtom *__restrict	atom_arr,
		const Force *__restrict	force_arr,
		const BigReal	dt,
		int	num_atoms
	)

Definition at line 3319 of file HomePatch.C.

References Node::Object(), Node::simParameters, and simParams.

Referenced by Sequencer::addForceToMomentum().

       {
   SimParameters *simParams = Node::Object()->simParameters;
 
   if ( simParams->fixedAtomsOn ) {
     for ( int i = 0; i < num_atoms; ++i ) {
       if ( atom_arr[i].atomFixed ) {
         atom_arr[i].velocity = 0;
       } else {
         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
       }
     }
   } else {
     for ( int i = 0; i < num_atoms; ++i ) {
       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
     }
   }
 }

◆ addForceToMomentum3()

void HomePatch::addForceToMomentum3	(	FullAtom *__restrict	atom_arr,
		const Force *__restrict	force_arr1,
		const Force *__restrict	force_arr2,
		const Force *__restrict	force_arr3,
		const BigReal	dt1,
		const BigReal	dt2,
		const BigReal	dt3,
		int	num_atoms
	)

Definition at line 3348 of file HomePatch.C.

References Node::Object(), Node::simParameters, and simParams.

Referenced by Sequencer::addForceToMomentum3().

       {
   SimParameters *simParams = Node::Object()->simParameters;
 
   if ( simParams->fixedAtomsOn ) {
     for ( int i = 0; i < num_atoms; ++i ) {
       if ( atom_arr[i].atomFixed ) {
         atom_arr[i].velocity = 0;
       } else {
         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
       }
     }
   } else {
     for ( int i = 0; i < num_atoms; ++i ) {
       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
     }
   }
 }

◆ addRattleForce()

void HomePatch::addRattleForce	(	const BigReal	invdt,
		Tensor &	wc
	)

Definition at line 3778 of file HomePatch.C.

References Patch::f, Results::normal, Patch::numAtoms, outer(), and velNew.

Referenced by rattle1().

                                                               {
   for (int ig = 0; ig < numAtoms; ++ig ) {
     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
     Tensor vir = outer(df, atom[ig].position);
     wc += vir;
     f[Results::normal][ig] += df;
     atom[ig].velocity = velNew[ig];
   }
 }

◆ addVelocityToPosition()

void HomePatch::addVelocityToPosition	(	FullAtom *__restrict	atom_arr,
		const BigReal	dt,
		int	num_atoms
	)

Definition at line 3387 of file HomePatch.C.

References Node::Object(), Node::simParameters, and simParams.

Referenced by Sequencer::addVelocityToPosition().

       {
   SimParameters *simParams = Node::Object()->simParameters;
   if ( simParams->fixedAtomsOn ) {
     for ( int i = 0; i < num_atoms; ++i ) {
       if ( ! atom_arr[i].atomFixed ) {
         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
       }
     }
   } else {
     for ( int i = 0; i < num_atoms; ++i ) {
       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
     }
   }
 }

◆ boxClosed()

void HomePatch::boxClosed ( int box )

protectedvirtual

Implements Patch.

Definition at line 370 of file HomePatch.C.

References Sequencer::awaken(), Patch::boxesOpen, DebugM, Patch::dEdaSumBox, Flags::doGBIS, Flags::doNonbonded, Patch::f, Patch::flags, ExtForce::force, OwnerBox< Owner, Data >::isOpen(), Results::maxNumForces, Results::normal, Patch::numAtoms, Node::Object(), Patch::patchID, Patch::psiSumBox, Node::simParameters, simParams, and ResizeArray< Elem >::size().

                                  {
   // begin gbis
   if (box == 5) {// end of phase 1
     phase1BoxClosedCalled = true;
     if (!psiSumBox.isOpen() && numGBISP1Arrived == proxy.size()) {
       if (flags.doGBIS && flags.doNonbonded) {
         // fprintf(stderr, "Calling awaken() on patch %d: 1\n", this->patchID);
         sequencer->awaken();
       }
     } else {
       //need to wait until proxies arrive before awakening
     }
   } else if (box == 6) {// intRad
     //do nothing
   } else if (box == 7) {// bornRad
     //do nothing
   } else if (box == 8) {// end of phase 2
     phase2BoxClosedCalled = true;
     //if no proxies, AfterP1 can't be called from receive
     //so it will be called from here
     if (!dEdaSumBox.isOpen() && numGBISP2Arrived == proxy.size()) {
       if (flags.doGBIS && flags.doNonbonded) {
         // fprintf(stderr, "Calling awaken() on patch %d: 2\n", this->patchID);
         sequencer->awaken();
       }
     } else {
       //need to wait until proxies arrive before awakening
     }
   } else if (box == 9) {
     //do nothing
   } else if (box == 10) {
     //lcpoType Box closed: do nothing
   } else {
     //do nothing
   }
   // end gbis
 
   if ( ! --boxesOpen ) {
     if ( replacementForces ) {
       for ( int i = 0; i < numAtoms; ++i ) {
         if ( replacementForces[i].replace ) {
           for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
           f[Results::normal][i] = replacementForces[i].force;
         }
       }
       replacementForces = 0;
     }
     DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
         << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
     // only awaken suspended threads.  Then say it is suspended.
 
     phase3BoxClosedCalled = true;
     if (flags.doGBIS) {
       if (flags.doNonbonded) {
         sequencer->awaken();
       } else {
         if (numGBISP1Arrived == proxy.size() &&
           numGBISP2Arrived == proxy.size() &&
           numGBISP3Arrived == proxy.size()) {
           sequencer->awaken();//all boxes closed and all proxies arrived
         }
       }
     } else {//non-gbis awaken
       SimParameters *simParams = Node::Object()->simParameters;
       if(!simParams->CUDASOAintegrate) {
         sequencer->awaken();
       }
     }
   } else {
     DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
   }
 }

◆ buildRattleList()

void HomePatch::buildRattleList ( )

Definition at line 3616 of file HomePatch.C.

References RattleParam::dsq, Patch::flags, getWaterModelGroupSize(), RattleParam::ia, RattleParam::ib, HomePatch::RattleList::icnt, HomePatch::RattleList::ig, NAMD_die(), noconstList, Patch::numAtoms, Node::Object(), Patch::patchID, posNew, rattleList, rattleParam, RattleParam::rma, RattleParam::rmb, settle1init(), settleList, Node::simParameters, simParams, Flags::step, SWM4, and velNew.

Referenced by rattle1().

                                 {
 #ifdef DEBUG_MINIMIZE
   if (patchID == 0) {
   printf("Step %d, patch %d: buildRattleList()\n",
       flags.step, (int)patchID);
   }
 #endif
   SimParameters *simParams = Node::Object()->simParameters;
   const int fixedAtomsOn = simParams->fixedAtomsOn;
   const int useSettle = simParams->useSettle;
 
   // Re-size to containt numAtoms elements
   velNew.resize(numAtoms);
   posNew.resize(numAtoms);
 
   // Size of a hydrogen group for water
   const WaterModel watmodel = simParams->watmodel;
   const int wathgsize = getWaterModelGroupSize(watmodel);
 
   // Initialize the settle algorithm with water parameters
   // settle1() assumes all waters are identical,
   // and will generate bad results if they are not.
   // XXX this will move to Molecule::build_atom_status when that 
   // version is debugged
   if ( ! settle_initialized ) {
     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
       // find a water
       if (atom[ig].rigidBondLength > 0) {
         int oatm;
         if (watmodel == WaterModel::SWM4) {
           oatm = ig+3;  // skip over Drude and Lonepair
           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
         }
         else {
           oatm = ig+1;
           // Avoid using the Om site to set this by mistake
           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
             oatm += 1;
           }
         }
 
         // initialize settle water parameters
         settle1init(atom[ig].mass, atom[oatm].mass, 
                     atom[ig].rigidBondLength,
                     atom[oatm].rigidBondLength,
                     settle_mO, settle_mH,
                     settle_mOrmT, settle_mHrmT, settle_ra,
                     settle_rb, settle_rc, settle_rra);
         settle_initialized = 1;
         break; // done with init
       }
     }
   }
   
   Vector ref[10];
   BigReal rmass[10];
   BigReal dsq[10];
   int fixed[10];
   int ial[10];
   int ibl[10];
 
   int numSettle = 0;
   int numRattle = 0;
   int posRattleParam = 0;
 
   settleList.clear();
   rattleList.clear();
   noconstList.clear();
   rattleParam.clear();
 
   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
     int hgs = atom[ig].hydrogenGroupSize;
     if ( hgs == 1 ) {
       // only one atom in group
       noconstList.push_back(ig);
       continue;
     }
     int anyfixed = 0;
     for (int i = 0; i < hgs; ++i ) {
       ref[i] = atom[ig+i].position;
       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
       if ( fixed[i] ) {
         anyfixed = 1;
         rmass[i] = 0.;
       }
     }
     int icnt = 0;
     BigReal tmp = atom[ig].rigidBondLength;
     if (tmp > 0.0) {  // for water
       if (hgs != wathgsize) {
         char errmsg[256];
         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
                          "but the specified water model requires %d atoms.\n",
                          atom[ig].id+1, hgs, wathgsize);
         NAMD_die(errmsg);
       }
       // Use SETTLE for water unless some of the water atoms are fixed,
       if (useSettle && !anyfixed) {
         // Store to Settle -list
         settleList.push_back(ig);
         continue;
       }
       if ( !(fixed[1] && fixed[2]) ) {
         dsq[icnt] = tmp * tmp;
         ial[icnt] = 1;
         ibl[icnt] = 2;
         ++icnt;
       }
     } // if (tmp > 0.0)
     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
         if ( !(fixed[0] && fixed[i]) ) {
           dsq[icnt] = tmp * tmp;
           ial[icnt] = 0;
           ibl[icnt] = i;
           ++icnt;
         }
       }
     }
     if ( icnt == 0 ) {
       // no constraints
       noconstList.push_back(ig);
       continue;  
     }
     // Store to Rattle -list
     RattleList rattleListElem;
     rattleListElem.ig  = ig;
     rattleListElem.icnt = icnt;
     rattleList.push_back(rattleListElem);
     for (int i = 0; i < icnt; ++i ) {
       int a = ial[i];
       int b = ibl[i];
       RattleParam rattleParamElem;
       rattleParamElem.ia = a;
       rattleParamElem.ib = b;
       rattleParamElem.dsq = dsq[i];
       rattleParamElem.rma = rmass[a];
       rattleParamElem.rmb = rmass[b];
       rattleParam.push_back(rattleParamElem);
     }
     //adding dummy atom in the hydrogen group
     for (int i = icnt; i < 4; ++i ) {
       RattleParam rattleParamElem;
       rattleParamElem.ia = 0;
       rattleParamElem.ib = 0;
       rattleParamElem.dsq = 0;
       rattleParamElem.rma = 0;
       rattleParamElem.rmb = 0;
       rattleParam.push_back(rattleParamElem);
     }
 #if 0
     for(int i = 0; i < 4; ++i) {
       std::cout << rattleParam[i].ia << " " << rattleParam[i].ib << std::endl;
     }
     std::cout << std::endl;
 #endif
   }
 
 }

◆ buildRattleList_SOA()

void HomePatch::buildRattleList_SOA ( )

Definition at line 4520 of file HomePatch.C.

References RattleParam::dsq, Patch::flags, getWaterModelGroupSize(), PatchDataSOA::hydrogenGroupSize, RattleParam::ia, RattleParam::ib, HomePatch::RattleList::icnt, HomePatch::RattleList::ig, PatchDataSOA::mass, NAMD_bug(), noconstList, Patch::numAtoms, Node::Object(), Patch::patchID, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, rattleList, rattleParam, PatchDataSOA::recipMass, PatchDataSOA::rigidBondLength, RattleParam::rma, RattleParam::rmb, settle1init(), Node::simParameters, simParams, Flags::step, and SWM4.

Referenced by depositMigration(), rattle1_SOA(), recvCheckpointLoad(), recvExchangeMsg(), and revert().

                                     {
 #ifdef DEBUG_MINIMIZE
   if (patchID == 0) {
   printf("Step %d, patch %d: buildRattleList_SOA()\n",
       flags.step, (int)patchID);
   }
 #endif
 
   // Called when rattleListValid_SOA is false.
   // List will stay valid until atom migration or some other event, 
   // such as exchanging replicas, SCRIPT_REVERT from Tcl, reinit atoms.
 
   const double * __restrict pos_x = patchDataSOA.pos_x;
   const double * __restrict pos_y = patchDataSOA.pos_y;
   const double * __restrict pos_z = patchDataSOA.pos_z;
   const float  * __restrict mass = patchDataSOA.mass;
   const double * __restrict recipMass = patchDataSOA.recipMass;
   const float  * __restrict rigidBondLength = patchDataSOA.rigidBondLength;
   const int    * __restrict hydrogenGroupSize = patchDataSOA.hydrogenGroupSize;
 
   SimParameters *simParams = Node::Object()->simParameters;
   const int fixedAtomsOn = simParams->fixedAtomsOn;
   const int useSettle = simParams->useSettle;
 
   // Size of a hydrogen group for water
   const WaterModel watmodel = simParams->watmodel;
   const int wathgsize = getWaterModelGroupSize(watmodel);
 
   // Initialize the settle algorithm with water parameters
   // settle1() assumes all waters are identical,
   // and will generate bad results if they are not.
   // XXX this will move to Molecule::build_atom_status when that 
   // version is debugged
   if ( ! settle_initialized ) {
     for (int ig = numSoluteAtoms;
          ig < numAtoms;
          ig += hydrogenGroupSize[ig]) {
       // find a water
       if (rigidBondLength[ig] > 0) {
         int oatm;
         if (watmodel == WaterModel::SWM4) {
           oatm = ig+3;  // skip over Drude and Lonepair
           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
         }
         else {
           oatm = ig+1;
           // Avoid using the Om site to set this by mistake
           if (mass[ig] < 0.5 || mass[ig+1] < 0.5) {
             oatm += 1;
           }
         }
 
         // initialize settle water parameters
         settle1init(mass[ig], mass[oatm], 
                     rigidBondLength[ig],
                     rigidBondLength[oatm],
                     settle_mO, settle_mH,
                     settle_mOrmT, settle_mHrmT, settle_ra,
                     settle_rb, settle_rc, settle_rra);
         settle_initialized = 1;
         break; // done with init
       }
     }
   }
   
   BigReal dsq[10];
   int ial[10];
   int ibl[10];
 
   rattleList.clear();
   noconstList.clear();
   rattleParam.clear();
 
   for (int ig = 0;  ig < numSoluteAtoms;  ig += hydrogenGroupSize[ig]) {
     int hgs = hydrogenGroupSize[ig];
     if (hgs == 0) {
       NAMD_bug("Hydrogen group size 0. Exiting to avoid infinite loop\n.");
     } else if ( hgs == 1 ) {
       // only one atom in group
       noconstList.push_back(ig);
       continue;
     }
     int icnt = 0;
     // XXX convert rigid bond length to double to square it?
     BigReal tmp = rigidBondLength[ig];
     if (tmp > 0.0) {  // for water
       dsq[icnt] = tmp * tmp;
       ial[icnt] = 1;
       ibl[icnt] = 2;
       ++icnt;
     }
     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
       if ( ( tmp = rigidBondLength[ig+i] ) > 0 ) {
         dsq[icnt] = tmp * tmp;
         ial[icnt] = 0;
         ibl[icnt] = i;
         ++icnt;
       }
     }
     if ( icnt == 0 ) {
       // no constraints
       noconstList.push_back(ig);
       continue;
     }
     // Store to Rattle -list
     RattleList rattleListElem;
     rattleListElem.ig  = ig;
     rattleListElem.icnt = icnt;
     rattleList.push_back(rattleListElem);
     for (int i = 0; i < icnt; ++i ) {
       int a = ial[i];
       int b = ibl[i];
       RattleParam rattleParamElem;
       rattleParamElem.ia = a;
       rattleParamElem.ib = b;
       rattleParamElem.dsq = dsq[i];
       rattleParamElem.rma = recipMass[ig+a];
       rattleParamElem.rmb = recipMass[ig+b];
       rattleParam.push_back(rattleParamElem);
     }
     //adding dummy atom in the hydrogen group
     for (int i = icnt; i < 4; ++i )
     {
         RattleParam rattleParamElem;
         rattleParamElem.ia = 0;
         rattleParamElem.ib = 0;
         rattleParamElem.dsq = 0;
         rattleParamElem.rma = 0;
         rattleParamElem.rmb = 0;
         rattleParam.push_back(rattleParamElem);
     }
 
   }
 }

◆ buildSpanningTree()

void HomePatch::buildSpanningTree ( void )

Definition at line 715 of file HomePatch.C.

References ResizeArray< Elem >::begin(), compDistance(), ResizeArray< Elem >::end(), ResizeArray< Elem >::find(), NAMD_bug(), PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), PatchMap::Object(), Patch::patchID, proxySpanDim, ResizeArray< Elem >::resize(), sendSpanningTree(), ResizeArray< Elem >::setall(), ResizeArray< Elem >::size(), and ResizeArray< Elem >::swap().

Referenced by ProxyMgr::buildProxySpanningTree().

 {
   nChild = 0;
   int psize = proxy.size();
   if (psize == 0) return;
   NodeIDList oldtree;  oldtree.swap(tree);
   int oldsize = oldtree.size();
   tree.resize(psize + 1);
   tree.setall(-1);
   tree[0] = CkMyPe();
   int s=1, e=psize+1;
   NodeIDList::iterator pli;
   int patchNodesLast =
     ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
   int nNonPatch = 0;
 #if 1
     // try to put it to the same old tree
   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
   {
     int oldindex = oldtree.find(*pli);
     if (oldindex != -1 && oldindex < psize) {
       tree[oldindex] = *pli;
     }
   }
   s=1; e=psize;
   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
   {
     if (tree.find(*pli) != -1) continue;    // already assigned
     if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
       while (tree[e] != -1) { e--; if (e==-1) e = psize; }
       tree[e] = *pli;
     } else {
       while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
       tree[s] = *pli;
       nNonPatch++;
     }
   }
 #if 1
   if (oldsize==0 && nNonPatch) {
     // first time, sort by distance
     qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
   }
 #endif
 
   //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
 
 #if USE_TOPOMAP && USE_SPANNING_TREE
   
   //Right now only works for spanning trees with two levels
   int *treecopy = new int [psize];
   int subroots[proxySpanDim];
   int subsizes[proxySpanDim];
   int subtrees[proxySpanDim][proxySpanDim];
   int idxes[proxySpanDim];
   int i = 0;
 
   for(i=0;i<proxySpanDim;i++){
     subsizes[i] = 0;
     idxes[i] = i;
   }
   
   for(i=0;i<psize;i++){
     treecopy[i] = tree[i+1];
   }
   
   TopoManager tmgr;
   tmgr.sortRanksByHops(treecopy,nNonPatch);
   tmgr.sortRanksByHops(treecopy+nNonPatch,
                                                 psize-nNonPatch);  
   
   /* build tree and subtrees */
   nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
   delete [] treecopy;
   
   for(int i=1;i<psize+1;i++){
     int isSubroot=0;
     for(int j=0;j<nChild;j++)
       if(tree[i]==subroots[j]){
         isSubroot=1;
         break;
       }
     if(isSubroot) continue;
     
     int bAdded = 0;
     tmgr.sortIndexByHops(tree[i], subroots,
                                                   idxes, proxySpanDim);
     for(int j=0;j<proxySpanDim;j++){
       if(subsizes[idxes[j]]<proxySpanDim){
         subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
         bAdded = 1; 
         break;
       }
     }
     if( psize > proxySpanDim && ! bAdded ) {
       NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
     }
   }
 
 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
   
   for (int i=1; i<=proxySpanDim; i++) {
     if (tree.size() <= i) break;
     child[i-1] = tree[i];
     nChild++;
   }
 #endif
 #endif
   
 #if 0
   // for debugging
   CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
   for (int i=0; i<psize+1; i++) {
     CkPrintf("%d ", tree[i]);
   }
   CkPrintf("\n");
 #endif
   // send to children nodes
   sendSpanningTree();
 }

◆ checkpoint()

void HomePatch::checkpoint ( void )

Definition at line 5222 of file HomePatch.C.

References ResizeArray< Elem >::copy(), and Patch::lattice.

Referenced by Sequencer::algorithm().

                                {
   checkpoint_atom.copy(atom);
   checkpoint_lattice = lattice;
 
   // DMK - Atom Separation (water vs. non-water)
   #if NAMD_SeparateWaters != 0
     checkpoint_numWaterAtoms = numWaterAtoms;
   #endif
 }

◆ depositMigration()

void HomePatch::depositMigration ( MigrateAtomsMsg * msg )

Definition at line 5962 of file HomePatch.C.

References ResizeArray< Elem >::add(), Lattice::apply_transform(), Sequencer::awaken(), ResizeArray< Elem >::begin(), buildRattleList_SOA(), DebugM, ResizeArray< Elem >::end(), CompAtom::hydrogenGroupSize, CompAtomExt::id, inMigration, Patch::lattice, FullAtom::migrationGroupSize, MigrateAtomsMsg::migrationList, msgbuf, Lattice::nearest(), Patch::numAtoms, numMlBuf, Node::Object(), Patch::patchID, CompAtom::position, rattleListValid_SOA, Lattice::reverse_transform(), RIGID_NONE, Node::simParameters, simParams, ResizeArray< Elem >::size(), and FullAtom::transform.

Referenced by MigrateAtomsCombinedMsg::distribute(), and doAtomMigration().

 {
 
   if (!inMigration) { // We have to buffer changes due to migration
                       // until our patch is in migration mode
     msgbuf[numMlBuf++] = msg;
     return;
   } 
 
   // DMK - Atom Separation (water vs. non-water)
   #if NAMD_SeparateWaters != 0
 
 
     // Merge the incoming list of atoms with the current list of
     //   atoms.  Note that mergeSeparatedAtomList() will apply any
     //   required transformations to the incoming atoms as it is
     //   separating them.
     mergeAtomList(msg->migrationList);
 
 
   #else
 
     // Merge the incoming list of atoms with the current list of
     // atoms.  Apply transformations to the incoming atoms as they are
     // added to this patch's list.
     {
       MigrationList &migrationList = msg->migrationList;
       MigrationList::iterator mi;
       Transform mother_transform;
       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
         DebugM(1,"Migrating atom " << mi->id << " to patch "
                   << patchID << " with position " << mi->position << "\n"); 
         if ( mi->migrationGroupSize ) {
           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
             Position pos = mi->position;
             int mgs = mi->migrationGroupSize;
             int c = 1;
             for ( int j=mi->hydrogenGroupSize; j<mgs;
                                 j+=(mi+j)->hydrogenGroupSize ) {
               pos += (mi+j)->position;
               ++c;
             }
             pos *= 1./c;
             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
             mother_transform = mi->transform;
             pos = lattice.nearest(pos,center,&mother_transform);
             mi->position = lattice.reverse_transform(mi->position,mi->transform);
             mi->position = lattice.apply_transform(mi->position,mother_transform);
             mi->transform = mother_transform;
           } else {
             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
             mother_transform = mi->transform;
           }
         } else {
           mi->position = lattice.reverse_transform(mi->position,mi->transform);
           mi->position = lattice.apply_transform(mi->position,mother_transform);
           mi->transform = mother_transform;
         }
         atom.add(*mi);
       }
     }
 
 
   #endif // if (NAMD_SeparateWaters != 0)
 
 
   numAtoms = atom.size();
   delete msg;
 
   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
   if (!--patchMigrationCounter) {
     // DH - All atoms are now incorporated from migration.
     // This is where we can separate waters from non-waters and
     // perhaps sort non-waters by hydrogen group size.
     SimParameters *simParams = Node::Object()->simParameters;
     if (simParams->SOAintegrateOn) {
       sort_solvent_atoms();
       copy_atoms_to_SOA();
 #if 0
       if (simParams->rigidBonds != RIGID_NONE) {
         buildRattleList_SOA();
         rattleListValid_SOA = true;
       }
 #endif
     }
     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
     allMigrationIn = true;
     patchMigrationCounter = numNeighbors;
     if (migrationSuspended) {
       DebugM(3,"patch "<<patchID<<" is being awakened\n");
       migrationSuspended = false;
       // fprintf(stderr, "Calling awaken() on patch %d: 8\n", this->patchID);
       sequencer->awaken();
       return;
     }
     else {
        DebugM(3,"patch "<<patchID<<" was not suspended\n");
     }
   }
 }

◆ doAtomMigration()

void HomePatch::doAtomMigration ( )

protected

Definition at line 5798 of file HomePatch.C.

References ResizeArray< Elem >::add(), Patch::atomMapper, ResizeArray< Elem >::begin(), DebugM, ResizeArray< Elem >::del(), depositMigration(), ResizeArray< Elem >::end(), Patch::getPatchID(), CompAtom::hydrogenGroupSize, CompAtomExt::id, inMigration, Patch::lattice, marginViolations, FullAtom::migrationGroupSize, MigrationInfo::mList, msgbuf, NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NamdProfileEventStr, Patch::numAtoms, numMlBuf, PatchMgr::Object(), Patch::patchID, CompAtom::position, ResizeArray< Elem >::resize(), Lattice::scale(), PatchMgr::sendMigrationMsgs(), ResizeArray< Elem >::size(), Sequencer::suspend(), AtomMapper::unregisterIDsFullAtom(), Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady(), and positionsReady_SOA().

 {
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
   char ambuf[32];
   sprintf(ambuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::ATOM_MIGRATIONS], this->getPatchID());
   NAMD_EVENT_START_EX(1, NamdProfileEvent::ATOM_MIGRATIONS, ambuf);
 #endif
 
   // every patch needs to call this once per migration step
   // XXX TODO: check if the cpu version also calls it once per tstep
   int i;
   for (i=0; i<numNeighbors; i++) {
     realInfo[i].mList.resize(0);
   }
 
   // Purge the AtomMap
   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
 
   // Determine atoms that need to migrate
 
   BigReal minx = min.x;
   BigReal miny = min.y;
   BigReal minz = min.z;
   BigReal maxx = max.x;
   BigReal maxy = max.y;
   BigReal maxz = max.z;
 
   int xdev, ydev, zdev;
   int delnum = 0;
 
   FullAtomList::iterator atom_i = atom.begin();
   FullAtomList::iterator atom_e = atom.end();
 
   // DMK - Atom Separation (water vs. non-water)
   #if NAMD_SeparateWaters != 0
     FullAtomList::iterator atom_first = atom_i;
     int numLostWaterAtoms = 0;
   #endif
 
   while ( atom_i != atom_e ) {
     // Even though this code iterates through all atoms successively
     // it moves entire hydrogen/migration groups as follows:
     // Only the parent atom of the hydrogen/migration group has
     // nonzero migrationGroupSize.  Values determined for xdev,ydev,zdev
     // will persist through the remaining group members so that each
     // following atom will again be added to the same mList.
     if ( atom_i->migrationGroupSize ) {
       Position pos = atom_i->position;
       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
         // If there are multiple hydrogen groups in a migration group
         // (e.g. for supporting lone pairs)
         // the following code takes the average position (midpoint)
         // of their parents.
         int mgs = atom_i->migrationGroupSize;
         int c = 1;
         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
                                 j+=(atom_i+j)->hydrogenGroupSize ) {
           pos += (atom_i+j)->position;
           ++c;
         }
         pos *= 1./c;
         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
       }
 
       // Scaling the position below transforms space within patch from
       // what could have been a rotated parallelepiped into
       // orthogonal coordinates, where we can use minmax comparison
       // to detect which of our nearest neighbors this
       // parent atom might have entered.
       ScaledPosition s = lattice.scale(pos);
 
       // check if atom is within bounds
       if (s.x < minx) xdev = 0;
       else if (maxx <= s.x) xdev = 2;
       else xdev = 1;
 
       if (s.y < miny) ydev = 0;
       else if (maxy <= s.y) ydev = 2;
       else ydev = 1;
 
       if (s.z < minz) zdev = 0;
       else if (maxz <= s.z) zdev = 2;
       else zdev = 1;
 
     }
 
     if (mInfo[xdev][ydev][zdev]) { // process atom for migration
                                     // Don't migrate if destination is myself
 
       // See if we have a migration list already
       MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
       DebugM(3,"Migrating atom " << atom_i->id << " from patch "
                 << patchID << " with position " << atom_i->position << "\n");
       mCur.add(*atom_i);
       ++delnum;
 
 
       // DMK - Atom Separation (water vs. non-water)
       #if NAMD_SeparateWaters != 0
         // Check to see if this atom is part of a water molecule.  If
         //   so, numWaterAtoms needs to be adjusted to refect the lost of
         //   this atom.
         // NOTE: The atom separation code assumes that if the oxygen
         //   atom of the hydrogen group making up the water molecule is
         //   migrated to another HomePatch, the hydrogens will also
         //   move!!!
         int atomIndex = atom_i - atom_first;
         if (atomIndex < numWaterAtoms)
           numLostWaterAtoms++;
       #endif
 
 
     } else {
       // By keeping track of delnum total being deleted from FullAtomList
       // the else clause allows us to fill holes as we visit each atom.
 
       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
 
     }
 
     ++atom_i;
   }
 
   // DMK - Atom Separation (water vs. non-water)
   #if NAMD_SeparateWaters != 0
     numWaterAtoms -= numLostWaterAtoms;
   #endif
 
 
   int delpos = numAtoms - delnum;
   DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
   atom.del(delpos,delnum);
 
   numAtoms = atom.size();
   // Calls sendMigrationMsgs to the manager.
   // the manager only
   // wait, does this work??????
   PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
 
   // signal depositMigration() that we are inMigration mode
   inMigration = true;
 
   // Drain the migration message buffer
   for (i=0; i<numMlBuf; i++) {
      DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
      depositMigration(msgbuf[i]);
   }
   numMlBuf = 0;
 
   NAMD_EVENT_STOP(1, NamdProfileEvent::ATOM_MIGRATIONS);
 
   if (!allMigrationIn) {
     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
     migrationSuspended = true;
     sequencer->suspend();
     migrationSuspended = false;
   }
   allMigrationIn = false;
 
   inMigration = false;
   marginViolations = 0;
 }

◆ doGroupSizeCheck()

void HomePatch::doGroupSizeCheck ( )

protected

Definition at line 5590 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Flags::doNonbonded, ResizeArray< Elem >::end(), Patch::flags, CompAtom::hydrogenGroupSize, Flags::maxGroupRadius, NAMD_bug(), CompAtom::nonbondedGroupSize, Node::Object(), CompAtom::position, Node::simParameters, simParams, Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

 {
   if ( ! flags.doNonbonded ) return;
 
   SimParameters *simParams = Node::Object()->simParameters;
   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
   BigReal maxrad2 = 0.;
 
   FullAtomList::iterator p_i = atom.begin();
   FullAtomList::iterator p_e = atom.end();
 
   while ( p_i != p_e ) {
     const int hgs = p_i->hydrogenGroupSize;
     if ( ! hgs ) break;  // avoid infinite loop on bug
     int ngs = hgs;
     if ( ngs > 5 ) ngs = 5;  // XXX why? limit to at most 5 atoms per group
     BigReal x = p_i->position.x;
     BigReal y = p_i->position.y;
     BigReal z = p_i->position.z;
     int i;
     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
       p_i[i].nonbondedGroupSize = 0;
       BigReal dx = p_i[i].position.x - x;
       BigReal dy = p_i[i].position.y - y;
       BigReal dz = p_i[i].position.z - z;
       BigReal r2 = dx * dx + dy * dy + dz * dz;
       if ( r2 > hgcut ) break;
       else if ( r2 > maxrad2 ) maxrad2 = r2;
     }
     p_i->nonbondedGroupSize = i;
     for ( ; i < hgs; ++i ) {
       p_i[i].nonbondedGroupSize = 1;
     }
     p_i += hgs;
   }
 
   if ( p_i != p_e ) {
     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
   }
 
   flags.maxGroupRadius = sqrt(maxrad2);
 
 }

◆ doGroupSizeCheck_SOA()

void HomePatch::doGroupSizeCheck_SOA ( )

protected

Definition at line 5537 of file HomePatch.C.

References Flags::doNonbonded, Patch::flags, PatchDataSOA::hydrogenGroupSize, Flags::maxGroupRadius, NAMD_bug(), PatchDataSOA::nonbondedGroupSize, Patch::numAtoms, Node::Object(), PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, Node::simParameters, and simParams.

Referenced by positionsReady_SOA().

 {
   if ( ! flags.doNonbonded ) return;
 
   SimParameters *simParams = Node::Object()->simParameters;
   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
   BigReal maxrad2 = 0.;
 
   double * __restrict pos_x = patchDataSOA.pos_x;
   double * __restrict pos_y = patchDataSOA.pos_y;
   double * __restrict pos_z = patchDataSOA.pos_z;
   int * __restrict hydrogenGroupSize = patchDataSOA.hydrogenGroupSize;
   int * __restrict nonbondedGroupSize = patchDataSOA.nonbondedGroupSize;
 
   int j=0;
   while (j < numAtoms) {
     const int hgs = hydrogenGroupSize[j];
     if ( ! hgs ) break;  // avoid infinite loop on bug
     int ngs = hgs;
     if ( ngs > 5 ) ngs = 5;  // XXX why? limit to at most 5 atoms per group
     BigReal x = pos_x[j];
     BigReal y = pos_y[j];
     BigReal z = pos_z[j];
     int i;
     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
       nonbondedGroupSize[j+i] = 0;
       BigReal dx = pos_x[j+i] - x;
       BigReal dy = pos_y[j+i] - y;
       BigReal dz = pos_z[j+i] - z;
       BigReal r2 = dx * dx + dy * dy + dz * dz;
       if ( r2 > hgcut ) break;
       else if ( r2 > maxrad2 ) maxrad2 = r2;
     }
     nonbondedGroupSize[j] = i;
     for ( ; i < hgs; ++i ) {
       nonbondedGroupSize[j+i] = 1;
     }
     j += hgs;
   }
 
   if (j != numAtoms) {
     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
   }
 
   flags.maxGroupRadius = sqrt(maxrad2);
 
 }

◆ doMarginCheck()

void HomePatch::doMarginCheck ( )

protected

Definition at line 5721 of file HomePatch.C.

References Lattice::a(), Lattice::a_r(), Lattice::b(), Lattice::b_r(), ResizeArray< Elem >::begin(), Lattice::c(), Lattice::c_r(), ResizeArray< Elem >::end(), Patch::lattice, marginViolations, NAMD_die(), Node::Object(), CompAtom::position, Lattice::scale(), Node::simParameters, simParams, Vector::unit(), Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

 {
   SimParameters *simParams = Node::Object()->simParameters;
 
   BigReal sysdima = lattice.a_r().unit() * lattice.a();
   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
 
   BigReal minSize = simParams->patchDimension - simParams->margin;
 
   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
 
     NAMD_die("Periodic cell has become too small for original patch grid!\n"
       "Possible solutions are to restart from a recent checkpoint,\n"
       "increase margin, or disable useFlexibleCell for liquid simulation.");
   }
 
   BigReal cutoff = simParams->cutoff;
 
   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
 
   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
     NAMD_die("Periodic cell has become too small for original patch grid!\n"
       "There are probably many margin violations already reported.\n"
       "Possible solutions are to restart from a recent checkpoint,\n"
       "increase margin, or disable useFlexibleCell for liquid simulation.");
   }
 
   BigReal minx = min.x - margina;
   BigReal miny = min.y - marginb;
   BigReal minz = min.z - marginc;
   BigReal maxx = max.x + margina;
   BigReal maxy = max.y + marginb;
   BigReal maxz = max.z + marginc;
 
   int xdev, ydev, zdev;
   int problemCount = 0;
 
   FullAtomList::iterator p_i = atom.begin();
   FullAtomList::iterator p_e = atom.end();
   for ( ; p_i != p_e; ++p_i ) {
 
     ScaledPosition s = lattice.scale(p_i->position);
 
     // check if atom is within bounds
     if (s.x < minx) xdev = 0;
     else if (maxx <= s.x) xdev = 2; 
     else xdev = 1;
 
     if (s.y < miny) ydev = 0;
     else if (maxy <= s.y) ydev = 2; 
     else ydev = 1;
 
     if (s.z < minz) zdev = 0;
     else if (maxz <= s.z) zdev = 2; 
     else zdev = 1;
 
     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
             ++problemCount;
     }
 
   }
 
   marginViolations = problemCount;
   // if ( problemCount ) {
   //     iout << iERROR <<
   //       "Found " << problemCount << " margin violations!\n" << endi;
   // } 
 
 }

◆ doMarginCheck_SOA()

void HomePatch::doMarginCheck_SOA ( )

protected

Definition at line 5639 of file HomePatch.C.

References Lattice::a(), Lattice::a_r(), Lattice::b(), Lattice::b_r(), Lattice::c(), Lattice::c_r(), Patch::lattice, marginViolations, NAMD_die(), Patch::numAtoms, Node::Object(), PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, Lattice::scale(), Node::simParameters, simParams, Vector::unit(), Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady_SOA().

 {
   SimParameters *simParams = Node::Object()->simParameters;
 
   BigReal sysdima = lattice.a_r().unit() * lattice.a();
   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
 
   BigReal minSize = simParams->patchDimension - simParams->margin;
 
   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
 
     NAMD_die("Periodic cell has become too small for original patch grid!\n"
       "Possible solutions are to restart from a recent checkpoint,\n"
       "increase margin, or disable useFlexibleCell for liquid simulation.");
   }
 
   BigReal cutoff = simParams->cutoff;
 
   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
 
   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
     NAMD_die("Periodic cell has become too small for original patch grid!\n"
       "There are probably many margin violations already reported.\n"
       "Possible solutions are to restart from a recent checkpoint,\n"
       "increase margin, or disable useFlexibleCell for liquid simulation.");
   }
 
   BigReal minx = min.x - margina;
   BigReal miny = min.y - marginb;
   BigReal minz = min.z - marginc;
   BigReal maxx = max.x + margina;
   BigReal maxy = max.y + marginb;
   BigReal maxz = max.z + marginc;
 
   int xdev, ydev, zdev;
   int problemCount = 0;
 
   double * __restrict pos_x = patchDataSOA.pos_x;
   double * __restrict pos_y = patchDataSOA.pos_y;
   double * __restrict pos_z = patchDataSOA.pos_z;
   for (int i=0;  i < numAtoms;  i++) {
     Vector pos(pos_x[i],pos_y[i],pos_z[i]);
 
     ScaledPosition s = lattice.scale(pos);
 
     // check if atom is within bounds
     if (s.x < minx) xdev = 0;
     else if (maxx <= s.x) xdev = 2; 
     else xdev = 1;
 
     if (s.y < miny) ydev = 0;
     else if (maxy <= s.y) ydev = 2; 
     else ydev = 1;
 
     if (s.z < minz) zdev = 0;
     else if (maxz <= s.z) zdev = 2; 
     else zdev = 1;
 
     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
         ++problemCount;
     }
 
   }
 
   marginViolations = problemCount;
   // if ( problemCount ) {
   //     iout << iERROR <<
   //       "Found " << problemCount << " margin violations!\n" << endi;
   // } 
 
 }

◆ doPairlistCheck()

void HomePatch::doPairlistCheck ( )

protected

Definition at line 5438 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Patch::flags, Patch::getPatchID(), Patch::lattice, Vector::length2(), Flags::maxAtomMovement, NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NamdProfileEventStr, Patch::numAtoms, Node::Object(), Patch::p, Flags::pairlistTolerance, CompAtom::position, ResizeArray< Elem >::resize(), Flags::savePairlists, Node::simParameters, simParams, Lattice::unscale(), Flags::usePairlists, Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady(), and positionsReady_SOA().

 {
 #if 0
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
   char dpcbuf[32];
   sprintf(dpcbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::DO_PAIRLIST_CHECK], this->getPatchID());
   NAMD_EVENT_START_EX(1, NamdProfileEvent::DO_PAIRLIST_CHECK, dpcbuf);
 #endif
 #endif
 
   SimParameters *simParams = Node::Object()->simParameters;
 
   if ( numAtoms == 0 || ! flags.usePairlists ) {
     flags.pairlistTolerance = 0.;
     flags.maxAtomMovement = 99999.;
 #if 0
     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
 #endif
     return;
   }
 
   int i; int n = numAtoms;
   CompAtom *p_i = p.begin();
 
   if ( flags.savePairlists ) {
     flags.pairlistTolerance = doPairlistCheck_newTolerance;
     flags.maxAtomMovement = 0.;
     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
     doPairlistCheck_lattice = lattice;
     doPairlistCheck_positions.resize(numAtoms);
     CompAtom *psave_i = doPairlistCheck_positions.begin();
     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
 #if 0
     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
 #endif
     return;
   }
 
   Lattice &lattice_old = doPairlistCheck_lattice;
   Position center_cur = lattice.unscale(center);
   Position center_old = lattice_old.unscale(center);
   Vector center_delta = center_cur - center_old;
   
   // find max deviation to corner (any neighbor shares a corner)
   BigReal max_cd = 0.;
   for ( i=0; i<2; ++i ) {
     for ( int j=0; j<2; ++j ) {
       for ( int k=0; k<2; ++k ) {
               ScaledPosition corner(    i ? min.x : max.x ,     j ? min.y : max.y ,     k ? min.z : max.z );
               Vector corner_delta =     lattice.unscale(corner) - lattice_old.unscale(corner);
         corner_delta -= center_delta;
               BigReal cd = corner_delta.length2();
         if ( cd > max_cd ) max_cd = cd;
       }
     }
   }
   max_cd = sqrt(max_cd);
 
   // find max deviation of atoms relative to center
   BigReal max_pd = 0.;
   CompAtom *p_old_i = doPairlistCheck_positions.begin();
   for ( i=0; i<n; ++i ) {
     // JM: Calculating position difference and making it patch-centered
     Vector p_delta = p_i[i].position - p_old_i[i].position;
     p_delta -= center_delta;
     BigReal pd = p_delta.length2();
     if ( pd > max_pd ) max_pd = pd;
   }
   max_pd = sqrt(max_pd);
 
   BigReal max_tol = max_pd + max_cd;
 
   flags.maxAtomMovement = max_tol;
 
   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
 
   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
                                 doPairlistCheck_newTolerance ) ) {
     //if(this->getPatchID() == 0) fprintf(stderr, "CPU: Increasing pairList tolerance(%lf %lf)\n", 
     //  max_tol, doPairlistCheck_newTolerance);
     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
   }
 
   if ( max_tol > doPairlistCheck_newTolerance ) {
     //if(this->getPatchID() == 0) fprintf(stderr, "CPU: Decreasing pairList tolerance(%lf %lf)\n", 
     //  max_tol, doPairlistCheck_newTolerance);
     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
   }
   
   //if(this->getPatchID() == 0) fprintf(stderr, "CPU: New patchTolerance: %lf\n", doPairlistCheck_newTolerance);
 
 //  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
 }

◆ exchangeAtoms()

void HomePatch::exchangeAtoms ( int scriptTask )

Definition at line 5370 of file HomePatch.C.

References ExchangeAtomsMsg::atoms, ResizeArray< Elem >::begin(), exchange_dst, exchange_msg, exchange_req, exchange_src, ExchangeAtomsMsg::lattice, Patch::lattice, ExchangeAtomsMsg::numAtoms, Patch::numAtoms, Node::Object(), PatchMgr::Object(), Patch::patchID, ExchangeAtomsMsg::pid, recvExchangeReq(), SCRIPT_ATOMRECV, SCRIPT_ATOMSEND, SCRIPT_ATOMSENDRECV, PatchMgr::sendExchangeReq(), Node::simParameters, and simParams.

Referenced by Sequencer::algorithm().

                                             {
   SimParameters *simParams = Node::Object()->simParameters;
   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
     exchange_dst = (int) simParams->scriptArg1;
     // create and save outgoing message
     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
     exchange_msg->lattice = lattice;
     exchange_msg->pid = patchID;
     exchange_msg->numAtoms = numAtoms;
     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
     if ( exchange_req >= 0 ) {
       recvExchangeReq(exchange_req);
     }
   }
   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
     exchange_src = (int) simParams->scriptArg2;
     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
   }
 }

◆ exchangeCheckpoint()

void HomePatch::exchangeCheckpoint	(	int	scriptTask,
		int &	bpc
	)

Definition at line 5263 of file HomePatch.C.

References checkpoint_task, Node::Object(), PatchMgr::Object(), Patch::patchID, PatchMgr::sendCheckpointReq(), Node::simParameters, and simParams.

Referenced by Sequencer::algorithm().

                                                            {  // initiating replica
   SimParameters *simParams = Node::Object()->simParameters;
   checkpoint_task = scriptTask;
   const int remote = simParams->scriptIntArg1;
   const char *key = simParams->scriptStringArg1;
   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
 }

◆ gbisComputeAfterP1()

void HomePatch::gbisComputeAfterP1 ( )

Definition at line 4915 of file HomePatch.C.

References Patch::bornRad, Patch::flags, gbisP2Ready(), Patch::intRad, Patch::numAtoms, Node::Object(), Patch::pExt, Patch::psiFin, Patch::psiSum, Flags::sequence, Node::simParameters, and simParams.

Referenced by Sequencer::runComputeObjects().

                                    {
 
   SimParameters *simParams = Node::Object()->simParameters;
   BigReal alphaMax = simParams->alpha_max;
   BigReal delta = simParams->gbis_delta;
   BigReal beta = simParams->gbis_beta;
   BigReal gamma = simParams->gbis_gamma;
   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
 
   BigReal rhoi;
   BigReal rhoi0;
   //calculate bornRad from psiSum
   for (int i = 0; i < numAtoms; i++) {
     rhoi0 = intRad[2*i];
     rhoi = rhoi0+coulomb_radius_offset;
     psiFin[i] += psiSum[i];
     psiFin[i] *= rhoi0;
     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
 #ifdef PRINT_COMP
     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
 #endif
   }
 
   gbisP2Ready();
 }

◆ gbisComputeAfterP2()

void HomePatch::gbisComputeAfterP2 ( )

Definition at line 4943 of file HomePatch.C.

References Patch::bornRad, COULOMB, Patch::dEdaSum, Patch::dHdrPrefix, Patch::flags, gbisP3Ready(), Patch::intRad, Patch::numAtoms, Node::Object(), Patch::pExt, Patch::psiFin, Flags::sequence, Node::simParameters, and simParams.

Referenced by Sequencer::runComputeObjects().

                                    {
 
   SimParameters *simParams = Node::Object()->simParameters;
   BigReal delta = simParams->gbis_delta;
   BigReal beta = simParams->gbis_beta;
   BigReal gamma = simParams->gbis_gamma;
   BigReal epsilon_s = simParams->solvent_dielectric;
   BigReal epsilon_p = simParams->dielectric;
   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
   BigReal epsilon_p_i = 1/simParams->dielectric;
   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
   BigReal kappa = simParams->kappa;
   BigReal fij, expkappa, Dij, dEdai, dedasum;
   BigReal rhoi, rhoi0, psii, nbetapsi;
   BigReal gammapsi2, tanhi, daidr;
   for (int i = 0; i < numAtoms; i++) {
     //add diagonal dEda term
     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
     fij = bornRad[i];//inf
     expkappa = exp(-kappa*fij);//0
     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
     //calculate dHij prefix
     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
     dHdrPrefix[i] += dEdai;
     dedasum = dHdrPrefix[i];
 
     rhoi0 = intRad[2*i];
     rhoi = rhoi0+coulomb_radius_offset;
     psii = psiFin[i];
     nbetapsi = -beta*psii;
     gammapsi2 = gamma*psii*psii;
     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
 #ifdef PRINT_COMP
     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
 #endif
   }
   gbisP3Ready();
 }

◆ gbisP2Ready()

void HomePatch::gbisP2Ready ( )

Definition at line 4987 of file HomePatch.C.

References ResizeArray< Elem >::begin(), ProxyGBISP2DataMsg::bornRad, Patch::bornRad, ProxyGBISP2DataMsg::destPe, Patch::flags, GB2_PROXY_DATA_PRIORITY, Patch::gbisP2Ready(), Patch::numAtoms, ProxyGBISP2DataMsg::origPe, ProxyGBISP2DataMsg::patch, PATCH_PRIORITY, Patch::patchID, PRIORITY_SIZE, Flags::sequence, SET_PRIORITY, and ResizeArray< Elem >::size().

Referenced by gbisComputeAfterP1().

                             {
   if (proxy.size() > 0) {
     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
     for (int i = 0; i < proxy.size(); i++) {
       int node = proxy[i];
       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
       msg->patch = patchID;
       msg->origPe = CkMyPe();
       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
       msg->destPe = node;
       int seq = flags.sequence;
       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
       SET_PRIORITY(msg,seq,priority);
       cp[node].recvData(msg);
     }
   }
   Patch::gbisP2Ready();
 }

◆ gbisP3Ready()

void HomePatch::gbisP3Ready ( )

Definition at line 5007 of file HomePatch.C.

References ResizeArray< Elem >::begin(), ProxyGBISP3DataMsg::destPe, ProxyGBISP3DataMsg::dHdrPrefix, Patch::dHdrPrefix, ProxyGBISP3DataMsg::dHdrPrefixLen, Flags::doFullElectrostatics, Patch::flags, GB3_PROXY_DATA_PRIORITY, Patch::gbisP3Ready(), Patch::numAtoms, ProxyGBISP3DataMsg::origPe, ProxyGBISP3DataMsg::patch, PATCH_PRIORITY, Patch::patchID, PRIORITY_SIZE, Flags::sequence, SET_PRIORITY, and ResizeArray< Elem >::size().

Referenced by gbisComputeAfterP2().

                             {
   if (proxy.size() > 0) {
     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
     //only nonzero message should be sent for doFullElec
     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
     for (int i = 0; i < proxy.size(); i++) {
       int node = proxy[i];
       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
       msg->patch = patchID;
       msg->dHdrPrefixLen = msgAtoms;
       msg->origPe = CkMyPe();
       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
       msg->destPe = node;
       int seq = flags.sequence;
       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
       SET_PRIORITY(msg,seq,priority);
       cp[node].recvData(msg);
     }
   }
   Patch::gbisP3Ready();
 }

◆ getAtomList()

FullAtomList& HomePatch::getAtomList ( )

inline

Definition at line 528 of file HomePatch.h.

Referenced by ComputeGridForce::checkGridForceRatio(), ComputeGridForce::createGridForcedIdxList(), ComputeHomePatch::doWork(), dumpbench(), and ComputeLonepairsCUDA::updateAtoms().

528 { return (atom); }

◆ getMax()

ScaledPosition HomePatch::getMax ( )

inline

Definition at line 530 of file HomePatch.h.

530 {return max;}

◆ getMin()

ScaledPosition HomePatch::getMin ( )

inline

Definition at line 529 of file HomePatch.h.

529 {return min;}

◆ hardWallDrude()

int HomePatch::hardWallDrude	(	const BigReal	timestep,
		Tensor *	virial,
		SubmitReduction *	ppreduction
	)

Definition at line 3410 of file HomePatch.C.

References BOLTZMANN, Lattice::c(), endi(), iERROR(), iout, SubmitReduction::item(), Patch::lattice, Node::molecule, Patch::numAtoms, Node::Object(), Lattice::origin(), outer(), partition(), Node::simParameters, simParams, TIMEFACTOR, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz.

Referenced by Sequencer::hardWallDrude().

 {
   Molecule *mol = Node::Object()->molecule;
   SimParameters *simParams = Node::Object()->simParameters;
   const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
   const int fixedAtomsOn = simParams->fixedAtomsOn;
   const BigReal dt = timestep / TIMEFACTOR;
   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
   int i, ia, ib, j;
   int dieOnError = simParams->rigidDie;
   Tensor wc;  // constraint virial
   BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
   int nslabs;
 
   // start data for hard wall boundary between drude and its host atom
   // static int Count=0;
   int Idx;
   double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
   Vector v_ab, vb_1, vp_1, vb_2, vp_2, new_vel_a, new_vel_b, new_pos_a, new_pos_b, *new_pos, *new_vel;
   double dot_v_r_1, dot_v_r_2;
   double vb_cm, dr_a, dr_b;
   // end data for hard wall boundary between drude and its host atom
 
   // start calculation of hard wall boundary between drude and its host atom
   if (simParams->drudeHardWallOn) {
     if (ppreduction) {
       nslabs = simParams->pressureProfileSlabs;
       idz = nslabs/lattice.c().z;
       zmin = lattice.origin().z - 0.5*lattice.c().z;
     }
 
     r_wall = simParams->drudeBondLen;
     r_wall_SQ = r_wall*r_wall;
     // Count++;
     for (i=1; i<numAtoms; i++)  {
       if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) { // drude particle
         ia = i-1;
         ib = i;
 
         v_ab = atom[ib].position - atom[ia].position;
         rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
 
         if (rab_SQ > r_wall_SQ) {  // to impose the hard wall constraint
           rab = sqrt(rab_SQ);
           if ( (rab > (2.0*r_wall)) && dieOnError ) {  // unexpected situation
             iout << iERROR << "HardWallDrude> "
               << "The drude is too far away from atom "
               << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
             return -1;  // triggers early exit
           }
 
           v_ab.x /= rab;
           v_ab.y /= rab;
           v_ab.z /= rab;
 
           if ( fixedAtomsOn && atom[ia].atomFixed ) {  // the heavy atom is fixed
             if (atom[ib].atomFixed) {  // the drude is fixed too
               continue;
             }
             else {  // only the heavy atom is fixed
               dot_v_r_2 = atom[ib].velocity.x*v_ab.x
                 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
               vb_2 = dot_v_r_2 * v_ab;
               vp_2 = atom[ib].velocity - vb_2;
 
               dr = rab - r_wall;
               if(dot_v_r_2 == 0.0) {
                 delta_T = maxtime;
               }
               else {
                 delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
                 if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
               }
 
               dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
 
               vb_2 = dot_v_r_2 * v_ab;
 
               new_vel_a = atom[ia].velocity;
               new_vel_b = vp_2 + vb_2;
 
               dr_b = -dr + delta_T*dot_v_r_2;  // L = L_0 + dT *v_new, v was flipped
 
               new_pos_a = atom[ia].position;
               new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
             }
           }
           else {
             mass_a = atom[ia].mass;
             mass_b = atom[ib].mass;
             mass_sum = mass_a+mass_b;
 
             dot_v_r_1 = atom[ia].velocity.x*v_ab.x
               + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
             vb_1 = dot_v_r_1 * v_ab;
             vp_1 = atom[ia].velocity - vb_1;
 
             dot_v_r_2 = atom[ib].velocity.x*v_ab.x
               + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
             vb_2 = dot_v_r_2 * v_ab;
             vp_2 = atom[ib].velocity - vb_2;
 
             vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
 
             dot_v_r_1 -= vb_cm;
             dot_v_r_2 -= vb_cm;
 
             dr = rab - r_wall;
 
             if(dot_v_r_2 == dot_v_r_1) {
                 delta_T = maxtime;
             }
             else {
               delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);  // the time since the collision occurs
               if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
             }
             
             // the relative velocity between ia and ib. Drawn according to T_Drude
             v_Bond = sqrt(kbt/mass_b);
 
             // reflect the velocity along bond vector and scale down
             dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
             dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
 
             dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
             dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
 
             new_pos_a = atom[ia].position + dr_a*v_ab;  // correct the position
             new_pos_b = atom[ib].position + dr_b*v_ab;
             // atom[ia].position += (dr_a*v_ab);  // correct the position
             // atom[ib].position += (dr_b*v_ab);
             
             dot_v_r_1 += vb_cm;
             dot_v_r_2 += vb_cm;
 
             vb_1 = dot_v_r_1 * v_ab;
             vb_2 = dot_v_r_2 * v_ab;
 
             new_vel_a = vp_1 + vb_1;
             new_vel_b = vp_2 + vb_2;
           }
 
           int ppoffset, partition;
           if ( invdt == 0 ) {
             atom[ia].position = new_pos_a;
             atom[ib].position = new_pos_b;
           }
           else if ( virial == 0 ) {
             atom[ia].velocity = new_vel_a;
             atom[ib].velocity = new_vel_b;
           }
           else {
             for ( j = 0; j < 2; j++ ) {
               if (j==0) {  // atom ia, heavy atom
                 Idx = ia;
                 new_pos = &new_pos_a;
                 new_vel = &new_vel_a;
               }
               else if (j==1) {  // atom ib, drude
                 Idx = ib;
                 new_pos = &new_pos_b;
                 new_vel = &new_vel_b;
               }
               Force df = (*new_vel - atom[Idx].velocity) *
                 ( atom[Idx].mass * invdt );
               Tensor vir = outer(df, atom[Idx].position);
               wc += vir;
               atom[Idx].velocity = *new_vel;
               atom[Idx].position = *new_pos;
 
               if (ppreduction) {
                 if (!j) {
                   BigReal z = new_pos->z;
                   int partition = atom[Idx].partition;
                   int slab = (int)floor((z-zmin)*idz);
                   if (slab < 0) slab += nslabs;
                   else if (slab >= nslabs) slab -= nslabs;
                   ppoffset = 3*(slab + nslabs*partition);
                 }
                 ppreduction->item(ppoffset  ) += vir.xx;
                 ppreduction->item(ppoffset+1) += vir.yy;
                 ppreduction->item(ppoffset+2) += vir.zz;
               }
 
             }
           }
         }                               
       }
     }
 
     // if ( (Count>10000) && (Count%10==0) ) {
     //   v_ab = atom[1].position - atom[0].position;
     //   rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
     //   iout << "DBG_R: " << Count << "  " << sqrt(rab_SQ) << "\n" << endi;
     // }
 
   }
 
   // end calculation of hard wall boundary between drude and its host atom
 
   if ( dt && virial ) *virial += wc;
 
   return 0;
 }

◆ loweAndersenFinish()

void HomePatch::loweAndersenFinish ( )

Definition at line 4881 of file HomePatch.C.

References DebugM, ResizeArray< Elem >::resize(), and Patch::v.

Referenced by Sequencer::runComputeObjects().

 {
     DebugM(2, "loweAndersenFinish\n");
     v.resize(0);
 }

◆ loweAndersenVelocities()

void HomePatch::loweAndersenVelocities ( )

Definition at line 4866 of file HomePatch.C.

References DebugM, Node::molecule, Patch::numAtoms, Node::Object(), ResizeArray< Elem >::resize(), Node::simParameters, simParams, and Patch::v.

Referenced by positionsReady(), and positionsReady_SOA().

 {
     DebugM(2, "loweAndersenVelocities\n");
     Molecule *mol = Node::Object()->molecule;
     SimParameters *simParams = Node::Object()->simParameters;
     v.resize(numAtoms);
     for (int i = 0; i < numAtoms; ++i) {
         //v[i] = p[i];
         // co-opt CompAtom structure to pass velocity and mass information
         v[i].position = atom[i].velocity;
         v[i].charge = atom[i].mass;
     }
     DebugM(2, "loweAndersenVelocities\n");
 }

◆ minimize_rattle2()

void HomePatch::minimize_rattle2	(	const BigReal	timestep,
		Tensor *	virial,
		bool	forces = `false`
	)

Definition at line 4382 of file HomePatch.C.

References ResizeArray< Elem >::begin(), endi(), Patch::f, getWaterModelGroupSize(), iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Results::normal, Patch::numAtoms, Node::Object(), outer(), settle2(), Node::simParameters, simParams, SWM4, TIMEFACTOR, TIP4, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::newMinimizeDirection(), and Sequencer::submitMinimizeReductions().

 {
   Molecule *mol = Node::Object()->molecule;
   SimParameters *simParams = Node::Object()->simParameters;
   Force *f1 = f[Results::normal].begin();
   const int fixedAtomsOn = simParams->fixedAtomsOn;
   const int useSettle = simParams->useSettle;
   const BigReal dt = timestep / TIMEFACTOR;
   Tensor wc;  // constraint virial
   BigReal tol = simParams->rigidTol;
   int maxiter = simParams->rigidIter;
   int dieOnError = simParams->rigidDie;
   int i, iter;
   BigReal dsqi[10], tmp;
   int ial[10], ibl[10];
   Vector ref[10];  // reference position
   Vector refab[10];  // reference vector
   Vector vel[10];  // new velocity
   BigReal rmass[10];  // 1 / mass
   BigReal redmass[10];  // reduced mass
   int fixed[10];  // is atom fixed?
   
   // Size of a hydrogen group for water
   const WaterModel watmodel = simParams->watmodel;
   const int wathgsize = getWaterModelGroupSize(watmodel);
 
   //  CkPrintf("In rattle2!\n");
   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
     //    CkPrintf("ig=%d\n",ig);
     int hgs = atom[ig].hydrogenGroupSize;
     if ( hgs == 1 ) continue;  // only one atom in group
     // cache data in local arrays and integrate positions normally
     int anyfixed = 0;
     for ( i = 0; i < hgs; ++i ) {
       ref[i] = atom[ig+i].position;
       vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
       rmass[i] = 1.0;
       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
     }
     int icnt = 0;
     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
       if (hgs != wathgsize) {
         NAMD_bug("Hydrogen group error caught in rattle2().");
       }
       // Use SETTLE for water unless some of the water atoms are fixed,
       if (useSettle && !anyfixed) {
         if (watmodel == WaterModel::SWM4) {
           // SWM4 ordering:  O D LP H1 H2
           // do swap(O,LP) and call settle with subarray O H1 H2
           // swap back after we return
           Vector lp_ref = ref[2];
           // Vector lp_vel = vel[2];
           ref[2] = ref[0];
           vel[2] = vel[0];
           settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
           ref[0] = ref[2];
           vel[0] = vel[2];
           ref[2] = lp_ref;
           // vel[2] = vel[0];  // set LP vel to O vel
         }
         else {
           settle2(1.0, 1.0, ref, vel, dt, virial);
           if (watmodel == WaterModel::TIP4) {
             vel[3] = vel[0];
           }
         }
         for (i=0; i<hgs; i++) {
           ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
         }
         continue;
       }
       if ( !(fixed[1] && fixed[2]) ) {
         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
       }
     }
     //    CkPrintf("Loop 2\n");
     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
         if ( !(fixed[0] && fixed[i]) ) {
           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
           ibl[icnt] = i;  ++icnt;
         }
       }
     }
     if ( icnt == 0 ) continue;  // no constraints
     //    CkPrintf("Loop 3\n");
     for ( i = 0; i < icnt; ++i ) {
       refab[i] = ref[ial[i]] - ref[ibl[i]];
     }
     //    CkPrintf("Loop 4\n");
     int done;
     for ( iter = 0; iter < maxiter; ++iter ) {
       done = 1;
       for ( i = 0; i < icnt; ++i ) {
         int a = ial[i];  int b = ibl[i];
         Vector vab = vel[a] - vel[b];
         Vector &rab = refab[i];
         BigReal rabsqi = dsqi[i];
         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
           wc += outer(dp,rab);
           vel[a] += rmass[a] * dp;
           vel[b] -= rmass[b] * dp;
           done = 0;
         }
       }
       if ( done ) break;
       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
     }
     if ( ! done ) {
       if ( dieOnError ) {
         NAMD_die("Exceeded maximum number of iterations in rattle2().");
       } else {
         iout << iWARN <<
           "Exceeded maximum number of iterations in rattle2().\n" << endi;
       }
     }
     // store data back to patch
     for ( i = 0; i < hgs; ++i ) {
       ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
     }
   }
   //  CkPrintf("Leaving rattle2!\n");
   // check that there isn't a constant needed here!
   *virial += wc / ( 0.5 * dt );
 
 }

◆ mollyAverage()

void HomePatch::mollyAverage ( )

Definition at line 5085 of file HomePatch.C.

References average(), endi(), iout, iWARN(), Node::molecule, NAMD_die(), Patch::numAtoms, Node::Object(), Patch::p, Patch::p_avg, ResizeArray< Elem >::resize(), Node::simParameters, and simParams.

Referenced by positionsReady(), and positionsReady_SOA().

 {
   Molecule *mol = Node::Object()->molecule;
   SimParameters *simParams = Node::Object()->simParameters;
   BigReal tol = simParams->mollyTol;
   int maxiter = simParams->mollyIter;
   int i, iter;
   HGArrayBigReal dsq;
   BigReal tmp;
   HGArrayInt ial, ibl;
   HGArrayVector ref;  // reference position
   HGArrayVector refab;  // reference vector
   HGArrayBigReal rmass;  // 1 / mass
   BigReal *lambda;  // Lagrange multipliers
   CompAtom *avg;  // averaged position
   int numLambdas = 0;
   HGArrayInt fixed;  // is atom fixed?
 
   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
   p_avg.resize(numAtoms);
   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
 
   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
     int hgs = atom[ig].hydrogenGroupSize;
     if ( hgs == 1 ) continue;  // only one atom in group
         for ( i = 0; i < hgs; ++i ) {
           ref[i] = atom[ig+i].position;
           rmass[i] = 1. / atom[ig+i].mass;
           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
           if ( fixed[i] ) rmass[i] = 0.;
         }
         avg = &(p_avg[ig]);
         int icnt = 0;
 
   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
           if ( hgs != 3 ) {
             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
           }
           if ( !(fixed[1] && fixed[2]) ) {
             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
           }
         }
         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
             if ( !(fixed[0] && fixed[i]) ) {
               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
             }
           }
         }
         if ( icnt == 0 ) continue;  // no constraints
         numLambdas += icnt;
         molly_lambda.resize(numLambdas);
         lambda = &(molly_lambda[numLambdas - icnt]);
         for ( i = 0; i < icnt; ++i ) {
           refab[i] = ref[ial[i]] - ref[ibl[i]];
         }
         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
         if ( iter == maxiter ) {
           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
         }
   }
 
   // for ( i=0; i<numAtoms; ++i ) {
   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
   //    }
   // }
 
 }

◆ mollyMollify()

void HomePatch::mollyMollify ( Tensor * virial )

Definition at line 5159 of file HomePatch.C.

References Patch::f, Node::molecule, mollify(), NAMD_die(), Patch::numAtoms, Node::Object(), outer(), Patch::p_avg, ResizeArray< Elem >::resize(), Node::simParameters, simParams, and Results::slow.

Referenced by Sequencer::runComputeObjects().

 {
   Molecule *mol = Node::Object()->molecule;
   SimParameters *simParams = Node::Object()->simParameters;
   Tensor wc;  // constraint virial
   int i;
   HGArrayInt ial, ibl;
   HGArrayVector ref;  // reference position
   CompAtom *avg;  // averaged position
   HGArrayVector refab;  // reference vector
   HGArrayVector force;  // new force
   HGArrayBigReal rmass;  // 1 / mass
   BigReal *lambda;  // Lagrange multipliers
   int numLambdas = 0;
   HGArrayInt fixed;  // is atom fixed?
 
   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
     int hgs = atom[ig].hydrogenGroupSize;
     if (hgs == 1 ) continue;  // only one atom in group
         for ( i = 0; i < hgs; ++i ) {
           ref[i] = atom[ig+i].position;
           force[i] = f[Results::slow][ig+i];
           rmass[i] = 1. / atom[ig+i].mass;
           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
           if ( fixed[i] ) rmass[i] = 0.;
         }
         int icnt = 0;
         // c-ji I'm only going to mollify water for now
   if ( atom[ig].rigidBondLength > 0 ) {  // for water
           if ( hgs != 3 ) {
             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
           }
           if ( !(fixed[1] && fixed[2]) ) {
             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
           }
         }
         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
     if ( atom[ig+i].rigidBondLength > 0 ) {
             if ( !(fixed[0] && fixed[i]) ) {
               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
             }
           }
         }
 
         if ( icnt == 0 ) continue;  // no constraints
         lambda = &(molly_lambda[numLambdas]);
         numLambdas += icnt;
         for ( i = 0; i < icnt; ++i ) {
           refab[i] = ref[ial[i]] - ref[ibl[i]];
         }
         avg = &(p_avg[ig]);
         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
         // store data back to patch
         for ( i = 0; i < hgs; ++i ) {
           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
           f[Results::slow][ig+i] = force[i];
         }
   }
   // check that there isn't a constant needed here!
   *virial += wc;
   p_avg.resize(0);
 }

◆ positionsReady()

void HomePatch::positionsReady ( int doMigration = 0 )

Definition at line 1895 of file HomePatch.C.

References ProxyDataMsg::avgPlLen, ProxyDataMsg::avgPositionList, Patch::avgPositionPtrBegin, Patch::avgPositionPtrEnd, ResizeArray< Elem >::begin(), CompAtom::charge, COULOMB, ProxyDataMsg::cudaAtomList, Patch::cudaAtomPtr, DebugM, ComputeNonbondedUtil::dielectric_1, doAtomMigration(), Flags::doGBIS, doGroupSizeCheck(), Flags::doLCPO, Flags::doLoweAndersen, doMarginCheck(), Flags::doMolly, doPairlistCheck(), ResizeArray< Elem >::end(), ProxyDataMsg::flags, Patch::flags, Patch::getPatchID(), gridForceIdxChecked, Patch::intRad, ProxyDataMsg::intRadList, Patch::lattice, Patch::lcpoType, ProxyDataMsg::lcpoTypeList, loweAndersenVelocities(), mollyAverage(), NAMD_EVENT_START, NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NamdProfileEventStr, Patch::numAtoms, PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), PatchMap::Object(), Sync::Object(), Node::Object(), ProxyMgr::Object(), Patch::p, Patch::p_avg, ProxyDataMsg::patch, PATCH_PRIORITY, Patch::patchID, Sync::PatchReady(), Patch::pExt, ProxyDataMsg::plExtLen, ProxyDataMsg::plLen, CompAtom::position, ProxyDataMsg::positionExtList, ProxyDataMsg::positionList, Patch::positionsReady(), PRIORITY_SIZE, PROXY_DATA_PRIORITY, proxySendSpanning, proxySpanDim, qmSwapAtoms(), rattleListValid, ResizeArray< Elem >::resize(), ComputeNonbondedUtil::scaling, ProxyMgr::sendProxyAll(), ProxyMgr::sendProxyData(), Flags::sequence, SET_PRIORITY, setGBISIntrinsicRadii(), setLcpoType(), Node::simParameters, simParams, ResizeArray< Elem >::size(), sortAtomsForCUDA(), CompAtomExt::sortOrder, Lattice::unscale(), Patch::v, ProxyDataMsg::velocityList, Patch::velocityPtrBegin, Patch::velocityPtrEnd, ProxyDataMsg::vlLen, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::runComputeObjects().

 {    
   SimParameters *simParams = Node::Object()->simParameters;
 
   flags.sequence++;
   if (doMigration)
     gridForceIdxChecked=false;
   if (!patchMapRead) { readPatchMap(); }
       
   if (numNeighbors && ! simParams->staticAtomAssignment) {
     if (doMigration) {
       rattleListValid = false;
       doAtomMigration();
     } else {
       doMarginCheck();
     }
   }
 
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)  || defined(NAMD_ROCTX_ENABLED)
   char prbuf[32];
   sprintf(prbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::POSITIONS_READY], this->getPatchID());
   NAMD_EVENT_START_EX(1, NamdProfileEvent::POSITIONS_READY, prbuf);
 #endif
 
   if (doMigration && simParams->qmLSSOn)
       qmSwapAtoms();
 
 #if defined(NAMD_CUDA) || defined(NAMD_MIC) || defined(NAMD_AVXTILES) || defined(NAMD_HIP)
   #ifdef NAMD_AVXTILES
   if ( simParams->useAVXTiles ) {
   #endif
   if ( doMigration ) {
     int n = numAtoms;
     FullAtom *a_i = atom.begin();
     #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_AVXTILES) || \
       (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
     int *ao = new int[n];
     int nfree;
     if ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces ) {
       int k = 0;
       int k2 = n;
       for ( int j=0; j<n; ++j ) {
         // put fixed atoms at end
         if ( a_i[j].atomFixed ) ao[--k2] = j;
         else ao[k++] = j;
       }
       nfree = k;
     } else {
       nfree = n;
       for ( int j=0; j<n; ++j ) {
         ao[j] = j;
       }
     }
     sortAtomsForCUDA(ao,a_i,nfree,n);
     for ( int i=0; i<n; ++i ) { 
       a_i[i].sortOrder = ao[i];
     }
     delete [] ao;
     #else
       for (int i = 0; i < n; ++i) {
         a_i[i].sortOrder = i;
       }
     #endif
   }
 
   { 
     const double charge_scaling = sqrt(COULOMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1);
     const Vector ucenter = lattice.unscale(center);
     const BigReal ucenter_x = ucenter.x;
     const BigReal ucenter_y = ucenter.y;
     const BigReal ucenter_z = ucenter.z;
     const int n = numAtoms;
     #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
       int n_16 = n;
       n_16 = (n + 15) & (~15);
       cudaAtomList.resize(n_16);
       CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
       mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
       mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
       mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
       mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
     #elif defined(NAMD_AVXTILES)
     int n_avx = (n + 15) & (~15);
     cudaAtomList.resize(n_avx);
     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
     tiles.realloc(n, ac);
     #else
     #if 1
     cudaAtomList.resize(n);
     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
     #else
     reallocate_host<CudaAtom>(&cudaAtomList, &sizeCudaAtomList, n);
     CudaAtom *ac = cudaAtomPtr = &cudaAtomList[0];
     #endif
     #endif
     const FullAtom *a = atom.begin();
     for ( int k=0; k<n; ++k ) {
       #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
         int j = a[k].sortOrder;
         atom_x[k] = a[j].position.x - ucenter_x;
         atom_y[k] = a[j].position.y - ucenter_y;
         atom_z[k] = a[j].position.z - ucenter_z;
         atom_q[k] = charge_scaling * a[j].charge;
       #else
       int j = a[k].sortOrder;
       ac[k].x = a[j].position.x - ucenter_x;
       ac[k].y = a[j].position.y - ucenter_y;
       ac[k].z = a[j].position.z - ucenter_z;
       ac[k].q = charge_scaling * a[j].charge;
       #endif
     }
     #ifdef NAMD_AVXTILES
     {
       if (n > 0) {
         int j = a[n-1].sortOrder;
         for ( int k=n; k<n_avx; ++k ) {
           ac[k].x = a[j].position.x - ucenter_x;
           ac[k].y = a[j].position.y - ucenter_y;
           ac[k].z = a[j].position.z - ucenter_z;
         }
       }
     }
     #endif
   }
   #ifdef NAMD_AVXTILES
   // If "Tiles" mode disabled, no use of the CUDA data structures
   } else doMigration = doMigration && numNeighbors;
   #endif
 #else
   doMigration = doMigration && numNeighbors;
 #endif
   doMigration = doMigration || ! patchMapRead;
 
   doMigration = doMigration || doAtomUpdate;
   doAtomUpdate = false;
 
   // Workaround for oversize groups:
   // reset nonbondedGroupSize (ngs) before force calculation,
   // making sure that subset of hydrogen group starting with 
   // parent atom are all within 0.5 * hgroupCutoff.
   // XXX hydrogentGroupSize remains constant but is checked for nonzero
   // XXX should be skipped for CUDA, ngs not used by CUDA kernels
   // XXX should this also be skipped for KNL kernels?
   // ngs used by ComputeNonbondedBase.h - CPU nonbonded kernels
   // ngs used by ComputeGBIS.C - CPU GB nonbonded kernels
 #if ! (defined(NAMD_CUDA) || defined(NAMD_HIP))
 #if defined(NAMD_AVXTILES)
   if (!simParams->useAVXTiles)
 #endif
   doGroupSizeCheck();
 #endif
 
   // Copy information needed by computes and proxys to Patch::p.
   // Resize only if atoms were migrated
   if (doMigration) {
     p.resize(numAtoms);
     pExt.resize(numAtoms);
   }
   CompAtom *p_i = p.begin();
   CompAtomExt *pExt_i = pExt.begin();
   FullAtom *a_i = atom.begin();
   int i; int n = numAtoms;
   for ( i=0; i<n; ++i ) { 
     p_i[i] = a_i[i]; 
     pExt_i[i] = a_i[i];
   }
 
   // Measure atom movement to test pairlist validity
   doPairlistCheck();
 
   if (flags.doMolly) mollyAverage();
   // BEGIN LA
   if (flags.doLoweAndersen) loweAndersenVelocities();
   // END LA
 
     if (flags.doGBIS) {
       //reset for next time step
       numGBISP1Arrived = 0;
       phase1BoxClosedCalled = false;
       numGBISP2Arrived = 0;
       phase2BoxClosedCalled = false;
       numGBISP3Arrived = 0;
       phase3BoxClosedCalled = false;
       if (doMigration || isNewProxyAdded)
         setGBISIntrinsicRadii();
     }
 
   if (flags.doLCPO) {
     if (doMigration || isNewProxyAdded) {
       setLcpoType();
     }
   }
 
   // Must Add Proxy Changes when migration completed!
   NodeIDList::iterator pli;
   int *pids = NULL;
   int pidsPreAllocated = 1;
   int npid;
   if (proxySendSpanning == 0) {
     npid = proxy.size();
     pids = new int[npid];
     pidsPreAllocated = 0;
     int *pidi = pids;
     int *pide = pids + proxy.size();
     int patchNodesLast =
       ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
     for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
     {
       if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
         *(--pide) = *pli;
       } else {
         *(pidi++) = *pli;
       }
     }
   }
   else {
 #ifdef NODEAWARE_PROXY_SPANNINGTREE
     #ifdef USE_NODEPATCHMGR
     npid = numNodeChild;
     pids = nodeChildren;
     #else
     npid = nChild;
     pids = child;
     #endif
 #else
     npid = nChild;
     pidsPreAllocated = 0;
     pids = new int[proxySpanDim];
     for (int i=0; i<nChild; i++) pids[i] = child[i];
 #endif
   }
   if (npid) { //have proxies
     int seq = flags.sequence;
     int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
     //begin to prepare proxy msg and send it
     int pdMsgPLLen = p.size();
     int pdMsgAvgPLLen = 0;
     if(flags.doMolly) {
         pdMsgAvgPLLen = p_avg.size();
     }
     // BEGIN LA
     int pdMsgVLLen = 0;
     if (flags.doLoweAndersen) {
         pdMsgVLLen = v.size();
     }
     // END LA
 
     int intRadLen = 0;
     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
             intRadLen = numAtoms * 2;
     }
 
     //LCPO
     int lcpoTypeLen = 0;
     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
             lcpoTypeLen = numAtoms;
     }
 
     int pdMsgPLExtLen = 0;
     if(doMigration || isNewProxyAdded) {
         pdMsgPLExtLen = pExt.size();
     }
 
     int cudaAtomLen = 0;
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
     cudaAtomLen = numAtoms;
 #elif defined(NAMD_AVXTILES)
     if (simParams->useAVXTiles)
       cudaAtomLen = (numAtoms + 15) & (~15);
 #elif defined(NAMD_MIC)
     #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
       cudaAtomLen = numAtoms;
     #else
       cudaAtomLen = (numAtoms + 15) & (~15);
     #endif
 #endif
 
    ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
       lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
 
     SET_PRIORITY(nmsg,seq,priority);
     nmsg->patch = patchID;
     nmsg->flags = flags;
     nmsg->plLen = pdMsgPLLen;                
     //copying data to the newly created msg
     NAMD_EVENT_START(1, NamdProfileEvent::MEMCPY);
     memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
     NAMD_EVENT_STOP(1, NamdProfileEvent::MEMCPY);
     nmsg->avgPlLen = pdMsgAvgPLLen;        
     if(flags.doMolly) {
         memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
     }
     // BEGIN LA
     nmsg->vlLen = pdMsgVLLen;
     if (flags.doLoweAndersen) {
         memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
     }
     // END LA
 
     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
       for (int i = 0; i < numAtoms * 2; i++) {
         nmsg->intRadList[i] = intRad[i];
       }
     }
 
     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
       for (int i = 0; i < numAtoms; i++) {
         nmsg->lcpoTypeList[i] = lcpoType[i];
       }
     }
 
     nmsg->plExtLen = pdMsgPLExtLen;
     if(doMigration || isNewProxyAdded){     
         memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
     }
 
 // DMK
 #if defined(NAMD_CUDA) || defined(NAMD_MIC) || defined(NAMD_AVXTILES) || defined(NAMD_HIP)
     #ifdef NAMD_AVXTILES
     if (simParams->useAVXTiles)
     #endif
     {
     NAMD_EVENT_START(1, NamdProfileEvent::MEMCPY);
     memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
     NAMD_EVENT_STOP(1, NamdProfileEvent::MEMCPY);
     }
 #endif
     
 #if NAMD_SeparateWaters != 0
     //DMK - Atom Separation (water vs. non-water)
     nmsg->numWaterAtoms = numWaterAtoms;
 #endif
 
 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
     nmsg->isFromImmMsgCall = 0;
 #endif
     
     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
     DebugFileTrace *dft = DebugFileTrace::Object();
     dft->openTrace();
     dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
     for(int i=0; i<npid; i++) {
         dft->writeTrace("%d ", pids[i]);
     }
     dft->writeTrace("\n");
     dft->closeTrace();
     #endif
 
 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
     if (isProxyChanged || localphs == NULL)
     {
 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
      //CmiAssert(isProxyChanged);
      if (nphs) {
        for (int i=0; i<nphs; i++) {
          CmiDestoryPersistent(localphs[i]);
        }
        delete [] localphs;
      }
      localphs = new PersistentHandle[npid];
      int persist_size = sizeof(envelope) + sizeof(ProxyDataMsg) + sizeof(CompAtom)*(pdMsgPLLen+pdMsgAvgPLLen+pdMsgVLLen) + intRadLen*sizeof(Real) + lcpoTypeLen*sizeof(int) + sizeof(CompAtomExt)*pdMsgPLExtLen + sizeof(CudaAtom)*cudaAtomLen + PRIORITY_SIZE/8 + 2048;
      for (int i=0; i<npid; i++) {
 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
        if (proxySendSpanning)
            localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
        else
 #endif
        localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
      }
      nphs = npid;
     }
     CmiAssert(nphs == npid && localphs != NULL);
     CmiUsePersistentHandle(localphs, nphs);
 #endif
     if(doMigration || isNewProxyAdded) {
         ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
     }else{
         ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
     }
 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
     CmiUsePersistentHandle(NULL, 0);
 #endif
     isNewProxyAdded = 0;
   }
   isProxyChanged = 0;
   if(!pidsPreAllocated) delete [] pids;
   DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
 
 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
   positionPtrBegin = p.begin();
   positionPtrEnd = p.end();
 #endif
 
   if(flags.doMolly) {
       avgPositionPtrBegin = p_avg.begin();
       avgPositionPtrEnd = p_avg.end();
   }
   // BEGIN LA
   if (flags.doLoweAndersen) {
       velocityPtrBegin = v.begin();
       velocityPtrEnd = v.end();
   }
   // END LA
 
   Patch::positionsReady(doMigration);
 
   patchMapRead = 1;
 
   // gzheng
   Sync::Object()->PatchReady();
 
   NAMD_EVENT_STOP(1, NamdProfileEvent::POSITIONS_READY);
 
 }

◆ positionsReady_GPU()

void HomePatch::positionsReady_GPU	(	int	doMigration = `0`,
		int	startup = `0`
	)

◆ positionsReady_SOA()

void HomePatch::positionsReady_SOA ( int doMigration = 0 )

Definition at line 971 of file HomePatch.C.

References ASSERT, PatchDataSOA::atomFixed, ProxyDataMsg::avgPlLen, ProxyDataMsg::avgPositionList, Patch::avgPositionPtrBegin, Patch::avgPositionPtrEnd, ResizeArray< Elem >::begin(), PatchDataSOA::charge, COULOMB, ProxyDataMsg::cudaAtomList, Patch::cudaAtomPtr, DebugM, ComputeNonbondedUtil::dielectric_1, doAtomMigration(), Flags::doGBIS, doGroupSizeCheck_SOA(), Flags::doLCPO, Flags::doLoweAndersen, doMarginCheck_SOA(), Flags::doMolly, doPairlistCheck(), ResizeArray< Elem >::end(), PatchDataSOA::exclId, ProxyDataMsg::flags, Patch::flags, Patch::getPatchID(), PatchDataSOA::groupFixed, PatchDataSOA::hydrogenGroupSize, CompAtomExtCopy::id, PatchDataSOA::id, LocalID::index, Patch::intRad, ProxyDataMsg::intRadList, PatchDataSOA::isWater, Patch::lattice, Patch::lcpoType, ProxyDataMsg::lcpoTypeList, AtomMap::localID(), loweAndersenVelocities(), mollyAverage(), NAMD_ATOM_FIXED_MASK, NAMD_EVENT_START, NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NAMD_GROUP_FIXED_MASK, NamdProfileEventStr, PatchDataSOA::nonbondedGroupSize, Patch::numAtoms, PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), PatchMap::Object(), AtomMap::Object(), Sync::Object(), Node::Object(), ProxyMgr::Object(), Patch::p, Patch::p_avg, partition(), PatchDataSOA::partition, ProxyDataMsg::patch, PATCH_PRIORITY, Patch::patchID, Sync::PatchReady(), Patch::pExt, LocalID::pid, ProxyDataMsg::plExtLen, ProxyDataMsg::plLen, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, ProxyDataMsg::positionExtList, ProxyDataMsg::positionList, Patch::positionsReady(), PRIORITY_SIZE, PROXY_DATA_PRIORITY, proxySendSpanning, proxySpanDim, CudaAtom::q, qmSwapAtoms(), rattleListValid, rattleListValid_SOA, ResizeArray< Elem >::resize(), ComputeNonbondedUtil::scaling, ProxyMgr::sendProxyAll(), ProxyMgr::sendProxyData(), Flags::sequence, SET_PRIORITY, setGBISIntrinsicRadii(), setLcpoType(), PatchDataSOA::sigId, Node::simParameters, simParams, ResizeArray< Elem >::size(), sortAtomsForCUDA_SOA(), CompAtomExtCopy::sortOrder, PatchDataSOA::sortOrder, Lattice::unscale(), PatchDataSOA::unsortOrder, Patch::v, PatchDataSOA::vdwType, ProxyDataMsg::velocityList, Patch::velocityPtrBegin, Patch::velocityPtrEnd, ProxyDataMsg::vlLen, CudaAtom::x, Vector::x, CudaAtom::y, Vector::y, CudaAtom::z, and Vector::z.

Referenced by Sequencer::runComputeObjects_SOA().

 {
   SimParameters *simParams = Node::Object()->simParameters;
   // char prbuf[32];
   // sprintf(prbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::POSITIONS_READY_SOA], this->getPatchID());
   // NAMD_EVENT_START_EX(1, NamdProfileEvent::POSITIONS_READY_SOA, prbuf);
   if(!simParams->CUDASOAintegrate) flags.sequence++;
   // flags.sequence++;
   if (!patchMapRead) { readPatchMap(); }
   if (numNeighbors && ! simParams->staticAtomAssignment) {
     if (doMigration) {
       // copy SOA updates to AOS
       // XXX TODO: 
       copy_updates_to_AOS();
       // make sure to invalidate RATTLE lists when atoms move
       rattleListValid_SOA = false;
       rattleListValid = false;
       // this has a suspend
       doAtomMigration();
     } else {
       // XXX TODO: Get rid of this marginCheck for every tstep
       // move this to the GPU afterwards
       if(!simParams->CUDASOAintegrate) doMarginCheck_SOA();
     }
   }
 
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
   char prbuf[32];
   sprintf(prbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::POSITIONS_READY], this->getPatchID());
   NAMD_EVENT_START_EX(1, NamdProfileEvent::POSITIONS_READY, prbuf);
 #endif
 
 #if 0
   if (doMigration && simParams->qmLSSOn)
       qmSwapAtoms();
 #endif
 
 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
   if ( doMigration ) {
     // XXX Since we just migrated FullAtom array is up-to-date
     int n = numAtoms;
     int * __restrict ao_SOA = patchDataSOA.sortOrder;
     int * __restrict unsort_SOA = patchDataSOA.unsortOrder;
     #if defined(NAMD_CUDA) || defined(NAMD_HIP) || (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
     //#if 0
     int nfree_SOA;
     if ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces ) {
       int k = 0;
       int k2 = n;
       int * __restrict atomFixed = patchDataSOA.atomFixed;
       for ( int j=0; j<n; ++j ) {
         // put fixed atoms at end
         if ( atomFixed[j] ) ao_SOA[--k2] = j;
         else ao_SOA[k++] = j;
       }
       nfree_SOA = k;
     } else {
       nfree_SOA = n;
       for ( int j=0; j<n; ++j ) {
         ao_SOA[j] = j;
       }
     }
 #if 1
     sortAtomsForCUDA_SOA(ao_SOA, unsort_SOA,
         patchDataSOA.pos_x, patchDataSOA.pos_y, patchDataSOA.pos_z,
         nfree_SOA, n);
 #endif
 #if 0
     for (int i = 0; i < n; ++i) {
       ao_SOA[i] = i;
       printf("ao_SOA[%d] = %d\n", i, ao_SOA[i]);
     }
 #endif
 
 #else
     for (int i = 0; i < n; ++i) {
       ao_SOA[i] = i;
     }
 #endif
   }
   { 
     const double charge_scaling = sqrt(COULOMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1);
     const Vector ucenter = lattice.unscale(center);
     const BigReal ucenter_x = ucenter.x;
     const BigReal ucenter_y = ucenter.y;
     const BigReal ucenter_z = ucenter.z;
     const int n = numAtoms;
     #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
       int n_16 = n;
       n_16 = (n + 15) & (~15);
       cudaAtomList.resize(n_16);
       CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
       mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
       mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
       mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
       mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
     #elif defined(NAMD_AVXTILES)
     int n_avx = (n + 15) & (~15);
     cudaAtomList.resize(n_avx);
     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
     tiles.realloc(n, ac);
     #else
     if(doMigration) cudaAtomList.resize(n);
     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
     #endif
     double * __restrict pos_x = patchDataSOA.pos_x;
     double * __restrict pos_y = patchDataSOA.pos_y;
     double * __restrict pos_z = patchDataSOA.pos_z;
     float * __restrict charge = patchDataSOA.charge;
     int * __restrict ao_SOA = patchDataSOA.sortOrder;
 #ifndef NODEGROUP_FORCE_REGISTER
 //#if 1
     for ( int k=0; k<n; ++k ) {
       #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
         int j = ao_SOA[k];
         atom_x[k] = pos_x[j] - ucenter_x;
         atom_y[k] = pos_y[j] - ucenter_y;
         atom_z[k] = pos_z[j] - ucenter_z;
         atom_q[k] = charge_scaling * charge[j];
       #else
       // XXX TODO: This has to go
       int j = ao_SOA[k];
       // JM: Calculating single-precision patch-centered atomic coordinates as
       //     offsets. By adding single-precision offsets to double-precision 
       //     patch-patch center distances, we maintain full precision
       //     XXX NOTE: check where I can use this to use float instead of double in NAMD
       ac[k].x = pos_x[j] - ucenter_x;
       ac[k].y = pos_y[j] - ucenter_y;
       ac[k].z = pos_z[j] - ucenter_z;
       // XXX TODO: Compute charge scaling on GPUs and not here to avoid a copy 
       //           for every timestep
       // XXX TODO: Check when do we have to update this value and do it 
       //           on the gpu
       ac[k].q = charge_scaling * charge[j];
       #endif
     }
 #else
     if(!simParams->CUDASOAintegrate || doMigration){
       for ( int k=0; k<n; ++k ) {
         #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
           int j = ao_SOA[k];
           atom_x[k] = pos_x[j] - ucenter_x;
           atom_y[k] = pos_y[j] - ucenter_y;
           atom_z[k] = pos_z[j] - ucenter_z;
           atom_q[k] = charge_scaling * charge[j];
         #else
         // XXX TODO: This has to go
         int j = ao_SOA[k];
         // JM: Calculating single-precision patch-centered atomic coordinates as
         //     offsets. By adding single-precision offsets to double-precision 
         //     patch-patch center distances, we maintain full precision
         ac[k].x = pos_x[j] - ucenter_x;
         ac[k].y = pos_y[j] - ucenter_y;
         ac[k].z = pos_z[j] - ucenter_z;
         // XXX TODO: Compute charge scaling on GPUs and not here to avoid a copy 
         //           for every timestep
         // XXX TODO: Check when do we have to update this value and do it 
         //           on the gpu
         ac[k].q = charge_scaling * charge[j];
         #endif
       }
     }
 #endif
   }
 #else
   doMigration = doMigration && numNeighbors;
 #endif
   doMigration = doMigration || ! patchMapRead;
 
   doMigration = doMigration || doAtomUpdate;
   doAtomUpdate = false;
 
   // Workaround for oversize groups:
   // reset nonbondedGroupSize (ngs) before force calculation,
   // making sure that subset of hydrogen group starting with 
   // parent atom are all within 0.5 * hgroupCutoff.
   // XXX hydrogentGroupSize remains constant but is checked for nonzero
   // XXX should be skipped for CUDA, ngs not used by CUDA kernels
   // XXX should this also be skipped for KNL kernels?
   // ngs used by ComputeNonbondedBase.h - CPU nonbonded kernels
   // ngs used by ComputeGBIS.C - CPU GB nonbonded kernels
 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP)
   doGroupSizeCheck_SOA();
 #endif
 
   // Copy information needed by computes and proxys to Patch::p.
   if (doMigration) {
     // Atom data changes after migration, so must copy all of it.
     // Must resize CompAtom and CompAtomExt arrays to new atom count.
     p.resize(numAtoms);
     pExt.resize(numAtoms);
     CompAtom * __restrict p_i = p.begin();
     CompAtomExt * __restrict pExt_i = pExt.begin();
     const double * __restrict pos_x = patchDataSOA.pos_x;
     const double * __restrict pos_y = patchDataSOA.pos_y;
     const double * __restrict pos_z = patchDataSOA.pos_z;
     const float * __restrict charge = patchDataSOA.charge;
     const int * __restrict vdwType = patchDataSOA.vdwType;
     const int * __restrict partition = patchDataSOA.partition;
 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP)
     const int * __restrict nonbondedGroupSize = patchDataSOA.nonbondedGroupSize;
 #endif
     const int * __restrict hydrogenGroupSize = patchDataSOA.hydrogenGroupSize;
     const int * __restrict isWater = patchDataSOA.isWater;
     int n = numAtoms;
     for (int i=0;  i < n;  i++) { 
       p_i[i].position.x = pos_x[i];
       p_i[i].position.y = pos_y[i];
       p_i[i].position.z = pos_z[i];
       p_i[i].charge = charge[i];
       p_i[i].vdwType = vdwType[i];
       p_i[i].partition = partition[i];
 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP)
       p_i[i].nonbondedGroupSize = nonbondedGroupSize[i];
 #endif
       p_i[i].hydrogenGroupSize = hydrogenGroupSize[i];
       p_i[i].isWater = isWater[i];
     }
 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
     const int * __restrict sortOrder = patchDataSOA.sortOrder;
 #endif
     const int * __restrict id = patchDataSOA.id;
 #if defined(MEM_OPT_VERSION)
     const int * __restrict exclId = patchDataSOA.exclId;
     const int * __restrict sigId = patchDataSOA.sigId;
 #endif
     const int * __restrict atomFixed = patchDataSOA.atomFixed;
     const int * __restrict groupFixed = patchDataSOA.groupFixed;
     // Copy into CompAtomExt using typecast to temporary CompAtomExtCopy
     // to avoid loop vectorization bug in Intel 2018 compiler.
 #ifndef USE_NO_BITFIELDS
     CompAtomExtCopy *pExtCopy_i = (CompAtomExtCopy *) pExt_i;
 #endif // USE_NO_BITFIELDS
     for (int i=0;  i < n;  i++) {
 #ifndef USE_NO_BITFIELDS
 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
       pExtCopy_i[i].sortOrder = sortOrder[i];
 #endif
 #if defined(MEM_OPT_VERSION)
       pExtCopy_i[i].id = id[i];
       pExtCopy_i[i].exclId = exclId[i];
       ASSERT(atomFixed[i] == 0 || atomFixed[i] == 1);
       ASSERT(groupFixed[i] == 0 || groupFixed[i] == 1);
       uint32 atomFixedBit = NAMD_ATOM_FIXED_MASK * atomFixed[i];
       uint32 groupFixedBit = NAMD_GROUP_FIXED_MASK * groupFixed[i];
       pExtCopy_i[i].sigId = (sigId[i] | atomFixedBit | groupFixedBit);
 #else
       ASSERT(atomFixed[i] == 0 || atomFixed[i] == 1);
       ASSERT(groupFixed[i] == 0 || groupFixed[i] == 1);
       uint32 atomFixedBit = NAMD_ATOM_FIXED_MASK * atomFixed[i];
       uint32 groupFixedBit = NAMD_GROUP_FIXED_MASK * groupFixed[i];
       pExtCopy_i[i].id = (id[i] | atomFixedBit | groupFixedBit);
 #endif // if defined(MEM_OPT_VERSION)
 #else
 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
       pExt_i[i].sortOrder = sortOrder[i];
 #endif
       pExt_i[i].id = id[i];
 #if defined(MEM_OPT_VERSION)
       pExt_i[i].exclId = exclId[i];
       pExt_i[i].sigId = sigId[i];
 #endif
       pExt_i[i].atomFixed = atomFixed[i];
       pExt_i[i].groupFixed = groupFixed[i];
 #endif // USE_NO_BITFIELDS
     }
   }
   else {
     // JM: This is done for every timestep
     // Only need to copy positions, nonbondedGroupSize, and sortOrder.
     // Other data remains unchanged.
     CompAtom     * __restrict p_i = p.begin();
     CompAtomExt  * __restrict pExt_i = pExt.begin();
     const double * __restrict pos_x = patchDataSOA.pos_x;
     const double * __restrict pos_y = patchDataSOA.pos_y;
     const double * __restrict pos_z = patchDataSOA.pos_z;
 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP)
     const int * __restrict nonbondedGroupSize = patchDataSOA.nonbondedGroupSize;
 #endif
     int n = numAtoms;
 #ifndef NODEGROUP_FORCE_REGISTER
 //#if 1
     for (int i=0;  i < n;  i++) { 
       p_i[i].position.x = pos_x[i];
       p_i[i].position.y = pos_y[i];
       p_i[i].position.z = pos_z[i];
 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP)
       p_i[i].nonbondedGroupSize = nonbondedGroupSize[i];
 #endif
     }
 #else
     if(!simParams->CUDASOAintegrate || doMigration){
       for (int i=0;  i < n;  i++) { 
         p_i[i].position.x = pos_x[i];
         p_i[i].position.y = pos_y[i];
         p_i[i].position.z = pos_z[i];
 #if !defined(NAMD_CUDA) && !defined(NAMD_HIP)
         p_i[i].nonbondedGroupSize = nonbondedGroupSize[i];
 #endif
       }
     }
 #endif
 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
     const int * __restrict sortOrder = patchDataSOA.sortOrder;
 #ifdef NODEGROUP_FORCE_REGISTER
 //#if 0
     if(!simParams->CUDASOAintegrate || doMigration){
       for (int i=0;  i < n;  i++) { 
         pExt_i[i].sortOrder = sortOrder[i];
       }
     }
 #else
   //if(!simParams->CUDASOAintegrate || doMigration){
     for (int i=0;  i < n;  i++) { 
       pExt_i[i].sortOrder = sortOrder[i];
     }
 #endif
 #endif
   } // end Copy information to Patch::p
 
   // Measure atom movement to test pairlist validity
   // XXX TODO: Check if this ever needs to be done if CUDASOAintegrate
 #ifdef NODEGROUP_FORCE_REGISTER
   if(!simParams->CUDASOAintegrate || flags.sequence == 0){ 
     doPairlistCheck();
   }
 #else
   doPairlistCheck();
 #endif
 
 #if 0
   if (flags.doMolly) mollyAverage();
   // BEGIN LA
   if (flags.doLoweAndersen) loweAndersenVelocities();
   // END LA
 
     if (flags.doGBIS) {
       //reset for next time step
       numGBISP1Arrived = 0;
       phase1BoxClosedCalled = false;
       numGBISP2Arrived = 0;
       phase2BoxClosedCalled = false;
       numGBISP3Arrived = 0;
       phase3BoxClosedCalled = false;
       if (doMigration || isNewProxyAdded)
         setGBISIntrinsicRadii();
     }
 
   if (flags.doLCPO) {
     if (doMigration || isNewProxyAdded) {
       setLcpoType();
     }
   }
 #endif
 
   // Must Add Proxy Changes when migration completed!
   NodeIDList::iterator pli;
   int *pids = NULL;
   int pidsPreAllocated = 1;
   int npid;
   if (proxySendSpanning == 0) {
     npid = proxy.size();
     pids = new int[npid];
     pidsPreAllocated = 0;
     int *pidi = pids;
     int *pide = pids + proxy.size();
     int patchNodesLast =
       ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
     for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
     {
       if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
         *(--pide) = *pli;
       } else {
         *(pidi++) = *pli;
       }
     }
   }
   else {
 #ifdef NODEAWARE_PROXY_SPANNINGTREE
     #ifdef USE_NODEPATCHMGR
     npid = numNodeChild;
     pids = nodeChildren;
     #else
     npid = nChild;
     pids = child;
     #endif
 #else
     npid = nChild;
     pidsPreAllocated = 0;
     pids = new int[proxySpanDim];
     for (int i=0; i<nChild; i++) pids[i] = child[i];
 #endif
   }
   if (npid) { //have proxies
     int seq = flags.sequence;
     int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
     //begin to prepare proxy msg and send it
     int pdMsgPLLen = p.size();
     int pdMsgAvgPLLen = 0;
 #if 0
     if(flags.doMolly) {
         pdMsgAvgPLLen = p_avg.size();
     }
 #endif
     // BEGIN LA
     int pdMsgVLLen = 0;
 #if 0
     if (flags.doLoweAndersen) {
       pdMsgVLLen = v.size();
     }
 #endif
     // END LA
 
     int intRadLen = 0;
 #if 0
     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
             intRadLen = numAtoms * 2;
     }
 #endif
 
     //LCPO
     int lcpoTypeLen = 0;
 #if 0
     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
             lcpoTypeLen = numAtoms;
     }
 #endif
 
     int pdMsgPLExtLen = 0;
     if(doMigration || isNewProxyAdded) {
         pdMsgPLExtLen = pExt.size();
     }
 
     int cudaAtomLen = 0;
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
     cudaAtomLen = numAtoms;
 #endif
 
     #ifdef NAMD_MIC
       #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
         cudaAtomLen = numAtoms;
       #else
         cudaAtomLen = (numAtoms + 15) & (~15);
       #endif
     #endif
     ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
       lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
 
     SET_PRIORITY(nmsg,seq,priority);
     nmsg->patch = patchID;
     nmsg->flags = flags;
     nmsg->plLen = pdMsgPLLen;                
     //copying data to the newly created msg
     NAMD_EVENT_START(1, NamdProfileEvent::MEMCPY);
     memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
     NAMD_EVENT_STOP(1, NamdProfileEvent::MEMCPY);
     nmsg->avgPlLen = pdMsgAvgPLLen;        
 #if 0
     if(flags.doMolly) {
         memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
     }
 #endif
     // BEGIN LA
     nmsg->vlLen = pdMsgVLLen;
 #if 0
     if (flags.doLoweAndersen) {
         memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
     }
 #endif
     // END LA
 
 #if 0
     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
       for (int i = 0; i < numAtoms * 2; i++) {
         nmsg->intRadList[i] = intRad[i];
       }
     }
 #endif
 
 #if 0
     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
       for (int i = 0; i < numAtoms; i++) {
         nmsg->lcpoTypeList[i] = lcpoType[i];
       }
     }
 #endif
     nmsg->plExtLen = pdMsgPLExtLen;
     if(doMigration || isNewProxyAdded){     
         memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
     }
 
 // DMK
 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
     NAMD_EVENT_START(1, NamdProfileEvent::MEMCPY);
     memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
     NAMD_EVENT_STOP(1, NamdProfileEvent::MEMCPY);
 #endif
     
 #if NAMD_SeparateWaters != 0
     //DMK - Atom Separation (water vs. non-water)
     nmsg->numWaterAtoms = numWaterAtoms;
 #endif
 
 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
     nmsg->isFromImmMsgCall = 0;
 #endif
     
     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
     DebugFileTrace *dft = DebugFileTrace::Object();
     dft->openTrace();
     dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
     for(int i=0; i<npid; i++) {
         dft->writeTrace("%d ", pids[i]);
     }
     dft->writeTrace("\n");
     dft->closeTrace();
     #endif
 
 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
     if (isProxyChanged || localphs == NULL)
     {
 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
      //CmiAssert(isProxyChanged);
      if (nphs) {
        for (int i=0; i<nphs; i++) {
          CmiDestoryPersistent(localphs[i]);
        }
        delete [] localphs;
      }
      localphs = new PersistentHandle[npid];
      int persist_size = sizeof(envelope) + sizeof(ProxyDataMsg) + sizeof(CompAtom)*(pdMsgPLLen+pdMsgAvgPLLen+pdMsgVLLen) + intRadLen*sizeof(Real) + lcpoTypeLen*sizeof(int) + sizeof(CompAtomExt)*pdMsgPLExtLen + sizeof(CudaAtom)*cudaAtomLen + PRIORITY_SIZE/8 + 2048;
      for (int i=0; i<npid; i++) {
 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
        if (proxySendSpanning)
            localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
        else
 #endif
        localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
      }
      nphs = npid;
     }
     CmiAssert(nphs == npid && localphs != NULL);
     CmiUsePersistentHandle(localphs, nphs);
 #endif
     if(doMigration || isNewProxyAdded) {
         ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
     }else{
         ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
     }
 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
     CmiUsePersistentHandle(NULL, 0);
 #endif
     isNewProxyAdded = 0;
   }
   isProxyChanged = 0;
   if(!pidsPreAllocated) delete [] pids;
   DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
 
 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
   positionPtrBegin = p.begin();
   positionPtrEnd = p.end();
 #endif
 
 #if 0 
   if(flags.doMolly) {
       avgPositionPtrBegin = p_avg.begin();
       avgPositionPtrEnd = p_avg.end();
   }
   // BEGIN LA
   if (flags.doLoweAndersen) {
       velocityPtrBegin = v.begin();
       velocityPtrEnd = v.end();
   }
   // END LA
 #endif
   // fprintf(stderr, "(Pe[%d] tstep %d)calling positionsReady on patch %d\n", 
   //      CkMyPe(), this->flags.step, this->patchID);
   Patch::positionsReady(doMigration);
 
   patchMapRead = 1;
 
   // gzheng
   NAMD_EVENT_STOP(1, NamdProfileEvent::POSITIONS_READY_SOA);
 
   Sync::Object()->PatchReady();
 
 #if 0
   fprintf(stderr, "Patch %d atom IDS\n", this->patchID);
   AtomMap* mapper = AtomMap::Object();
   for(int i = 0 ; i < numAtoms; i++){
      fprintf(stderr, "atom[%d] = %d %d %d\n", i, atom[i].id, 
       mapper->localID(atom[i].id).pid, mapper->localID(atom[i].id).index);
   }
 #endif
 
 }

◆ qmSwapAtoms()

void HomePatch::qmSwapAtoms ( )

Definition at line 907 of file HomePatch.C.

References ResizeArray< Elem >::begin(), CompAtom::charge, SortedArray< Elem >::find(), Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), CompAtomExt::id, lssSubs(), Node::molecule, LSSSubsDat::newCharge, LSSSubsDat::newID, LSSSubsDat::newVdWType, Patch::numAtoms, Node::Object(), Node::simParameters, simParams, and CompAtom::vdwType.

Referenced by positionsReady(), and positionsReady_SOA().

 {
     // This is used for LSS in QM/MM simulations.
     // Changes atom labels so that we effectively exchange solvent
     // molecules between classical and quantum modes.
     
     SimParameters *simParams = Node::Object()->simParameters;
     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
     
     ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
             CkpvAccess(BOCclass_group).computeQMMgr) ;
     
     FullAtom *a_i = atom.begin();
     
     for (int i=0; i<numAtoms; ++i ) { 
         
         LSSSubsDat *subP = lssSubs(mgrP).find( LSSSubsDat(a_i[i].id) ) ;
         
         if ( subP != NULL ) {
             a_i[i].id = subP->newID ;
             a_i[i].vdwType = subP->newVdWType ;
             
             // If we are swappign a classical atom with a QM one, the charge
             // will need extra handling.
             if (qmAtomGroup[subP->newID] > 0 && simParams->PMEOn) {
                 // We make sure that the last atom charge calculated for the 
                 // QM atom being transfered here is available for PME
                 // in the next step.
                 
                 // Loops over all QM atoms (in all QM groups) comparing their 
                 // global indices (the sequential atom ID from NAMD).
                 for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
                     
                     if (qmAtmIndx[qmIter] == subP->newID) {
                         qmAtmChrg[qmIter] = subP->newCharge;
                         break;
                     }
                     
                 }
                 
                 // For QM atoms, the charge in the full atom structure is zero.
                 // Electrostatic interactions between QM atoms and their 
                 // environment is handled in ComputeQM.
                 a_i[i].charge = 0;
             }
             else {
                 // If we are swappign a QM atom with a Classical one, only the charge
                 // in the full atomstructure needs updating, since it used to be zero.
                 a_i[i].charge = subP->newCharge ;
             }
         }
     }
     
     return;
 }

◆ rattle1()

int HomePatch::rattle1	(	const BigReal	timestep,
		Tensor *	virial,
		SubmitReduction *	ppreduction
	)

Definition at line 3788 of file HomePatch.C.

References addRattleForce(), buildRattleList(), endi(), iERROR(), iout, iWARN(), LINCS(), MSHAKEIterate(), noconstList, Patch::numAtoms, Node::Object(), posNew, rattle1old(), rattleList, rattleListValid, rattleN(), rattlePair< 1 >(), rattleParam, settle1_SIMD< 1 >(), settle1_SIMD< 2 >(), settleList, Node::simParameters, simParams, TIMEFACTOR, TIP3, velNew, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::minimize(), Sequencer::minimizeMoveDownhill(), Sequencer::newMinimizePosition(), and Sequencer::rattle1().

                                 {
 
   SimParameters *simParams = Node::Object()->simParameters;
   if (simParams->watmodel != WaterModel::TIP3 || ppreduction) {
     // Call old rattle1 -method instead
     return rattle1old(timestep, virial, ppreduction);
   }
 
   if (!rattleListValid) {
     buildRattleList();
     rattleListValid = true;
   }
 
   const int fixedAtomsOn = simParams->fixedAtomsOn;
   const int useSettle = simParams->useSettle;
   const BigReal dt = timestep / TIMEFACTOR;
   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
   const BigReal tol2 = 2.0 * simParams->rigidTol;
   int maxiter = simParams->rigidIter;
   int dieOnError = simParams->rigidDie;
 
   Vector ref[10];  // reference position
   Vector pos[10];  // new position
   Vector vel[10];  // new velocity
 
   // Manual un-roll
   int n = (settleList.size()/2)*2;
   for (int j=0;j < n;j+=2) {
     int ig;
     ig = settleList[j];
     for (int i = 0; i < 3; ++i ) {
       ref[i] = atom[ig+i].position;
       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
     }
     ig = settleList[j+1];
     for (int i = 0; i < 3; ++i ) {
       ref[i+3] = atom[ig+i].position;
       pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
     }
     settle1_SIMD<2>(ref, pos,
       settle_mOrmT, settle_mHrmT, settle_ra,
       settle_rb, settle_rc, settle_rra);
 
     ig = settleList[j];
     for (int i = 0; i < 3; ++i ) {
       velNew[ig+i] = (pos[i] - ref[i])*invdt;
       posNew[ig+i] = pos[i];
     }
     ig = settleList[j+1];
     for (int i = 0; i < 3; ++i ) {
       velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
       posNew[ig+i] = pos[i+3];
     }
 
   }
 
   if (settleList.size() % 2) {
     int ig = settleList[settleList.size()-1];
     for (int i = 0; i < 3; ++i ) {
       ref[i] = atom[ig+i].position;
       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
     }
     settle1_SIMD<1>(ref, pos,
             settle_mOrmT, settle_mHrmT, settle_ra,
             settle_rb, settle_rc, settle_rra);        
     for (int i = 0; i < 3; ++i ) {
       velNew[ig+i] = (pos[i] - ref[i])*invdt;
       posNew[ig+i] = pos[i];
     }
   }
 
   int posParam = 0;
   for (int j=0;j < rattleList.size();++j) {
 
     BigReal refx[10];
     BigReal refy[10];
     BigReal refz[10];
 
     BigReal posx[10];
     BigReal posy[10];
     BigReal posz[10];
 
     int ig = rattleList[j].ig;
     int icnt = rattleList[j].icnt;
     int hgs = atom[ig].hydrogenGroupSize;
     for (int i = 0; i < hgs; ++i ) {
       ref[i] = atom[ig+i].position;
       pos[i] = atom[ig+i].position;
       if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
         pos[i] += atom[ig+i].velocity * dt;
       }
       refx[i] = ref[i].x;
       refy[i] = ref[i].y;
       refz[i] = ref[i].z;
       posx[i] = pos[i].x;
       posy[i] = pos[i].y;
       posz[i] = pos[i].z;
     }
 
     bool done;
     bool consFailure;
     if (icnt == 1) {
       rattlePair<1>(&rattleParam[posParam],
         refx, refy, refz,
         posx, posy, posz,
         consFailure);
       done = true;
     } else {
       if (simParams->mshakeOn) {
        MSHAKEIterate(icnt, &rattleParam[posParam],
             refx, refy, refz,
             posx, posy, posz,
             tol2, maxiter,
             done, consFailure);
       }
       else if(simParams->lincsOn) {
                             LINCS(icnt, &rattleParam[posParam],
                             refx, refy, refz,
                             posx, posy, posz,
                             tol2, maxiter, done, consFailure);
       }
       else
         rattleN(icnt, &rattleParam[posParam],
             refx, refy, refz,
             posx, posy, posz,
             tol2, maxiter,
             done, consFailure);
     }
 
     // Advance position in rattleParam
     //posParam += icnt;
     posParam += 4;
     for (int i = 0; i < hgs; ++i ) {
       pos[i].x = posx[i];
       pos[i].y = posy[i];
       pos[i].z = posz[i];
     }
 
     for (int i = 0; i < hgs; ++i ) {
       velNew[ig+i] = (pos[i] - ref[i])*invdt;
       posNew[ig+i] = pos[i];
     }
 
     if ( consFailure ) {
       if ( dieOnError ) {
         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
         return -1;  // triggers early exit
       } else {
         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
       }
     } else if ( ! done ) {
       if ( dieOnError ) {
         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
         return -1;  // triggers early exit
       } else {
         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
       }
     }
   }
   // Finally, we have to go through atoms that are not involved in rattle just so that we have
   // their positions and velocities up-to-date in posNew and velNew
   for (int j=0;j < noconstList.size();++j) {
     int ig = noconstList[j];
     int hgs = atom[ig].hydrogenGroupSize;
     for (int i = 0; i < hgs; ++i ) {
       velNew[ig+i] = atom[ig+i].velocity;
       posNew[ig+i] = atom[ig+i].position;
     }
   }
 
   if ( invdt == 0 ) {
     for (int ig = 0; ig < numAtoms; ++ig )
       atom[ig].position = posNew[ig];
   } else if ( virial == 0 ) {
     for (int ig = 0; ig < numAtoms; ++ig )
       atom[ig].velocity = velNew[ig];
   } else {
     Tensor wc;  // constraint virial
     addRattleForce(invdt, wc);
     *virial += wc;
   }
 
   return 0;
 }

◆ rattle1_SOA()

int HomePatch::rattle1_SOA	(	const BigReal	dt,
		Tensor *	virial,
		SubmitReduction *	ppreduction
	)

Definition at line 4659 of file HomePatch.C.

References buildRattleList_SOA(), endi(), PatchDataSOA::f_normal_x, PatchDataSOA::f_normal_y, PatchDataSOA::f_normal_z, PatchDataSOA::hydrogenGroupSize, iERROR(), iout, iWARN(), LINCS(), PatchDataSOA::mass, MSHAKEIterate(), noconstList, Patch::numAtoms, Node::Object(), PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, PatchDataSOA::posNew_x, PatchDataSOA::posNew_y, PatchDataSOA::posNew_z, rattleList, rattleListValid_SOA, rattleN(), rattlePair< 1 >(), rattleParam, settle1_SOA(), Node::simParameters, simParams, PatchDataSOA::vel_x, PatchDataSOA::vel_y, PatchDataSOA::vel_z, PatchDataSOA::velNew_x, PatchDataSOA::velNew_y, PatchDataSOA::velNew_z, Tensor::xx, Tensor::xy, Tensor::xz, Tensor::yx, Tensor::yy, Tensor::yz, Tensor::zx, Tensor::zy, and Tensor::zz.

Referenced by Sequencer::rattle1_SOA().

                                 {
 
   SimParameters *simParams = Node::Object()->simParameters;
 
 #if 1
   if (!rattleListValid_SOA) {
     buildRattleList_SOA();
     rattleListValid_SOA = true;
   }
 #endif
 
   double * __restrict pos_x = patchDataSOA.pos_x;
   double * __restrict pos_y = patchDataSOA.pos_y;
   double * __restrict pos_z = patchDataSOA.pos_z;
   double * __restrict vel_x = patchDataSOA.vel_x;
   double * __restrict vel_y = patchDataSOA.vel_y;
   double * __restrict vel_z = patchDataSOA.vel_z;
   double * __restrict posNew_x = patchDataSOA.posNew_x;
   double * __restrict posNew_y = patchDataSOA.posNew_y;
   double * __restrict posNew_z = patchDataSOA.posNew_z;
   double * __restrict velNew_x = patchDataSOA.velNew_x;
   double * __restrict velNew_y = patchDataSOA.velNew_y;
   double * __restrict velNew_z = patchDataSOA.velNew_z;
   const int * __restrict hydrogenGroupSize = patchDataSOA.hydrogenGroupSize;
 #ifdef __INTEL_COMPILER
   __assume_aligned(pos_x,64);
   __assume_aligned(pos_y,64);
   __assume_aligned(pos_z,64);
   __assume_aligned(vel_x,64);
   __assume_aligned(vel_y,64);
   __assume_aligned(vel_z,64);
   __assume_aligned(posNew_x,64);
   __assume_aligned(posNew_y,64);
   __assume_aligned(posNew_z,64);
   __assume_aligned(velNew_x,64);
   __assume_aligned(velNew_y,64);
   __assume_aligned(velNew_z,64);
   __assume_aligned(hydrogenGroupSize,64);
 #endif
 
   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
   const BigReal tol2 = 2.0 * simParams->rigidTol;
   int maxiter = simParams->rigidIter;
   int dieOnError = simParams->rigidDie;
 
   // calculate full step update to all positions
   for (int i=0;  i < numAtoms;  i++) {
     posNew_x[i] = pos_x[i] + vel_x[i] * dt;
     posNew_y[i] = pos_y[i] + vel_y[i] * dt;
     posNew_z[i] = pos_z[i] + vel_z[i] * dt;
   }
 
   // call settle to process all waters at once
   // XXX this assumes sorting waters into consecutive part of list
   //int numWaters = settleList.size();
   if (numSolventAtoms > 0) {
     int n = numSoluteAtoms;  // index of first water in list is past solute
     settle1_SOA(&pos_x[n], &pos_y[n], &pos_z[n],
         &posNew_x[n], &posNew_y[n], &posNew_z[n],
         numWaters,
         settle_mOrmT, settle_mHrmT, settle_ra,
         settle_rb, settle_rc, settle_rra);        
   }
   int posParam = 0;
   for (int j=0;j < rattleList.size();++j) {
     int ig = rattleList[j].ig;
     int icnt = rattleList[j].icnt;
     bool done;
     bool consFailure;
     if (icnt == 1) {
       rattlePair<1>(&rattleParam[posParam],
         &pos_x[ig], &pos_y[ig], &pos_z[ig],
         &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
         consFailure
         );
       done = true;
     } else {
       if (simParams->mshakeOn) {
         //buildConstantMatrix();
         MSHAKEIterate(icnt, &rattleParam[posParam],
             &pos_x[ig], &pos_y[ig], &pos_z[ig],
             &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
             tol2, maxiter,
             done, consFailure);
       }
       else if(simParams->lincsOn) {
                             LINCS(icnt, &rattleParam[posParam],
                             &pos_x[ig], &pos_y[ig], &pos_z[ig],
                             &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
                             tol2, maxiter, done, consFailure);
       }
 
       else
         rattleN(icnt, &rattleParam[posParam],
             &pos_x[ig], &pos_y[ig], &pos_z[ig],
             &posNew_x[ig], &posNew_y[ig], &posNew_z[ig],
             tol2, maxiter,
             done, consFailure);
     }
 
     // Advance position in rattleParam
 //    posParam += icnt;
     posParam += 4;
     if ( consFailure ) {
       if ( dieOnError ) {
         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
         return -1;  // triggers early exit
       } else {
         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
       }
     } else if ( ! done ) {
       if ( dieOnError ) {
         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
         return -1;  // triggers early exit
       } else {
         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
       }
     }
 
   } // end rattle
   // Now that all new positions are known, determine new velocities
   // needed to reach new position.
   for (int i=0;  i < numAtoms;  i++) {
     velNew_x[i] = (posNew_x[i] - pos_x[i]) * invdt;
     velNew_y[i] = (posNew_y[i] - pos_y[i]) * invdt;
     velNew_z[i] = (posNew_z[i] - pos_z[i]) * invdt;
   }
 
   // Bring new positions and velocities back to reference for noconstList.
   // No need to check hydrogen group size, since no fixed atoms.
   int numNoconst = noconstList.size();
   for (int j=0;  j < numNoconst;  j++) {
     int ig = noconstList[j];
     posNew_x[ig] = pos_x[ig];
     posNew_y[ig] = pos_y[ig];
     posNew_z[ig] = pos_z[ig];
     velNew_x[ig] = vel_x[ig];
     velNew_y[ig] = vel_y[ig];
     velNew_z[ig] = vel_z[ig];
   }
 
   if ( invdt == 0 ) {
     for (int ig = 0; ig < numAtoms; ++ig ) {
       pos_x[ig] = posNew_x[ig];
       pos_y[ig] = posNew_y[ig];
       pos_z[ig] = posNew_z[ig];
     }
   }
   else if ( virial == 0 ) {
     for (int ig = 0; ig < numAtoms; ++ig ) {
       vel_x[ig] = velNew_x[ig];
       vel_y[ig] = velNew_y[ig];
       vel_z[ig] = velNew_z[ig];
     }
   }
   else {
     const float  * __restrict mass = patchDataSOA.mass;
     double * __restrict f_normal_x = patchDataSOA.f_normal_x;
     double * __restrict f_normal_y = patchDataSOA.f_normal_y;
     double * __restrict f_normal_z = patchDataSOA.f_normal_z;
 #ifdef __INTEL_COMPILER
     __assume_aligned(mass,64);
     __assume_aligned(f_normal_x,64);
     __assume_aligned(f_normal_y,64);
     __assume_aligned(f_normal_z,64);
 #endif
     Tensor vir;  // = 0
     for (int ig = 0;  ig < numAtoms;  ig++) {
       BigReal df_x = (velNew_x[ig] - vel_x[ig]) * ( mass[ig] * invdt );
       BigReal df_y = (velNew_y[ig] - vel_y[ig]) * ( mass[ig] * invdt );
       BigReal df_z = (velNew_z[ig] - vel_z[ig]) * ( mass[ig] * invdt );
       f_normal_x[ig] += df_x;
       f_normal_y[ig] += df_y;
       f_normal_z[ig] += df_z;
       vir.xx += df_x * pos_x[ig];
       vir.xy += df_x * pos_y[ig];
       vir.xz += df_x * pos_z[ig];
       vir.yx += df_y * pos_x[ig];
       vir.yy += df_y * pos_y[ig];
       vir.yz += df_y * pos_z[ig];
       vir.zx += df_z * pos_x[ig];
       vir.zy += df_z * pos_y[ig];
       vir.zz += df_z * pos_z[ig];
     }
     *virial += vir;
     for (int ig = 0;  ig < numAtoms;  ig++) {
       vel_x[ig] = velNew_x[ig];
       vel_y[ig] = velNew_y[ig];
       vel_z[ig] = velNew_z[ig];
     }
   }
 
   return 0;
 }

◆ rattle1old()

int HomePatch::rattle1old	(	const BigReal	timestep,
		Tensor *	virial,
		SubmitReduction *	ppreduction
	)

Definition at line 3979 of file HomePatch.C.

References Lattice::c(), endi(), Patch::f, getWaterModelGroupSize(), iERROR(), iout, SubmitReduction::item(), iWARN(), Patch::lattice, Node::molecule, NAMD_die(), Results::normal, Patch::numAtoms, Node::Object(), Lattice::origin(), outer(), partition(), settle1(), settle1init(), Node::simParameters, simParams, SWM4, TIMEFACTOR, TIP4, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz.

Referenced by rattle1().

 {
   // CkPrintf("Call HomePatch::rattle1old\n");
   Molecule *mol = Node::Object()->molecule;
   SimParameters *simParams = Node::Object()->simParameters;
   const int fixedAtomsOn = simParams->fixedAtomsOn;
   const int useSettle = simParams->useSettle;
   const BigReal dt = timestep / TIMEFACTOR;
   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
   BigReal tol2 = 2.0 * simParams->rigidTol;
   int maxiter = simParams->rigidIter;
   int dieOnError = simParams->rigidDie;
   int i, iter;
   BigReal dsq[10], tmp;
   int ial[10], ibl[10];
   Vector ref[10];  // reference position
   Vector refab[10];  // reference vector
   Vector pos[10];  // new position
   Vector vel[10];  // new velocity
   Vector netdp[10];  // total momentum change from constraint
   BigReal rmass[10];  // 1 / mass
   int fixed[10];  // is atom fixed?
   Tensor wc;  // constraint virial
   BigReal idz, zmin;
   int nslabs;
 
   // Size of a hydrogen group for water
   const WaterModel watmodel = simParams->watmodel;
   const int wathgsize = getWaterModelGroupSize(watmodel);
 
   // Initialize the settle algorithm with water parameters
   // settle1() assumes all waters are identical,
   // and will generate bad results if they are not.
   // XXX this will move to Molecule::build_atom_status when that 
   // version is debugged
   if ( ! settle_initialized ) {
     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
       // find a water
       if (atom[ig].rigidBondLength > 0) {
         int oatm;
         if (watmodel == WaterModel::SWM4) {
           oatm = ig+3;  // skip over Drude and Lonepair
           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
         }
         else {
           oatm = ig+1;
           // Avoid using the Om site to set this by mistake
           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
             oatm += 1;
           }
         }
 
         // initialize settle water parameters
         settle1init(atom[ig].mass, atom[oatm].mass, 
                     atom[ig].rigidBondLength,
                     atom[oatm].rigidBondLength,
                     settle_mO, settle_mH,
                     settle_mOrmT, settle_mHrmT, settle_ra,
                     settle_rb, settle_rc, settle_rra);
         settle_initialized = 1;
         break; // done with init
       }
     }
   }
 
   if (ppreduction) {
     nslabs = simParams->pressureProfileSlabs;
     idz = nslabs/lattice.c().z;
     zmin = lattice.origin().z - 0.5*lattice.c().z;
   }
   
   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
     int hgs = atom[ig].hydrogenGroupSize;
     if ( hgs == 1 ) continue;  // only one atom in group
     // cache data in local arrays and integrate positions normally
     int anyfixed = 0;
     for ( i = 0; i < hgs; ++i ) {
       ref[i] = atom[ig+i].position;
       pos[i] = atom[ig+i].position;
       vel[i] = atom[ig+i].velocity;
       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
       //printf("rmass of %i is %f\n", ig+i, rmass[i]);
       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
       //printf("fixed status of %i is %i\n", i, fixed[i]);
       // undo addVelocityToPosition to get proper reference coordinates
       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
     }
     int icnt = 0;
     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
       if (hgs != wathgsize) {
         char errmsg[256];
         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
                          "but the specified water model requires %d atoms.\n",
                          atom[ig].id+1, hgs, wathgsize);
         NAMD_die(errmsg);
       }
       // Use SETTLE for water unless some of the water atoms are fixed,
       if (useSettle && !anyfixed) {
         if (watmodel == WaterModel::SWM4) {
           // SWM4 ordering:  O D LP H1 H2
           // do swap(O,LP) and call settle with subarray O H1 H2
           // swap back after we return
           Vector lp_ref = ref[2];
           Vector lp_pos = pos[2];
           Vector lp_vel = vel[2];
           ref[2] = ref[0];
           pos[2] = pos[0];
           vel[2] = vel[0];
           settle1(ref+2, pos+2, vel+2, invdt,
                   settle_mOrmT, settle_mHrmT, settle_ra,
                   settle_rb, settle_rc, settle_rra);
           ref[0] = ref[2];
           pos[0] = pos[2];
           vel[0] = vel[2];
           ref[2] = lp_ref;
           pos[2] = lp_pos;
           vel[2] = lp_vel;
           // determine for LP updated pos and vel
           swm4_omrepos(ref, pos, vel, invdt);
         }
         else {
           settle1(ref, pos, vel, invdt,
                   settle_mOrmT, settle_mHrmT, settle_ra,
                   settle_rb, settle_rc, settle_rra);
           if (watmodel == WaterModel::TIP4) {
             tip4_omrepos(ref, pos, vel, invdt);
           }
         }
 
         // which slab the hydrogen group will belong to
         // for pprofile calculations.
         int ppoffset, partition;
         if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
           atom[ig+i].position = pos[i];
         } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
           atom[ig+i].velocity = vel[i];
         } else for ( i = 0; i < wathgsize; ++i ) {
           Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
           Tensor vir = outer(df, ref[i]);
           wc += vir;
           f[Results::normal][ig+i] += df;
           atom[ig+i].velocity = vel[i];
           if (ppreduction) {
             // put all the atoms from a water in the same slab.  Atom 0
             // should be the parent atom.
             if (!i) {
               BigReal z = pos[i].z;
               partition = atom[ig].partition;
               int slab = (int)floor((z-zmin)*idz);
               if (slab < 0) slab += nslabs;
               else if (slab >= nslabs) slab -= nslabs;
               ppoffset = 3*(slab + nslabs*partition);
             }
             ppreduction->item(ppoffset  ) += vir.xx;
             ppreduction->item(ppoffset+1) += vir.yy;
             ppreduction->item(ppoffset+2) += vir.zz;
           }
         }
         continue;
       }
       if ( !(fixed[1] && fixed[2]) ) {
         dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
       }
     }
     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
         if ( !(fixed[0] && fixed[i]) ) {
           dsq[icnt] = tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
         }
       }
     }
     if ( icnt == 0 ) continue;  // no constraints
     for ( i = 0; i < icnt; ++i ) {
       refab[i] = ref[ial[i]] - ref[ibl[i]];
     }
     for ( i = 0; i < hgs; ++i ) {
       netdp[i] = 0.;
     }
     int done;
     int consFailure;
     for ( iter = 0; iter < maxiter; ++iter ) {
 //if (iter > 0) CkPrintf("iteration %d\n", iter);
       done = 1;
       consFailure = 0;
       for ( i = 0; i < icnt; ++i ) {
         int a = ial[i];  int b = ibl[i];
         Vector pab = pos[a] - pos[b];
               BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
         BigReal rabsq = dsq[i];
         BigReal diffsq = rabsq - pabsq;
         if ( fabs(diffsq) > (rabsq * tol2) ) {
           Vector &rab = refab[i];
           BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
           if ( rpab < ( rabsq * 1.0e-6 ) ) {
             done = 0;
             consFailure = 1;
             continue;
           }
           BigReal rma = rmass[a];
           BigReal rmb = rmass[b];
           BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
           Vector dp = rab * gab;
           pos[a] += rma * dp;
           pos[b] -= rmb * dp;
           if ( invdt != 0. ) {
             dp *= invdt;
             netdp[a] += dp;
             netdp[b] -= dp;
           }
           done = 0;
         }
       }
       if ( done ) break;
     }
 
     if ( consFailure ) {
       if ( dieOnError ) {
         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
         return -1;  // triggers early exit
       } else {
         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
       }
     } else if ( ! done ) {
       if ( dieOnError ) {
         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
         return -1;  // triggers early exit
       } else {
         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
         << (atom[ig].id + 1) << "!\n" << endi;
       }
     }
 
     // store data back to patch
     int ppoffset, partition;
     if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
       atom[ig+i].position = pos[i];
     } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
     } else for ( i = 0; i < hgs; ++i ) {
       Force df = netdp[i] * invdt;
       Tensor vir = outer(df, ref[i]);
       wc += vir;
       f[Results::normal][ig+i] += df;
       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
       if (ppreduction) {
         if (!i) {
           BigReal z = pos[i].z;
           int partition = atom[ig].partition;
           int slab = (int)floor((z-zmin)*idz);
           if (slab < 0) slab += nslabs;
           else if (slab >= nslabs) slab -= nslabs;
           ppoffset = 3*(slab + nslabs*partition);
         }
         ppreduction->item(ppoffset  ) += vir.xx;
         ppreduction->item(ppoffset+1) += vir.yy;
         ppreduction->item(ppoffset+2) += vir.zz;
       }
     }
   }
   if ( dt && virial ) *virial += wc;
 
   return 0;
 }

◆ rattle2()

void HomePatch::rattle2	(	const BigReal	timestep,
		Tensor *	virial
	)

Definition at line 4249 of file HomePatch.C.

References endi(), getWaterModelGroupSize(), iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Patch::numAtoms, Node::Object(), outer(), settle2(), Node::simParameters, simParams, SWM4, TIMEFACTOR, TIP4, Vector::x, Vector::y, and Vector::z.

 {
   Molecule *mol = Node::Object()->molecule;
   SimParameters *simParams = Node::Object()->simParameters;
   const int fixedAtomsOn = simParams->fixedAtomsOn;
   const int useSettle = simParams->useSettle;
   const BigReal dt = timestep / TIMEFACTOR;
   Tensor wc;  // constraint virial
   BigReal tol = simParams->rigidTol;
   int maxiter = simParams->rigidIter;
   int dieOnError = simParams->rigidDie;
   int i, iter;
   BigReal dsqi[10], tmp;
   int ial[10], ibl[10];
   Vector ref[10];  // reference position
   Vector refab[10];  // reference vector
   Vector vel[10];  // new velocity
   BigReal rmass[10];  // 1 / mass
   BigReal redmass[10];  // reduced mass
   int fixed[10];  // is atom fixed?
   
   // Size of a hydrogen group for water
   const WaterModel watmodel = simParams->watmodel;
   const int wathgsize = getWaterModelGroupSize(watmodel);
 
   //  CkPrintf("In rattle2!\n");
   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
     //    CkPrintf("ig=%d\n",ig);
     int hgs = atom[ig].hydrogenGroupSize;
     if ( hgs == 1 ) continue;  // only one atom in group
     // cache data in local arrays and integrate positions normally
     int anyfixed = 0;
     for ( i = 0; i < hgs; ++i ) {
       ref[i] = atom[ig+i].position;
       vel[i] = atom[ig+i].velocity;
       rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
     }
     int icnt = 0;
     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
       if (hgs != wathgsize) {
         NAMD_bug("Hydrogen group error caught in rattle2().");
       }
       // Use SETTLE for water unless some of the water atoms are fixed,
       if (useSettle && !anyfixed) {
         if (watmodel == WaterModel::SWM4) {
           // SWM4 ordering:  O D LP H1 H2
           // do swap(O,LP) and call settle with subarray O H1 H2
           // swap back after we return
           Vector lp_ref = ref[2];
           // Vector lp_vel = vel[2];
           ref[2] = ref[0];
           vel[2] = vel[0];
           settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
           ref[0] = ref[2];
           vel[0] = vel[2];
           ref[2] = lp_ref;
           // vel[2] = vel[0];  // set LP vel to O vel
         }
         else {
           settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
           if (watmodel == WaterModel::TIP4) {
             vel[3] = vel[0];
           }
         }
         for (i=0; i<hgs; i++) {
           atom[ig+i].velocity = vel[i];
         }
         continue;
       }
       if ( !(fixed[1] && fixed[2]) ) {
         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
       }
     }
     //    CkPrintf("Loop 2\n");
     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
         if ( !(fixed[0] && fixed[i]) ) {
           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
           ibl[icnt] = i;  ++icnt;
         }
       }
     }
     if ( icnt == 0 ) continue;  // no constraints
     //    CkPrintf("Loop 3\n");
     for ( i = 0; i < icnt; ++i ) {
       refab[i] = ref[ial[i]] - ref[ibl[i]];
     }
     //    CkPrintf("Loop 4\n");
     int done;
     for ( iter = 0; iter < maxiter; ++iter ) {
       done = 1;
       for ( i = 0; i < icnt; ++i ) {
         int a = ial[i];  int b = ibl[i];
         Vector vab = vel[a] - vel[b];
         Vector &rab = refab[i];
         BigReal rabsqi = dsqi[i];
         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
           wc += outer(dp,rab);
           vel[a] += rmass[a] * dp;
           vel[b] -= rmass[b] * dp;
           done = 0;
         }
       }
       if ( done ) break;
       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
     }
     if ( ! done ) {
       if ( dieOnError ) {
         NAMD_die("Exceeded maximum number of iterations in rattle2().");
       } else {
         iout << iWARN <<
           "Exceeded maximum number of iterations in rattle2().\n" << endi;
       }
     }
     // store data back to patch
     for ( i = 0; i < hgs; ++i ) {
       atom[ig+i].velocity = vel[i];
     }
   }
   //  CkPrintf("Leaving rattle2!\n");
   // check that there isn't a constant needed here!
   *virial += wc / ( 0.5 * dt );
 
 }

◆ receiveResult() [1/2]

void HomePatch::receiveResult ( ProxyGBISP1ResultMsg * msg )

Definition at line 5030 of file HomePatch.C.

References Sequencer::awaken(), Patch::boxesOpen, Flags::doNonbonded, Patch::flags, Patch::psiFin, ProxyGBISP1ResultMsg::psiSum, ProxyGBISP1ResultMsg::psiSumLen, and ResizeArray< Elem >::size().

Referenced by ProxyMgr::recvResult().

                                                        {
   ++numGBISP1Arrived;
     for ( int i = 0; i < msg->psiSumLen; ++i ) {
       psiFin[i] += msg->psiSum[i];
     }
   delete msg;
 
   if (flags.doNonbonded) {
     //awaken if phase 1 done
     if (phase1BoxClosedCalled == true &&
         numGBISP1Arrived==proxy.size() ) {
         // fprintf(stderr, "Calling awaken() on patch %d: 4\n", this->patchID);
         sequencer->awaken();
     }
   } else {
     //awaken if all phases done on noWork step
     if (boxesOpen == 0 &&
         numGBISP1Arrived == proxy.size() &&
         numGBISP2Arrived == proxy.size() &&
         numGBISP3Arrived == proxy.size()) {
       // fprintf(stderr, "Calling awaken() on patch %d: 5\n", this->patchID);
       sequencer->awaken();
     }
   }
 }

◆ receiveResult() [2/2]

void HomePatch::receiveResult ( ProxyGBISP2ResultMsg * msg )

Definition at line 5057 of file HomePatch.C.

References Sequencer::awaken(), Patch::boxesOpen, ProxyGBISP2ResultMsg::dEdaSum, ProxyGBISP2ResultMsg::dEdaSumLen, Patch::dHdrPrefix, Flags::doNonbonded, Patch::flags, and ResizeArray< Elem >::size().

                                                        {
   ++numGBISP2Arrived;
   //accumulate dEda
   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
     dHdrPrefix[i] += msg->dEdaSum[i];
   }
   delete msg;
 
   if (flags.doNonbonded) {
     //awaken if phase 2 done
     if (phase2BoxClosedCalled == true &&
         numGBISP2Arrived==proxy.size() ) {
       // fprintf(stderr, "Calling awaken() on patch %d: 6\n", this->patchID);
       sequencer->awaken();
     }
   } else {
     //awaken if all phases done on noWork step
     if (boxesOpen == 0 &&
         numGBISP1Arrived == proxy.size() &&
         numGBISP2Arrived == proxy.size() &&
         numGBISP3Arrived == proxy.size()) {
       // fprintf(stderr, "Calling awaken() on patch %d: 7\n", this->patchID);
       sequencer->awaken();
     }
   }
 }

◆ receiveResults() [1/3]

void HomePatch::receiveResults ( ProxyResultVarsizeMsg * msg )

Definition at line 837 of file HomePatch.C.

References OwnerBox< Owner, Data >::clientClose(), OwnerBox< Owner, Data >::clientOpen(), DebugM, Results::f, ProxyResultVarsizeMsg::flLen, ProxyResultVarsizeMsg::forceArr, Patch::forceBox, ProxyResultVarsizeMsg::isZero, Results::maxNumForces, ProxyResultVarsizeMsg::node, and Patch::patchID.

Referenced by ProxyMgr::recvResults().

                                                          {
 
     numGBISP3Arrived++;
     DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
     int n = msg->node;
     Results *r = forceBox.clientOpen();
 
     char *iszeroPtr = msg->isZero;
     Force *msgFPtr = msg->forceArr;
 
     for ( int k = 0; k < Results::maxNumForces; ++k )
     {
       Force *rfPtr = r->f[k];
       for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
           if((*iszeroPtr)!=1) {
               *rfPtr += *msgFPtr;
               msgFPtr++;
           }
       }      
     }
     forceBox.clientClose();
     delete msg;
 }

◆ receiveResults() [2/3]

void HomePatch::receiveResults ( ProxyResultMsg * msg )

Definition at line 861 of file HomePatch.C.

References ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::clientClose(), OwnerBox< Owner, Data >::clientOpen(), DebugM, ResizeArray< Elem >::end(), Results::f, Patch::f, Patch::forceBox, ProxyResultMsg::forceList, Results::maxNumForces, ProxyResultMsg::node, and Patch::patchID.

                                                   {
   numGBISP3Arrived++;
   DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
   int n = msg->node;
   Results *r = forceBox.clientOpen();
   for ( int k = 0; k < Results::maxNumForces; ++k )
   {
     Force *f = r->f[k];
     register ForceList::iterator f_i, f_e;
     f_i = msg->forceList[k]->begin();
     f_e = msg->forceList[k]->end();
     for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
   }
   forceBox.clientClose();
   delete msg;
 }

◆ receiveResults() [3/3]

void HomePatch::receiveResults ( ProxyCombinedResultRawMsg * msg )

Definition at line 878 of file HomePatch.C.

References OwnerBox< Owner, Data >::clientClose(), OwnerBox< Owner, Data >::clientOpen(), DebugM, Results::f, Patch::f, ProxyCombinedResultRawMsg::flLen, ProxyCombinedResultRawMsg::forceArr, Patch::forceBox, ProxyCombinedResultRawMsg::isForceNonZero, Results::maxNumForces, ProxyCombinedResultRawMsg::nodeSize, Patch::patchID, Vector::x, Vector::y, and Vector::z.

 {
     numGBISP3Arrived++;
   DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodeSize<<")\n");
     Results *r = forceBox.clientOpen(msg->nodeSize);
       register char* isNonZero = msg->isForceNonZero;
       register Force* f_i = msg->forceArr;
       for ( int k = 0; k < Results::maxNumForces; ++k )
       {
         Force *f = r->f[k];
         int nf = msg->flLen[k];
 #ifdef ARCH_POWERPC
 #pragma disjoint (*f_i, *f)
 #endif
         for (int count = 0; count < nf; count++) {
           if(*isNonZero){
             f[count].x += f_i->x;
             f[count].y += f_i->y;
             f[count].z += f_i->z;
             f_i++;
           }
           isNonZero++;
         }
       }
     forceBox.clientClose(msg->nodeSize);
 
     delete msg;
 }

◆ recvCheckpointAck()

void HomePatch::recvCheckpointAck ( )

Definition at line 5365 of file HomePatch.C.

Referenced by PatchMgr::recvCheckpointAck(), and recvCheckpointLoad().

                                   {  // initiating replica
   CkpvAccess(_qd)->process();
 }

◆ recvCheckpointLoad()

void HomePatch::recvCheckpointLoad ( CheckpointAtomsMsg * msg )

Definition at line 5301 of file HomePatch.C.

Referenced by PatchMgr::recvCheckpointLoad().

                                                           {  // initiating replica
   SimParameters *simParams = Node::Object()->simParameters;
   const int remote = simParams->scriptIntArg1;
   const char *key = simParams->scriptStringArg1;
   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
   }
   if ( msg->replica != remote ) {
     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
   }
   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
     strcpy(newmsg->key,key);
     newmsg->lattice = lattice;
     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
     newmsg->pid = patchID;
     newmsg->pe = CkMyPe();
     newmsg->replica = CmiMyPartition();
     newmsg->numAtoms = numAtoms;
     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
   }
   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
     lattice = msg->lattice;
     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
     numAtoms = msg->numAtoms;
     atom.resize(numAtoms);
     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
     doAtomUpdate = true;
     rattleListValid = false;
     rattleListValid_SOA = false;
     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
     if (simParams->SOAintegrateOn) {
       sort_solvent_atoms();
       copy_atoms_to_SOA();
 #if 0
       if (simParams->rigidBonds != RIGID_NONE) {
         buildRattleList_SOA();
         rattleListValid_SOA = true;
       }
 #endif
     }
   }
   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
     recvCheckpointAck();
   }
   delete msg;
 }

◆ recvCheckpointReq()

void HomePatch::recvCheckpointReq	(	int	task,
		const char *	key,
		int	replica,
		int	pe
	)

Definition at line 5271 of file HomePatch.C.

References CheckpointAtomsMsg::atoms, HomePatch::checkpoint_t::atoms, ResizeArray< Elem >::begin(), CheckpointAtomsMsg::berendsenPressure_count, HomePatch::checkpoint_t::berendsenPressure_count, checkpoints, CheckpointAtomsMsg::lattice, HomePatch::checkpoint_t::lattice, NAMD_die(), CheckpointAtomsMsg::numAtoms, HomePatch::checkpoint_t::numAtoms, PatchMgr::Object(), Patch::patchID, CheckpointAtomsMsg::pe, CheckpointAtomsMsg::pid, CheckpointAtomsMsg::replica, SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_SWAP, PatchMgr::sendCheckpointAck(), and PatchMgr::sendCheckpointLoad().

Referenced by PatchMgr::recvCheckpointReq().

                                                                                 {  // responding replica
   if ( task == SCRIPT_CHECKPOINT_FREE ) {
     if ( ! checkpoints.count(key) ) {
       NAMD_die("Unable to free checkpoint, requested key was never stored.");
     }
     delete checkpoints[key];
     checkpoints.erase(key);
     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
     return;
   }
   CheckpointAtomsMsg *msg;
   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
     if ( ! checkpoints.count(key) ) {
       NAMD_die("Unable to load checkpoint, requested key was never stored.");
     }
     checkpoint_t &cp = *checkpoints[key];
     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
     msg->lattice = cp.lattice;
     msg->berendsenPressure_count = cp.berendsenPressure_count;
     msg->numAtoms = cp.numAtoms;
     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
   } else {
     msg = new (0,1,0) CheckpointAtomsMsg;
   }
   msg->pid = patchID;
   msg->replica = CmiMyPartition();
   msg->pe = CkMyPe();
   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
 }

◆ recvCheckpointStore()

void HomePatch::recvCheckpointStore ( CheckpointAtomsMsg * msg )

Definition at line 5351 of file HomePatch.C.

References CheckpointAtomsMsg::atoms, HomePatch::checkpoint_t::atoms, ResizeArray< Elem >::begin(), CheckpointAtomsMsg::berendsenPressure_count, HomePatch::checkpoint_t::berendsenPressure_count, checkpoints, CheckpointAtomsMsg::key, CheckpointAtomsMsg::lattice, HomePatch::checkpoint_t::lattice, CheckpointAtomsMsg::numAtoms, HomePatch::checkpoint_t::numAtoms, PatchMgr::Object(), Patch::patchID, CheckpointAtomsMsg::pe, CheckpointAtomsMsg::replica, ResizeArray< Elem >::resize(), and PatchMgr::sendCheckpointAck().

Referenced by PatchMgr::recvCheckpointStore().

                                                            {  // responding replica
   if ( ! checkpoints.count(msg->key) ) {
     checkpoints[msg->key] = new checkpoint_t;
   }
   checkpoint_t &cp = *checkpoints[msg->key];
   cp.lattice = msg->lattice;
   cp.berendsenPressure_count = msg->berendsenPressure_count;
   cp.numAtoms = msg->numAtoms;
   cp.atoms.resize(cp.numAtoms);
   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
   delete msg;
 }

◆ recvExchangeMsg()

void HomePatch::recvExchangeMsg ( ExchangeAtomsMsg * msg )

Definition at line 5402 of file HomePatch.C.

References Patch::atomMapper, ExchangeAtomsMsg::atoms, ResizeArray< Elem >::begin(), buildRattleList_SOA(), ResizeArray< Elem >::end(), ExchangeAtomsMsg::lattice, Patch::lattice, ExchangeAtomsMsg::numAtoms, Patch::numAtoms, Node::Object(), rattleListValid, rattleListValid_SOA, AtomMapper::registerIDsFullAtom(), ResizeArray< Elem >::resize(), RIGID_NONE, Node::simParameters, simParams, and AtomMapper::unregisterIDsFullAtom().

Referenced by PatchMgr::recvExchangeMsg().

                                                      {
   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
   lattice = msg->lattice;
   numAtoms = msg->numAtoms;
   atom.resize(numAtoms);
   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
   delete msg;
   CkpvAccess(_qd)->process();
   doAtomUpdate = true;
   rattleListValid = false;
   rattleListValid_SOA = false;
   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
   SimParameters *simParams = Node::Object()->simParameters;
   if (simParams->SOAintegrateOn) {
     sort_solvent_atoms();
     copy_atoms_to_SOA();
 #if 0
     if (simParams->rigidBonds != RIGID_NONE) {
       buildRattleList_SOA();
       rattleListValid_SOA = true;
     }
 #endif
   }
 }

◆ recvExchangeReq()

void HomePatch::recvExchangeReq ( int req )

Definition at line 5391 of file HomePatch.C.

References exchange_dst, exchange_msg, exchange_req, PatchMgr::Object(), and PatchMgr::sendExchangeMsg().

Referenced by exchangeAtoms(), and PatchMgr::recvExchangeReq().

                                        {
   exchange_req = req;
   if ( exchange_msg ) {
     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
     exchange_msg = 0;
     exchange_req = -1;
     CkpvAccess(_qd)->process();
   }
 }

◆ recvNodeAwareSpanningTree()

void HomePatch::recvNodeAwareSpanningTree ( ProxyNodeAwareSpanningTreeMsg * msg )

Definition at line 671 of file HomePatch.C.

Referenced by ProxyMgr::recvNodeAwareSpanningTreeOnHomePatch().

671 {}

◆ recvSpanningTree()

void HomePatch::recvSpanningTree	(	int *	t,
		int	n
	)

Definition at line 677 of file HomePatch.C.

References proxySpanDim, ResizeArray< Elem >::resize(), sendSpanningTree(), and ResizeArray< Elem >::size().

Referenced by ProxyMgr::recvSpanningTreeOnHomePatch().

 {
   int i;
   nChild=0;
   tree.resize(n);
   for (i=0; i<n; i++) {
     tree[i] = t[i];
   }
 
   for (i=1; i<=proxySpanDim; i++) {
     if (tree.size() <= i) break;
     child[i-1] = tree[i];
     nChild++;
   }
 
 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
   isProxyChanged = 1;
 #endif
 
   // send down to children
   sendSpanningTree();
 }

◆ registerProxy()

void HomePatch::registerProxy ( RegisterProxyMsg * msg )

Definition at line 443 of file HomePatch.C.

References ResizeArray< Elem >::add(), ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::clientAdd(), DebugM, Patch::forceBox, RegisterProxyMsg::node, Patch::patchID, Random::reorder(), and ResizeArray< Elem >::size().

Referenced by ProxyMgr::recvRegisterProxy().

                                                    {
   DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
   proxy.add(msg->node);
   forceBox.clientAdd();
 
   isNewProxyAdded = 1;
 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
   isProxyChanged = 1;
 #endif
 
   Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
   delete msg;
 }

◆ replaceForces()

void HomePatch::replaceForces ( ExtForce * f )

Definition at line 2310 of file HomePatch.C.

References Patch::f.

 {
   replacementForces = f;
 }

◆ revert()

void HomePatch::revert ( void )

Definition at line 5232 of file HomePatch.C.

References Patch::atomMapper, ResizeArray< Elem >::begin(), buildRattleList_SOA(), ResizeArray< Elem >::copy(), ResizeArray< Elem >::end(), Patch::lattice, Patch::numAtoms, Node::Object(), rattleListValid, rattleListValid_SOA, AtomMapper::registerIDsFullAtom(), RIGID_NONE, Node::simParameters, simParams, ResizeArray< Elem >::size(), and AtomMapper::unregisterIDsFullAtom().

Referenced by Sequencer::algorithm().

                            {
   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
 
   atom.copy(checkpoint_atom);
   numAtoms = atom.size();
   lattice = checkpoint_lattice;
 
   doAtomUpdate = true;
   rattleListValid = false;
   rattleListValid_SOA = false;
 
   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
 
   SimParameters *simParams = Node::Object()->simParameters;
   if (simParams->SOAintegrateOn) {
     sort_solvent_atoms();
     copy_atoms_to_SOA();
 #if 0
     if (simParams->rigidBonds != RIGID_NONE) {
       buildRattleList_SOA();
       rattleListValid_SOA = true;
     }
 #endif
   }
 
   // DMK - Atom Separation (water vs. non-water)
   #if NAMD_SeparateWaters != 0
     numWaterAtoms = checkpoint_numWaterAtoms;
   #endif
 }

◆ runSequencer()

void HomePatch::runSequencer ( void )

Definition at line 305 of file HomePatch.C.

References Sequencer::run().

Referenced by Node::run().

306 { sequencer->run(); }

Sequencer::run

void run(void)

Definition: Sequencer.C:269

◆ saveForce()

void HomePatch::saveForce ( const int ftag = Results::normal )

Definition at line 2315 of file HomePatch.C.

References Patch::f, Patch::numAtoms, and ResizeArray< Elem >::resize().

Referenced by Sequencer::saveForce().

 {
   f_saved[ftag].resize(numAtoms);
   for ( int i = 0; i < numAtoms; ++i )
   {
     f_saved[ftag][i] = f[ftag][i];
   }
 }

◆ sendNodeAwareSpanningTree()

void HomePatch::sendNodeAwareSpanningTree ( )

Definition at line 672 of file HomePatch.C.

672 {}

◆ sendProxies()

void HomePatch::sendProxies ( )

Definition at line 509 of file HomePatch.C.

References ResizeArray< Elem >::begin(), ProxyMgr::Object(), Patch::patchID, ProxyMgr::sendProxies(), NodeProxyMgr::sendProxyList(), and ResizeArray< Elem >::size().

Referenced by ProxyMgr::buildProxySpanningTree2().

 {
 #if USE_NODEPATCHMGR
         CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
         NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
         npm->sendProxyList(patchID, proxy.begin(), proxy.size());
 #else
         ProxyMgr::Object()->sendProxies(patchID, proxy.begin(), proxy.size());
 #endif
 }

◆ sendSpanningTree()

void HomePatch::sendSpanningTree ( )

Definition at line 700 of file HomePatch.C.

References ResizeArray< Elem >::copy(), ProxySpanningTreeMsg::node, ProxyMgr::Object(), ProxySpanningTreeMsg::patch, Patch::patchID, ProxyMgr::sendSpanningTree(), ResizeArray< Elem >::size(), and ProxySpanningTreeMsg::tree.

Referenced by buildSpanningTree(), and recvSpanningTree().

 {
   if (tree.size() == 0) return;
   ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
   msg->patch = patchID;
   msg->node = CkMyPe();
   msg->tree.copy(tree);  // copy data for thread safety
   ProxyMgr::Object()->sendSpanningTree(msg);  
 }

◆ setGBISIntrinsicRadii()

void HomePatch::setGBISIntrinsicRadii ( )

Definition at line 4900 of file HomePatch.C.

References Patch::intRad, MassToRadius(), MassToScreen(), Node::molecule, Patch::numAtoms, Node::Object(), ResizeArray< Elem >::resize(), ResizeArray< Elem >::setall(), Node::simParameters, and simParams.

Referenced by positionsReady(), and positionsReady_SOA().

                                       {
   intRad.resize(numAtoms*2);
   intRad.setall(0);
   Molecule *mol = Node::Object()->molecule;
   SimParameters *simParams = Node::Object()->simParameters;
   Real offset = simParams->coulomb_radius_offset;
   for (int i = 0; i < numAtoms; i++) {
     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
     Real screen = MassToScreen(atom[i].mass);//same
     intRad[2*i+0] = rad - offset;//r0
     intRad[2*i+1] = screen*(rad - offset);//s0
   }
 }

◆ setLcpoType()

void HomePatch::setLcpoType ( )

Definition at line 4889 of file HomePatch.C.

References Molecule::getLcpoParamType(), Patch::lcpoType, Node::molecule, Patch::numAtoms, Node::Object(), Patch::pExt, and ResizeArray< Elem >::resize().

Referenced by positionsReady(), and positionsReady_SOA().

                             {
   Molecule *mol = Node::Object()->molecule;
   const int *lcpoParamType = mol->getLcpoParamType();
 
   lcpoType.resize(numAtoms);
   for (int i = 0; i < numAtoms; i++) {
     lcpoType[i] = lcpoParamType[pExt[i].id];
   }
 }

◆ submitLoadStats()

void HomePatch::submitLoadStats ( int timestep )

Definition at line 5428 of file HomePatch.C.

References Patch::numAtoms, LdbCoordinator::Object(), Patch::patchID, and LdbCoordinator::patchLoad().

Referenced by Sequencer::rebalanceLoad().

 {
   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
 }

◆ unregisterProxy()

void HomePatch::unregisterProxy ( UnregisterProxyMsg * msg )

Definition at line 457 of file HomePatch.C.

References ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::clientRemove(), ResizeArray< Elem >::del(), Patch::forceBox, and UnregisterProxyMsg::node.

Referenced by ProxyMgr::recvUnregisterProxy().

                                                        {
 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
   isProxyChanged = 1;
 #endif
   int n = msg->node;
   NodeID *pe = proxy.begin();
   for ( ; *pe != n; ++pe );
   forceBox.clientRemove();
   proxy.del(pe - proxy.begin());
   delete msg;
 }

◆ updateAtomBuffers()

void HomePatch::updateAtomBuffers ( )

◆ updateAtomCount()

void HomePatch::updateAtomCount	(	const int	n,
		const int	reallocate
	)

◆ useSequencer()

void HomePatch::useSequencer ( Sequencer * sequencerPtr )

Definition at line 301 of file HomePatch.C.

302 { sequencer=sequencerPtr; }

Friends And Related Function Documentation

◆ ComputeGlobal

friend class ComputeGlobal

friend

Definition at line 336 of file HomePatch.h.

◆ PatchMgr

friend class PatchMgr

friend

Definition at line 330 of file HomePatch.h.

◆ Sequencer

friend class Sequencer

friend

Definition at line 331 of file HomePatch.h.

◆ SequencerCUDA

friend class SequencerCUDA

friend

Definition at line 334 of file HomePatch.h.

Member Data Documentation

◆ checkpoint_task

int HomePatch::checkpoint_task

Definition at line 501 of file HomePatch.h.

Referenced by exchangeCheckpoint(), and recvCheckpointLoad().

◆ checkpoints

std::map<std::string,checkpoint_t*> HomePatch::checkpoints

Definition at line 508 of file HomePatch.h.

Referenced by recvCheckpointReq(), and recvCheckpointStore().

◆ exchange_dst

int HomePatch::exchange_dst

Definition at line 514 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

◆ exchange_msg

ExchangeAtomsMsg* HomePatch::exchange_msg

Definition at line 517 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

◆ exchange_req

int HomePatch::exchange_req

Definition at line 516 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

◆ exchange_src

int HomePatch::exchange_src

Definition at line 515 of file HomePatch.h.

Referenced by exchangeAtoms().

◆ gridForceIdxChecked

bool HomePatch::gridForceIdxChecked

Definition at line 402 of file HomePatch.h.

Referenced by ComputeGridForce::doForce(), and positionsReady().

◆ inMigration

int HomePatch::inMigration

protected

Definition at line 569 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

◆ ldObjHandle

LDObjHandle HomePatch::ldObjHandle

Definition at line 554 of file HomePatch.h.

Referenced by LdbCoordinator::initialize(), Sequencer::Sequencer(), Sequencer::suspend(), and Sequencer::terminate().

◆ marginViolations

int HomePatch::marginViolations

Definition at line 401 of file HomePatch.h.

Referenced by doAtomMigration(), doMarginCheck(), doMarginCheck_SOA(), Sequencer::multigratorPressure(), Sequencer::submitReductions(), and Sequencer::submitReductions_SOA().

◆ msgbuf

MigrateAtomsMsg* HomePatch::msgbuf[PatchMap::MaxOneAway]

protected

Definition at line 571 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

◆ noconstList

std::vector<int> HomePatch::noconstList

Definition at line 451 of file HomePatch.h.

Referenced by buildRattleList(), buildRattleList_SOA(), rattle1(), and rattle1_SOA().

◆ numMlBuf

int HomePatch::numMlBuf

protected

Definition at line 570 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

◆ posNew

std::vector<Vector> HomePatch::posNew

Definition at line 458 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

◆ rattleList

std::vector<RattleList> HomePatch::rattleList

Definition at line 449 of file HomePatch.h.

Referenced by buildRattleList(), buildRattleList_SOA(), rattle1(), and rattle1_SOA().

◆ rattleListValid

bool HomePatch::rattleListValid

Definition at line 453 of file HomePatch.h.

Referenced by positionsReady(), positionsReady_SOA(), rattle1(), recvCheckpointLoad(), recvExchangeMsg(), and revert().

◆ rattleListValid_SOA

bool HomePatch::rattleListValid_SOA

Definition at line 454 of file HomePatch.h.

Referenced by depositMigration(), positionsReady_SOA(), rattle1_SOA(), recvCheckpointLoad(), recvExchangeMsg(), and revert().

◆ rattleParam

std::vector<RattleParam> HomePatch::rattleParam

Definition at line 450 of file HomePatch.h.

Referenced by buildRattleList(), buildRattleList_SOA(), rattle1(), and rattle1_SOA().

◆ settleList

std::vector<int> HomePatch::settleList

Definition at line 448 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

◆ velNew

std::vector<Vector> HomePatch::velNew

Definition at line 457 of file HomePatch.h.

Referenced by addRattleForce(), buildRattleList(), and rattle1().

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

Classes

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

◆ ~HomePatch()

Member Function Documentation

◆ addForceToMomentum()

◆ addForceToMomentum3()

◆ addRattleForce()

◆ addVelocityToPosition()

◆ boxClosed()

◆ buildRattleList()

◆ buildRattleList_SOA()

◆ buildSpanningTree()

◆ checkpoint()

◆ depositMigration()

◆ doAtomMigration()

◆ doGroupSizeCheck()

◆ doGroupSizeCheck_SOA()

◆ doMarginCheck()

◆ doMarginCheck_SOA()

◆ doPairlistCheck()

◆ exchangeAtoms()

◆ exchangeCheckpoint()

◆ gbisComputeAfterP1()

◆ gbisComputeAfterP2()

◆ gbisP2Ready()

◆ gbisP3Ready()

◆ getAtomList()

◆ getMax()

◆ getMin()

◆ hardWallDrude()

◆ loweAndersenFinish()

◆ loweAndersenVelocities()

◆ minimize_rattle2()

◆ mollyAverage()

◆ mollyMollify()

◆ positionsReady()

◆ positionsReady_GPU()

◆ positionsReady_SOA()

◆ qmSwapAtoms()

◆ rattle1()

◆ rattle1_SOA()

◆ rattle1old()

◆ rattle2()

◆ receiveResult() [1/2]

◆ receiveResult() [2/2]

◆ receiveResults() [1/3]

◆ receiveResults() [2/3]

◆ receiveResults() [3/3]

◆ recvCheckpointAck()

◆ recvCheckpointLoad()

◆ recvCheckpointReq()

◆ recvCheckpointStore()

◆ recvExchangeMsg()

◆ recvExchangeReq()

◆ recvNodeAwareSpanningTree()

◆ recvSpanningTree()

◆ registerProxy()

◆ replaceForces()

◆ revert()

◆ runSequencer()

◆ saveForce()

◆ sendNodeAwareSpanningTree()

◆ sendProxies()

◆ sendSpanningTree()

◆ setGBISIntrinsicRadii()

◆ setLcpoType()

◆ submitLoadStats()

◆ unregisterProxy()

◆ updateAtomBuffers()

◆ updateAtomCount()

◆ useSequencer()

Friends And Related Function Documentation

◆ ComputeGlobal

◆ PatchMgr