HomePatch Class Reference

#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	saveForce (const int ftag=Results::normal)
void	addForceToMomentum (FullAtom *__restrict atom_arr, const Force *__restrict force_arr, const BigReal dt, int num_atoms)
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)
void	addVelocityToPosition (FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms)
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	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 ()
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)
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 ()
int	getNumFixedAtoms ()
void	setNumFixedAtoms (int numFixed)
PatchID	getPatchID ()
int	getNumComputes ()
CompAtomExt *	getCompAtomExtInfo ()
CudaAtom *	getCudaAtomList ()

Public Attributes
int	marginViolations
std::vector< int >	settleList
std::vector< RattleList >	rattleList
std::vector< RattleParam >	rattleParam
std::vector< int >	noconstList
bool	rattleListValid
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	doMarginCheck ()
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	ComputeGlobal

Detailed Description

Definition at line 273 of file HomePatch.h.

Constructor & Destructor Documentation

HomePatch::~HomePatch

(

)

Definition at line 321 of file HomePatch.C.

References Patch::atomMapper, ResizeArray< T >::begin(), ResizeArray< T >::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

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

Definition at line 1933 of file HomePatch.C.

References SimParameters::fixedAtomsOn, num_atoms, Node::Object(), and Node::simParameters.

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;
    }
  }
}

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 1962 of file HomePatch.C.

References SimParameters::fixedAtomsOn, num_atoms, Node::Object(), and Node::simParameters.

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;
    }
  }
}

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

Definition at line 2372 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];
  }
}

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

Definition at line 2001 of file HomePatch.C.

References SimParameters::fixedAtomsOn, num_atoms, Node::Object(), and Node::simParameters.

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;
    }
  }
}

void HomePatch::boxClosed ( int  box )

protectedvirtual

Implements Patch.

Definition at line 334 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, Patch::patchID, Patch::psiSumBox, and ResizeArray< T >::size().

                                 {
  // begin gbis
  if (box == 5) {// end of phase 1
    phase1BoxClosedCalled = true;
    if (!psiSumBox.isOpen() && numGBISP1Arrived == proxy.size()) {
      if (flags.doGBIS && flags.doNonbonded) {
        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) {
        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
      sequencer->awaken();
    }
  } else {
    DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
  }
}

void HomePatch::buildRattleList

(

)

Definition at line 2230 of file HomePatch.C.

References RattleParam::dsq, SimParameters::fixedAtomsOn, RattleParam::ia, RattleParam::ib, HomePatch::RattleList::icnt, HomePatch::RattleList::ig, NAMD_die(), noconstList, Patch::numAtoms, Node::Object(), posNew, rattleList, rattleParam, RattleParam::rma, RattleParam::rmb, settle1init(), settleList, Node::simParameters, SimParameters::useSettle, velNew, WAT_SWM4, WAT_TIP4, and SimParameters::watmodel.

Referenced by rattle1().

                                {

  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
  int wathgsize = 3;
  int watmodel = simParams->watmodel;
  if (watmodel == WAT_TIP4) wathgsize = 4;
  else if (watmodel == WAT_SWM4) wathgsize = 5;

  // 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 (simParams->watmodel == WAT_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_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);
    }
  }

}

void HomePatch::buildSpanningTree

(

void

)

Definition at line 674 of file HomePatch.C.

References ResizeArray< T >::begin(), compDistance(), ResizeArray< T >::end(), ResizeArray< T >::find(), NAMD_bug(), PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), PatchMap::Object(), Patch::patchID, proxySpanDim, ResizeArray< T >::resize(), sendSpanningTree(), ResizeArray< T >::setall(), ResizeArray< T >::size(), and ResizeArray< T >::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();
}

void HomePatch::checkpoint

(

void

)

Definition at line 3450 of file HomePatch.C.

References ResizeArray< T >::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
}

void HomePatch::depositMigration

(

MigrateAtomsMsg *

msg

)

Definition at line 3999 of file HomePatch.C.

References ResizeArray< T >::add(), Lattice::apply_transform(), Sequencer::awaken(), ResizeArray< T >::begin(), DebugM, ResizeArray< T >::end(), inMigration, Patch::lattice, MigrateAtomsMsg::migrationList, msgbuf, Lattice::nearest(), numMlBuf, Patch::patchID, Lattice::reverse_transform(), and ResizeArray< T >::size().

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) {
    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;
      sequencer->awaken();
      return;
    }
    else {
       DebugM(3,"patch "<<patchID<<" was not suspended\n");
    }
  }
}

void HomePatch::doAtomMigration ( )

protected

Definition at line 3839 of file HomePatch.C.

References ResizeArray< T >::add(), atomIndex, Patch::atomMapper, ResizeArray< T >::begin(), DebugM, ResizeArray< T >::del(), depositMigration(), ResizeArray< T >::end(), inMigration, Patch::lattice, marginViolations, MigrationInfo::mList, msgbuf, NAMD_EVENT_START, NAMD_EVENT_STOP, numMlBuf, PatchMgr::Object(), Patch::patchID, ResizeArray< T >::resize(), Lattice::scale(), PatchMgr::sendMigrationMsgs(), ResizeArray< T >::size(), Sequencer::suspend(), AtomMapper::unregisterIDsFullAtom(), Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

{
  int i;
#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
  NAMD_EVENT_START(1, NamdProfileEvent::ATOM_MIGRATIONS);
#endif
  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();

  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;

#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
  NAMD_EVENT_STOP(1, NamdProfileEvent::ATOM_MIGRATIONS);
#endif

  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;
}

void HomePatch::doGroupSizeCheck ( )

protected

Definition at line 3718 of file HomePatch.C.

References ResizeArray< T >::begin(), Flags::doNonbonded, ResizeArray< T >::end(), Patch::flags, SimParameters::hgroupCutoff, Flags::maxGroupRadius, NAMD_bug(), Node::Object(), Node::simParameters, x, y, and 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);

}

void HomePatch::doMarginCheck ( )

protected

Definition at line 3762 of file HomePatch.C.

References Lattice::a(), Lattice::a_r(), Lattice::b(), Lattice::b_r(), ResizeArray< T >::begin(), Lattice::c(), Lattice::c_r(), SimParameters::cutoff, ResizeArray< T >::end(), Patch::lattice, SimParameters::margin, marginViolations, NAMD_die(), Node::Object(), SimParameters::patchDimension, Lattice::scale(), Node::simParameters, 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;
  // } 

}

void HomePatch::doPairlistCheck ( )

protected

Definition at line 3629 of file HomePatch.C.

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

Referenced by positionsReady().

{
#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

  SimParameters *simParams = Node::Object()->simParameters;

  if ( numAtoms == 0 || ! flags.usePairlists ) {
    flags.pairlistTolerance = 0.;
    flags.maxAtomMovement = 99999.;
#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
    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 defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
    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 ) {
    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 ) ) {
    doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
  }

  if ( max_tol > doPairlistCheck_newTolerance ) {
    doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
  }
#if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
    NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
#endif
}

void HomePatch::exchangeAtoms

(

int

scriptTask

)

Definition at line 3574 of file HomePatch.C.

References ExchangeAtomsMsg::atoms, ResizeArray< T >::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, SimParameters::scriptArg1, SimParameters::scriptArg2, PatchMgr::sendExchangeReq(), and Node::simParameters.

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);
  }
}

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

Definition at line 3478 of file HomePatch.C.

References checkpoint_task, Node::Object(), PatchMgr::Object(), Patch::patchID, SimParameters::scriptIntArg1, SimParameters::scriptStringArg1, PatchMgr::sendCheckpointReq(), and Node::simParameters.

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);
}

void HomePatch::gbisComputeAfterP1

(

)

Definition at line 3147 of file HomePatch.C.

References SimParameters::alpha_max, Patch::bornRad, SimParameters::coulomb_radius_offset, Patch::flags, SimParameters::gbis_beta, SimParameters::gbis_delta, SimParameters::gbis_gamma, gbisP2Ready(), Patch::intRad, Patch::numAtoms, Node::Object(), Patch::pExt, Patch::psiFin, Patch::psiSum, Flags::sequence, and Node::simParameters.

                                   {

  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();
}

void HomePatch::gbisComputeAfterP2

(

)

Definition at line 3175 of file HomePatch.C.

References Patch::bornRad, COULOMB, SimParameters::coulomb_radius_offset, Patch::dEdaSum, Patch::dHdrPrefix, SimParameters::dielectric, Patch::flags, SimParameters::gbis_beta, SimParameters::gbis_delta, SimParameters::gbis_gamma, gbisP3Ready(), Patch::intRad, SimParameters::kappa, Patch::numAtoms, Node::Object(), Patch::pExt, Patch::psiFin, Flags::sequence, Node::simParameters, and SimParameters::solvent_dielectric.

                                   {

  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();
}

void HomePatch::gbisP2Ready

(

)

Definition at line 3219 of file HomePatch.C.

References ResizeArray< T >::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< T >::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();
}

void HomePatch::gbisP3Ready

(

)

Definition at line 3239 of file HomePatch.C.

References ResizeArray< T >::begin(), ProxyGBISP3DataMsg::destPe, ProxyGBISP3DataMsg::dHdrPrefix, Patch::dHdrPrefix, ProxyGBISP3DataMsg::dHdrPrefixLen, Flags::doFullElectrostatics, Patch::flags, GB3_PROXY_DATA_PRIORITY, Patch::gbisP3Ready(), ProxyGBISP3DataMsg::origPe, ProxyGBISP3DataMsg::patch, PATCH_PRIORITY, Patch::patchID, PRIORITY_SIZE, Flags::sequence, SET_PRIORITY, and ResizeArray< T >::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();
}

FullAtomList& HomePatch::getAtomList ( )

inline

Definition at line 455 of file HomePatch.h.

Referenced by ComputeHomePatch::doWork(), and dumpbench().

{ return (atom); }

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

Definition at line 2024 of file HomePatch.C.

References BOLTZMANN, Lattice::c(), SimParameters::drudeBondLen, SimParameters::drudeHardWallOn, SimParameters::drudeTemp, endi(), SimParameters::fixedAtomsOn, iERROR(), iout, SubmitReduction::item(), Patch::lattice, Node::molecule, Patch::numAtoms, Node::Object(), Lattice::origin(), outer(), partition(), SimParameters::pressureProfileSlabs, SimParameters::rigidDie, Node::simParameters, TIMEFACTOR, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, 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;
}

void HomePatch::loweAndersenFinish

(

)

Definition at line 3113 of file HomePatch.C.

References DebugM.

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

void HomePatch::loweAndersenVelocities

(

)

Definition at line 3098 of file HomePatch.C.

References DebugM, Node::molecule, Patch::numAtoms, Node::Object(), and Node::simParameters.

Referenced by positionsReady().

{
    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");
}

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

Definition at line 2963 of file HomePatch.C.

References ResizeArray< T >::begin(), endi(), Patch::f, SimParameters::fixedAtomsOn, iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Results::normal, Patch::numAtoms, Node::Object(), outer(), SimParameters::rigidDie, SimParameters::rigidIter, SimParameters::rigidTol, settle2(), Node::simParameters, TIMEFACTOR, SimParameters::useSettle, WAT_SWM4, WAT_TIP4, SimParameters::watmodel, 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
  int wathgsize = 3;
  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;

  //  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 (simParams->watmodel == WAT_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 (simParams->watmodel == WAT_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 );

}

void HomePatch::mollyAverage

(

)

Definition at line 3313 of file HomePatch.C.

References average(), endi(), SimParameters::fixedAtomsOn, iout, iWARN(), Node::molecule, SimParameters::mollyIter, SimParameters::mollyTol, NAMD_die(), Patch::numAtoms, Node::Object(), Patch::p_avg, ResizeArray< T >::resize(), and Node::simParameters.

Referenced by positionsReady().

{
  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;
  //    }
  // }

}

void HomePatch::mollyMollify

(

Tensor *

virial

)

Definition at line 3387 of file HomePatch.C.

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

{
  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);
}

void HomePatch::positionsReady

(

int

doMigration = 0

)

Definition at line 925 of file HomePatch.C.

References ProxyDataMsg::avgPlLen, ProxyDataMsg::avgPositionList, Patch::avgPositionPtrBegin, Patch::avgPositionPtrEnd, ResizeArray< T >::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< T >::end(), SimParameters::fixedAtomsForces, SimParameters::fixedAtomsOn, ProxyDataMsg::flags, Patch::flags, Patch::getPatchID(), Patch::intRad, ProxyDataMsg::intRadList, Patch::lattice, Patch::lcpoType, ProxyDataMsg::lcpoTypeList, loweAndersenVelocities(), mollyAverage(), NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NamdProfileEventStr, Patch::numAtoms, PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), PatchMap::Object(), Sync::Object(), Node::Object(), ProxyMgr::Object(), 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, CudaAtom::q, SimParameters::qmLSSOn, qmSwapAtoms(), rattleListValid, ResizeArray< T >::resize(), ComputeNonbondedUtil::scaling, ProxyMgr::sendProxyAll(), ProxyMgr::sendProxyData(), Flags::sequence, SET_PRIORITY, setGBISIntrinsicRadii(), setLcpoType(), Node::simParameters, ResizeArray< T >::size(), sortAtomsForCUDA(), CompAtomExt::sortOrder, SimParameters::staticAtomAssignment, Lattice::unscale(), SimParameters::useAVXTiles, Patch::v, ProxyDataMsg::velocityList, Patch::velocityPtrBegin, Patch::velocityPtrEnd, ProxyDataMsg::vlLen, Vector::x, CudaAtom::x, Vector::y, CudaAtom::y, Vector::z, and CudaAtom::z.

{    
  SimParameters *simParams = Node::Object()->simParameters;

  flags.sequence++;

  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_HIP) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
  #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
    cudaAtomList.resize(n);
    CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
    #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
    memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
    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_HIP) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
    #ifdef NAMD_AVXTILES
    if (simParams->useAVXTiles)
    #endif
    memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
#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);

}

void HomePatch::qmSwapAtoms

(

)

Definition at line 866 of file HomePatch.C.

References ResizeArray< T >::begin(), CompAtom::charge, SortedArray< Type >::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(), SimParameters::PMEOn, Node::simParameters, and CompAtom::vdwType.

Referenced by positionsReady().

{
    // 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;
}

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

Definition at line 2382 of file HomePatch.C.

References addRattleForce(), buildRattleList(), endi(), SimParameters::fixedAtomsOn, iERROR(), iout, iWARN(), noconstList, Patch::numAtoms, Node::Object(), posNew, rattle1old(), rattleList, rattleListValid, rattleN(), rattlePair< 1 >(), rattleParam, SimParameters::rigidDie, SimParameters::rigidIter, SimParameters::rigidTol, settle1_SIMD< 1 >(), settle1_SIMD< 2 >(), settleList, Node::simParameters, TIMEFACTOR, SimParameters::useSettle, velNew, WAT_TIP3, SimParameters::watmodel, 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 != WAT_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 {
      rattleN(icnt, &rattleParam[posParam],
        refx, refy, refz,
        posx, posy, posz,
        tol2, maxiter,
        done, consFailure);
    }

    // Advance position in rattleParam
    posParam += icnt;

    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;
}

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

Definition at line 2559 of file HomePatch.C.

References Lattice::c(), endi(), Patch::f, SimParameters::fixedAtomsOn, iERROR(), iout, SubmitReduction::item(), iWARN(), Patch::lattice, Node::molecule, NAMD_die(), Results::normal, Patch::numAtoms, Node::Object(), Lattice::origin(), outer(), partition(), SimParameters::pressureProfileSlabs, SimParameters::rigidDie, SimParameters::rigidIter, SimParameters::rigidTol, settle1(), settle1init(), Node::simParameters, TIMEFACTOR, SimParameters::useSettle, WAT_SWM4, WAT_TIP4, SimParameters::watmodel, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, z, and Tensor::zz.

Referenced by rattle1().

{
  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;

  // 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 (simParams->watmodel == WAT_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_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;
  }

  // Size of a hydrogen group for water
  int wathgsize = 3;
  int watmodel = simParams->watmodel;
  if (watmodel == WAT_TIP4) wathgsize = 4;
  else if (watmodel == WAT_SWM4) wathgsize = 5;
  
  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 (simParams->watmodel == WAT_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 (simParams->watmodel == WAT_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;
}

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

Definition at line 2829 of file HomePatch.C.

References endi(), SimParameters::fixedAtomsOn, iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Patch::numAtoms, Node::Object(), outer(), SimParameters::rigidDie, SimParameters::rigidIter, SimParameters::rigidTol, settle2(), Node::simParameters, TIMEFACTOR, SimParameters::useSettle, WAT_SWM4, WAT_TIP4, SimParameters::watmodel, 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
  int wathgsize = 3;
  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;

  //  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 (simParams->watmodel == WAT_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 (simParams->watmodel == WAT_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 );

}

void HomePatch::receiveResult

(

ProxyGBISP1ResultMsg *

msg

)

Definition at line 3262 of file HomePatch.C.

References Sequencer::awaken(), Patch::boxesOpen, Flags::doNonbonded, Patch::flags, Patch::psiFin, ProxyGBISP1ResultMsg::psiSum, ProxyGBISP1ResultMsg::psiSumLen, and ResizeArray< T >::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() ) {
        sequencer->awaken();
    }
  } else {
    //awaken if all phases done on noWork step
    if (boxesOpen == 0 &&
        numGBISP1Arrived == proxy.size() &&
        numGBISP2Arrived == proxy.size() &&
        numGBISP3Arrived == proxy.size()) {
      sequencer->awaken();
    }
  }
}

void HomePatch::receiveResult

(

ProxyGBISP2ResultMsg *

msg

)

Definition at line 3287 of file HomePatch.C.

References Sequencer::awaken(), Patch::boxesOpen, ProxyGBISP2ResultMsg::dEdaSum, ProxyGBISP2ResultMsg::dEdaSumLen, Patch::dHdrPrefix, Flags::doNonbonded, Patch::flags, and ResizeArray< T >::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() ) {
      sequencer->awaken();
    }
  } else {
    //awaken if all phases done on noWork step
    if (boxesOpen == 0 &&
        numGBISP1Arrived == proxy.size() &&
        numGBISP2Arrived == proxy.size() &&
        numGBISP3Arrived == proxy.size()) {
      sequencer->awaken();
    }
  }
}

void HomePatch::receiveResults

(

ProxyResultVarsizeMsg *

msg

)

Definition at line 796 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;
}

void HomePatch::receiveResults

(

ProxyResultMsg *

msg

)

Definition at line 820 of file HomePatch.C.

References ResizeArray< T >::begin(), OwnerBox< Owner, Data >::clientClose(), OwnerBox< Owner, Data >::clientOpen(), DebugM, ResizeArray< T >::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;
}

void HomePatch::receiveResults

(

ProxyCombinedResultRawMsg *

msg

)

Definition at line 837 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;
}

void HomePatch::recvCheckpointAck

(

)

Definition at line 3569 of file HomePatch.C.

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

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

void HomePatch::recvCheckpointLoad

(

CheckpointAtomsMsg *

msg

)

Definition at line 3516 of file HomePatch.C.

References Patch::atomMapper, CheckpointAtomsMsg::atoms, ResizeArray< T >::begin(), CheckpointAtomsMsg::berendsenPressure_count, Sequencer::berendsenPressure_count, checkpoint_task, ResizeArray< T >::end(), CheckpointAtomsMsg::key, CheckpointAtomsMsg::lattice, Patch::lattice, NAMD_bug(), CheckpointAtomsMsg::numAtoms, Patch::numAtoms, Node::Object(), PatchMgr::Object(), Patch::patchID, CheckpointAtomsMsg::pe, CheckpointAtomsMsg::pid, rattleListValid, recvCheckpointAck(), AtomMapper::registerIDsFullAtom(), CheckpointAtomsMsg::replica, ResizeArray< T >::resize(), SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_STORE, SCRIPT_CHECKPOINT_SWAP, SimParameters::scriptIntArg1, SimParameters::scriptStringArg1, PatchMgr::sendCheckpointStore(), Node::simParameters, and AtomMapper::unregisterIDsFullAtom().

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;
    if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
  }
  if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
    recvCheckpointAck();
  }
  delete msg;
}

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

Definition at line 3486 of file HomePatch.C.

References CheckpointAtomsMsg::atoms, HomePatch::checkpoint_t::atoms, ResizeArray< T >::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);
}

void HomePatch::recvCheckpointStore

(

CheckpointAtomsMsg *

msg

)

Definition at line 3555 of file HomePatch.C.

References CheckpointAtomsMsg::atoms, HomePatch::checkpoint_t::atoms, ResizeArray< T >::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< T >::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;
}

void HomePatch::recvExchangeMsg

(

ExchangeAtomsMsg *

msg

)

Definition at line 3606 of file HomePatch.C.

References Patch::atomMapper, ExchangeAtomsMsg::atoms, ResizeArray< T >::begin(), ResizeArray< T >::end(), ExchangeAtomsMsg::lattice, Patch::lattice, ExchangeAtomsMsg::numAtoms, rattleListValid, AtomMapper::registerIDsFullAtom(), ResizeArray< T >::resize(), 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;
  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
}

void HomePatch::recvExchangeReq

(

int

req

)

Definition at line 3595 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();
  }
}

void HomePatch::recvNodeAwareSpanningTree

(

ProxyNodeAwareSpanningTreeMsg *

msg

)

Definition at line 630 of file HomePatch.C.

Referenced by ProxyMgr::recvNodeAwareSpanningTreeOnHomePatch().

{}

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

Definition at line 636 of file HomePatch.C.

References proxySpanDim, ResizeArray< T >::resize(), sendSpanningTree(), and ResizeArray< T >::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();
}

void HomePatch::registerProxy

(

RegisterProxyMsg *

msg

)

Definition at line 402 of file HomePatch.C.

References ResizeArray< T >::add(), ResizeArray< T >::begin(), OwnerBox< Owner, Data >::clientAdd(), DebugM, Patch::forceBox, RegisterProxyMsg::node, Patch::patchID, Random::reorder(), and ResizeArray< T >::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;
}

void HomePatch::replaceForces

(

ExtForce *

f

)

Definition at line 1331 of file HomePatch.C.

References Patch::f.

{
  replacementForces = f;
}

void HomePatch::revert

(

void

)

Definition at line 3460 of file HomePatch.C.

References Patch::atomMapper, ResizeArray< T >::begin(), ResizeArray< T >::copy(), ResizeArray< T >::end(), Patch::lattice, rattleListValid, AtomMapper::registerIDsFullAtom(), ResizeArray< T >::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;

  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());

  // DMK - Atom Separation (water vs. non-water)
  #if NAMD_SeparateWaters != 0
    numWaterAtoms = checkpoint_numWaterAtoms;
  #endif
}

void HomePatch::runSequencer

(

void

)

Definition at line 269 of file HomePatch.C.

References Sequencer::run().

Referenced by Node::run().

{ sequencer->run(); }

void HomePatch::saveForce

(

const int

ftag = Results::normal

)

Definition at line 1336 of file HomePatch.C.

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

Referenced by Sequencer::saveForce().

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

void HomePatch::sendNodeAwareSpanningTree

(

)

Definition at line 631 of file HomePatch.C.

{}

void HomePatch::sendProxies

(

)

Definition at line 468 of file HomePatch.C.

References ResizeArray< T >::begin(), ProxyMgr::Object(), Patch::patchID, ProxyMgr::sendProxies(), NodeProxyMgr::sendProxyList(), and ResizeArray< T >::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
}

void HomePatch::sendSpanningTree

(

)

Definition at line 659 of file HomePatch.C.

References ResizeArray< T >::copy(), ProxySpanningTreeMsg::node, ProxyMgr::Object(), ProxySpanningTreeMsg::patch, Patch::patchID, ProxyMgr::sendSpanningTree(), ResizeArray< T >::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);  
}

void HomePatch::setGBISIntrinsicRadii

(

)

Definition at line 3132 of file HomePatch.C.

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

Referenced by positionsReady().

                                      {
  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
  }
}

void HomePatch::setLcpoType

(

)

Definition at line 3121 of file HomePatch.C.

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

Referenced by positionsReady().

                            {
  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];
  }
}

void HomePatch::submitLoadStats

(

int

timestep

)

Definition at line 3620 of file HomePatch.C.

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

Referenced by Sequencer::rebalanceLoad().

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

void HomePatch::unregisterProxy

(

UnregisterProxyMsg *

msg

)

Definition at line 416 of file HomePatch.C.

References ResizeArray< T >::begin(), OwnerBox< Owner, Data >::clientRemove(), ResizeArray< T >::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;
}

void HomePatch::useSequencer

(

Sequencer *

sequencerPtr

)

Definition at line 265 of file HomePatch.C.

{ sequencer=sequencerPtr; }

Friends And Related Function Documentation

friend class ComputeGlobal

friend

Definition at line 276 of file HomePatch.h.

friend class PatchMgr

friend

Definition at line 274 of file HomePatch.h.

friend class Sequencer

friend

Definition at line 275 of file HomePatch.h.

Member Data Documentation

int HomePatch::checkpoint_task

Definition at line 428 of file HomePatch.h.

Referenced by exchangeCheckpoint(), and recvCheckpointLoad().

<