ComputeLjPmeSerialMgr Class Reference

Inheritance diagram for ComputeLjPmeSerialMgr:

Public Member Functions
	ComputeLjPmeSerialMgr ()
	~ComputeLjPmeSerialMgr ()
void	setCompute (ComputeLjPmeSerial *c)
void	recvCoord (LjPmeSerialCoordMsg *)
void	recvForce (LjPmeSerialForceMsg *)
void	getLJparameters (const int &i, const int &j, Real &A, Real &B, Real &A14, Real &B14)
void	get_vdw_parameters (Real &sigma, Real &epsilon, Real &sigma14, Real &epsilon14, Index index)
void	setLJparameters ()

Detailed Description

Definition at line 58 of file ComputeLjPmeSerial.C.

Constructor & Destructor Documentation

ComputeLjPmeSerialMgr()

ComputeLjPmeSerialMgr::ComputeLjPmeSerialMgr

(

)

Definition at line 100 of file ComputeLjPmeSerial.C.

                                             :
  ljPmeProxy(thisgroup), ljPmeCompute(0), numSources(0), numArrived(0),
  coordMsgs(0), coord(0), forceNonbond(0), forceSlow(0), oldmsg(0),
  numAtoms(0), ljPmeCoord(0), ljForceSlow(0), ljForceNonbond(0),
  ljParameter(0), ljParameter14(0), atomType(0)
{
  initialized = false;
  CkpvAccess(BOCclass_group).computeLjPmeSerialMgr = thisgroup;
}

~ComputeLjPmeSerialMgr()

ComputeLjPmeSerialMgr::~ComputeLjPmeSerialMgr

(

)

Definition at line 110 of file ComputeLjPmeSerial.C.

{
  initialized = false;
  for (int i=0;  i < numSources;  i++)  delete coordMsgs[i];
  delete [] coordMsgs;
  delete [] coord;
  delete [] forceNonbond;
  delete [] forceSlow;
  delete oldmsg;
  // extra
  if (ljPmeCoord) delete[] ljPmeCoord;
  if (ljForceSlow) delete[] ljForceSlow;
  if (ljForceNonbond) delete[] ljForceNonbond;
  if (ljParameter) delete[] ljParameter;
  if (ljParameter14) delete[] ljParameter14;
  if (atomType) delete[] atomType;
}

Member Function Documentation

get_vdw_parameters()

void ComputeLjPmeSerialMgr::get_vdw_parameters	(	Real &	sigma,
		Real &	epsilon,
		Real &	sigma14,
		Real &	epsilon14,
		Index	index
	)

Definition at line 236 of file ComputeLjPmeSerial.C.

References cbrt, and getLJparameters().

Referenced by setLJparameters().

                                      {
  // Get the index index interaction
  Real A, B, A14, B14;
  this->getLJparameters(index, index, A, B, A14, B14);
  // Copied from get_vdw_params in Paramters.h
  sigma = sqrt(cbrt(A)) / sqrt(cbrt(B));
  epsilon = B*B / (4*A);

  sigma14 = sqrt(cbrt(A14)) / sqrt(cbrt(B14));
  epsilon14 = B14*B14 / (4*A14);
}

getLJparameters()

void ComputeLjPmeSerialMgr::getLJparameters	(	const int &	i,
		const int &	j,
		Real &	A,
		Real &	B,
		Real &	A14,
		Real &	B14
	)

Returns epsilon and sigma term for nonbonded and 1-4 interaction Copied from Parameters.h

Definition at line 138 of file ComputeLjPmeSerial.C.

References Parameters::get_num_vdw_params(), Parameters::get_table_pair_params(), Parameters::get_vdw_pair_params(), Parameters::get_vdw_params(), Node::molecule, NAMD_bug(), NAMD_die(), Node::Object(), Node::parameters, Node::simParameters, simParams, Molecule::ss_num_vdw_params, and Molecule::ss_vdw_type.

Referenced by get_vdw_parameters(), and setLJparameters().

                                          {

  Parameters *params = Node::Object()->parameters;
  SimParameters *simParams = Node::Object()->simParameters;
  int useGeom = simParams->vdwGeometricSigma;
  Bool tabulatedEnergies = simParams->tabulatedEnergies;
  Bool soluteScalingOn = simParams->soluteScalingOn;
  BigReal soluteScalingFactor = simParams->soluteScalingFactorVdw;
  unsigned int table_dim_org = params->get_num_vdw_params();
  int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
  int K = -1;
  int *ss_vdw_type = Node::Object()->molecule->ss_vdw_type;
  // BigReal sigma_max;
  //  We need the A and B parameters for the Van der Waals.  These can
  //  be explicitly be specified for this pair or calculated from the
  //  sigma and epsilon values for the two atom types
  //  printf("Looking at interaction of  %i with %i\n", i, j);
  if ( tabulatedEnergies && params->get_table_pair_params(i,j,&K)) {

    if ( K < 0 ) { NAMD_bug(
        "LJPME::compute_vdw_params: energy table index is negative");
    }
    A = -1 - K;
    B = 0;
    A14 = -1 - K;
    B14 = 0;
  } else if (params->get_vdw_pair_params(i,j, &A, &B, &A14, &B14)) {
    if ( tabulatedEnergies && ( A < 0 || A14 < 0 ) ) {
      NAMD_die("LJ A is negative with tabulatedEnergies enabled");
    }
  } else {
    //  We didn't find explicit parameters for this pair. So instead,
    //  get the parameters for each atom type separately and use them
    //  to calculate the values we need
    Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
    Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;

    if (!soluteScalingOn) {
      params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
                &epsilon_i14,i);
      params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14, 
                &epsilon_j14,j);
    } else {
      int i_type = (i >= table_dim_org)? ss_vdw_type[i-table_dim_org]:i;
      int j_type = (j >= table_dim_org)? ss_vdw_type[j-table_dim_org]:j;
      params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
                                          &epsilon_i14,i_type);
      params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14,
                                          &epsilon_j14,j_type);
    }   
  
    BigReal sigma_ij =
        useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
    BigReal sigma_ij14 =
        useGeom ? sqrt(sigma_i14*sigma_j14) : 0.5 * (sigma_i14+sigma_j14);
    BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
    BigReal epsilon_ij14 = sqrt(epsilon_i14*epsilon_j14);

    //  Calculate sigma^6
    sigma_ij *= sigma_ij*sigma_ij;
    sigma_ij *= sigma_ij;
    sigma_ij14 *= sigma_ij14*sigma_ij14;
    sigma_ij14 *= sigma_ij14;
      
    //  Calculate LJ constants A & B
    B = 4.0 * sigma_ij * epsilon_ij;
    A = B * sigma_ij;
    B14 = 4.0 * sigma_ij14 * epsilon_ij14;
    A14 = B14 * sigma_ij14;

    if (soluteScalingOn) {
      if (i >= table_dim_org && i < (table_dim_org+ss_dim) && j < table_dim_org) {
        A *= sqrt(soluteScalingFactor);
        B *= sqrt(soluteScalingFactor);
        A14 *= sqrt(soluteScalingFactor);
        B14 *= sqrt(soluteScalingFactor);
      }
      if (i < table_dim_org && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
        A *= sqrt(soluteScalingFactor);
        B *= sqrt(soluteScalingFactor);
        A14 *= sqrt(soluteScalingFactor);
        B14 *= sqrt(soluteScalingFactor);
      }
      if (i >=table_dim_org && i < (table_dim_org+ss_dim) && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
        A *= soluteScalingFactor;
        B *= soluteScalingFactor;
        A14 *= soluteScalingFactor;
        B14 *= soluteScalingFactor;
      }
    }

    if ( tabulatedEnergies && ( A < 0 || A14 < 0 ) ) {
      NAMD_die("LJPME: LJ A is negative with tabulatedEnergies enabled");
    }
  }
}

recvCoord()

void ComputeLjPmeSerialMgr::recvCoord

(

LjPmeSerialCoordMsg *

msg

)

Definition at line 347 of file ComputeLjPmeSerial.C.

References LjPmeCompute::computeLJpotential(), LjPmeSerialCoordMsg::coord, LjPmeSerialForceMsg::energyNonbond, LjPmeSerialForceMsg::energySlow, LjPmeSerialForceMsg::forceNonbond, LjPmeSerialForceMsg::forceSlow, Parameters::get_num_vdw_params(), LjPmeSerialAtom::id, LjPmeCompute::initialize(), LjPmeSerialCoordMsg::lattice, Node::molecule, NAMD_die(), LjPmeSerialCoordMsg::numAtoms, Molecule::numAtoms, PatchMap::Object(), Node::Object(), Node::parameters, LjPmeSerialAtom::position, setLJparameters(), Node::simParameters, simParams, LjPmeSerialCoordMsg::sourceNode, LjPmeSerialForceMsg::virialNonbond, LjPmeSerialForceMsg::virialSlow, Vector::x, Vector::y, and Vector::z.

                                                              {
  SimParameters *simParams = Node::Object()->simParameters;
  Parameters *params = Node::Object()->parameters;
  if ( ! numSources ) {
    numSources = (PatchMap::Object())->numNodesWithPatches();
    coordMsgs = new LjPmeSerialCoordMsg*[numSources];
    for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
    numArrived = 0;
    numAtoms = Node::Object()->molecule->numAtoms;
    // Allocate memory to receive coords and send force
    coord = new LjPmeSerialAtom[numAtoms];
    forceNonbond = new LjPmeSerialForce[numAtoms];
    forceSlow = new LjPmeSerialForce[numAtoms];
  }

  // Receive the message at store using global atom index
  for (int i=0; i < msg->numAtoms; ++i) {
    coord[msg->coord[i].id] = msg->coord[i];
  }
  coordMsgs[numArrived] = msg;
  ++numArrived;

  if ( numArrived < numSources ) return;
  numArrived = 0;

  //###########################################
  // ALL DATA ARRIVED --- CALCULATE FORCES
  Lattice lattice = msg->lattice;

  if(!initialized) {
    initialized = true;
    // Allocate memory for LJ-PME
    table_dim = params->get_num_vdw_params();
    ljParameter = new double[2*table_dim*table_dim];
    ljParameter14 = new double[2*table_dim*table_dim];
    ljPmeCoord = new double[4*numAtoms];
    ljForceSlow = new double[3*numAtoms];
    ljForceNonbond = new double[3*numAtoms];
    atomType = new int[numAtoms];

    if (ljPmeCoord==0 || ljForceSlow==0 || ljForceNonbond==0 || ljParameter==0 
        || ljParameter==0 || atomType == 0) {
      NAMD_die("can't allocate LJ-PME atom buffers");
    }

    // Set the data for LJPME compute. To save memory, all arrays in LjPmeCompute
    // are pointing to arrayes in ComputeLjPmeSerialMgr class
    LjPme.initialize(ljPmeCoord, ljParameter, ljParameter14, atomType,
                    ljForceNonbond, ljForceSlow, table_dim, numAtoms);
    // Set LJ parameters
    this->setLJparameters();
    oldNonbondedScaling = simParams->nonbondedScaling;
    oldSoluteScaling = simParams->soluteScalingFactorVdw;
  }

  // Nonbonded scaling or solute VDW scaling has changed, so update the LJ parameters
  if((simParams->nonbondedScaling != oldNonbondedScaling) ||
    (simParams->soluteScalingOn && (simParams->soluteScalingFactorVdw != oldSoluteScaling))) {
    this->setLJparameters();
    oldNonbondedScaling = simParams->nonbondedScaling;
    oldSoluteScaling = simParams->soluteScalingFactorVdw;
  }

  // Update the coordinates and initialize the forces to zero
  for (int i = 0;  i < numAtoms;  i++) {
    ljPmeCoord[4*i  ] = coord[i].position.x;
    ljPmeCoord[4*i+1] = coord[i].position.y;
    ljPmeCoord[4*i+2] = coord[i].position.z;
    // Initialize the force to zero
    ljForceSlow[3*i  ] = 0.0;
    ljForceSlow[3*i+1] = 0.0;
    ljForceSlow[3*i+2] = 0.0;
    ljForceNonbond[3*i  ] = 0.0;
    ljForceNonbond[3*i+1] = 0.0;
    ljForceNonbond[3*i+2] = 0.0;
  }

  double energySlow, energyNonbond;
  double virialSlow[3][3], virialNonbond[3][3];
  LjPme.computeLJpotential(simParams->LJPMEEwaldCoefficient,
        simParams->cutoff,
        lattice,
        virialNonbond,
        virialSlow,
        energyNonbond,
        energySlow,
        true,
        true);
  
  if (0) {
    printf("LJ-PME: E_nb: %6.12f, E_rec: %6.12f, E_sum: %6.12f! \n", 
    energyNonbond, energySlow, energyNonbond+energySlow);
    printf("LJ-PME: virial_nb: [%3.6f, %3.6f, %3.6f, %3.6f, %3.6f, %3.6f] \n",
        virialNonbond[0][0], virialNonbond[0][1], virialNonbond[0][2],
        virialNonbond[1][1], virialNonbond[1][2], virialNonbond[2][2]);
    printf("LJ-PME: virial_slow: [%3.6f, %3.6f, %3.6f, %3.6f, %3.6f, %3.6f] \n",
        virialSlow[0][0], virialSlow[0][1], virialSlow[0][2],
        virialSlow[1][1], virialSlow[1][2], virialSlow[2][2]);
  }

  for (int i = 0;  i < numAtoms;  i++) {
    forceSlow[i].x = ljForceSlow[3*i  ];
    forceSlow[i].y = ljForceSlow[3*i+1];
    forceSlow[i].z = ljForceSlow[3*i+2];

    forceNonbond[i].x = ljForceNonbond[3*i  ];
    forceNonbond[i].y = ljForceNonbond[3*i+1];
    forceNonbond[i].z = ljForceNonbond[3*i+2];

    if (0) {
      forceNonbond[i].x = 0.0;
      forceNonbond[i].y = 0.0;
      forceNonbond[i].z = 0.0;
      double fx = forceNonbond[i].x + forceSlow[i].x;
      double fy = forceNonbond[i].y + forceSlow[i].y;
      double fz = forceNonbond[i].z + forceSlow[i].z;
      printf("LJ-PME send: atom %5d<-> nb: (%3.6f,%3.6f,%3.6f), slow:(%3.6f,%3.6f,%3.6f), sum:(%3.6f,%3.6f,%3.6f)! \n", i,
      forceNonbond[i].x, forceNonbond[i].y, forceNonbond[i].z,
      forceSlow[i].x, forceSlow[i].y, forceSlow[i].z,
      fx, fy, fz);
    }
  }

  // distribute forces
  for (int j=0;  j < numSources;  j++) {
    LjPmeSerialCoordMsg *cmsg = coordMsgs[j];
    coordMsgs[j] = 0;
    LjPmeSerialForceMsg *fmsg = new (cmsg->numAtoms, cmsg->numAtoms, 0) LjPmeSerialForceMsg;

    for (int i=0;  i < cmsg->numAtoms;  i++) {
      fmsg->forceSlow[i] = forceSlow[cmsg->coord[i].id];
      fmsg->forceNonbond[i] = forceNonbond[cmsg->coord[i].id];
    }

    if ( ! j ) {  // set virial and energy only for first message
      fmsg->energyNonbond = energyNonbond;
      fmsg->energySlow = energySlow;
      for (int k=0;  k < 3;  k++) {
        for (int l=0;  l < 3;  l++) {
          fmsg->virialSlow[k][l] = virialSlow[k][l];
          fmsg->virialNonbond[k][l] = virialNonbond[k][l];
        }
      }
    }
    else {  // set other messages to zero, add into reduction only once
      fmsg->energyNonbond = 0;
      fmsg->energySlow = 0;
      for (int k=0;  k < 3;  k++) {
        for (int l=0;  l < 3;  l++) {
          fmsg->virialSlow[k][l] = 0;
          fmsg->virialNonbond[k][l] = 0;
        }
      }
    }

    ljPmeProxy[cmsg->sourceNode].recvForce(fmsg);
    delete cmsg;
  }
}

recvForce()

void ComputeLjPmeSerialMgr::recvForce

(

LjPmeSerialForceMsg *

msg

)

Returns Attraction and Repulsion term for nonbonded and 1-4 interaction Copied from LJTable.C

Definition at line 508 of file ComputeLjPmeSerial.C.

References ComputeLjPmeSerial::saveResults().

                                                              {
  ljPmeCompute->saveResults(msg);
  delete oldmsg;
  oldmsg = msg;
}

setCompute()

void ComputeLjPmeSerialMgr::setCompute ( ComputeLjPmeSerial *  c )

inline

Definition at line 63 of file ComputeLjPmeSerial.C.

{ ljPmeCompute = c; }

setLJparameters()

void ComputeLjPmeSerialMgr::setLJparameters

(

)

Definition at line 249 of file ComputeLjPmeSerial.C.

References get_vdw_parameters(), getLJparameters(), Node::Object(), Node::simParameters, simParams, and LjPmeSerialAtom::vdwType.

Referenced by recvCoord().

                                            {
  SimParameters *simParams = Node::Object()->simParameters;
  Real scaling = simParams->nonbondedScaling;
  // These won't change
  for (int i = 0; i < table_dim; i++) {
    for (int j = i; j < table_dim; j++) {
      int index_ij = 2*(i*table_dim+j);
      int index_ji = 2*(j*table_dim+i);
      Real A, B, A14, B14;
      this->getLJparameters(i, j, A, B, A14, B14);
      ljParameter[index_ij] = scaling * A; // repulsion
      ljParameter[index_ij+1] = scaling * B; // attraction
      ljParameter14[index_ij] = scaling * A14; // 1-4 repulsion
      ljParameter14[index_ij+1] = scaling * B14; // 1-4 attraction
      // Copy to transpose entry
      ljParameter[index_ji] = ljParameter[index_ij];
      ljParameter[index_ji+1] = ljParameter[index_ij+1];
      ljParameter14[index_ji] = ljParameter14[index_ij];
      ljParameter14[index_ji+1] = ljParameter14[index_ij+1];
  
      //printf("Table Dim: %d, i: %d, j: %d \n", table_dim, i, j);
      //printf("A: %3.6f, B: %3.6f, A14: %3.6f, B14: %3.6f \n", A, B, A14, B14);
    }
  }

  // Store the dispersion part of LJ using Geometric sigma combining rules
  for (int i = 0;  i < numAtoms;  i++) {
    Real sigma, sigma_14, epsilon, epsilon_14;
    atomType[i] = coord[i].vdwType;
    this->get_vdw_parameters(sigma, epsilon, sigma_14,
        epsilon_14, coord[i].vdwType);
    // printf("VDWtype: %d, eps: %3.6f, sigma: %3.6f, eps14: %3.6f, sigma14: %3.6f \n",
    //         atomType[i], epsilon, sigma, epsilon_14, sigma_14);

    // store sqrt(4*eps*sig^6) = 2*sqrt(eps)*sigma^3
    ljPmeCoord[4*i+3] = 2.0*sigma*sigma*sigma*sqrt(scaling * epsilon);
  }
}

---

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

• ComputeLjPmeSerial.C