more specific question regarding new potentials

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

Ok, so I am looking at this function here:

void PairCoulCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double rsq,r2inv,rinv,forcecoul,factor_coul;
int *ilist,*jlist,*numneigh,**firstneigh;

ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;

double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;

inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;

// loop over neighbors of my atoms

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;

delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
jtype = type[j];

if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
forcecoul = qqrd2e * scale[itype][jtype] * qtmp*q[j]rinv;
fpair = factor_coulforcecoul * r2inv;

f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delzfpair;
}

if (eflag)
ecoul = factor_coul * qqrd2e * scale[itype][jtype] * qtmp*q[j]*rinv;

if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz);
}
}
}

if (vflag_fdotr) virial_fdotr_compute();
}

There appears to be a distance-dielectric function in pair_cut_diel.cpp and this one is from pair_coul_cut.cpp. I want to add new functions for non-polar solvation and possibly electrostatics in the future… I was looking at this reference implementation which may need to be modified:

http://www.ncbi.nlm.nih.gov/pubmed/17918282

My specific question is that if I only update the energy and not the forces, would this work with a hybrid MD/MC that uses MD to propose moves and MC to accept? I was going to start by just computing the energies and doing MC and then put in the full forces later. From a theoretical perspective the answer is yes (since the MD moves fulfill detailed balance), but I am asking from the perspective of LAMMPS. So I would just not update the forces yet and use zeros, but add to either the vdw or coulomb energies and use the GCMC in LAMMPS along with fix_nvt.

Thanks,
Luke

2

Paul can answer, but I think for a pair style to be
usable with fix gcmc it needs both energy and forces.
Also, I don’t think fix gcmc is the same as what you
may be referring to as hybrid MD/MC where you use
MD to equilibrate before (possibly with a big timestep)
before making the MC accept/reject decision. At
one point we had an initial implementation of a new
command like that, but never released it in LAMMPS.

Steve

3

Thanks Steve,

Yes, I learned a bit more when I set up the GCMC calculation last night. The NVT with Verlet is just used to do regular MD and occasionally only insertion and deletion of particles/molecules is attempted. So the only MC that is implemented is just the translation and rotation of molecules, not a hybrid MD… so I have to write it myself.

So what I would like to do is create a new fix that reassigns velocities from a normal distribution at the target temperature at every step, then does a short verlet integration (it could specified as options), and then accepts or rejects based on the energy of the states.

That would be the first step, and then my second step would be writing a new pair style based on one of the ones like that example I pasted there, because this part is quite understandable to me. Writing the new fix will be more difficult but I looked at the MC/fix_gcmc.cpp and it is fairly comprehensible.

So I would appreciate any suggestions you guys have as to the best way to code this would be. If there is code that never made it into LAMMPS that’s available I would like to take a look at it, and I am just copying the templates from existing source folders (.h and .cpp).

My goal here is to first implement a new energy function for non-polar solvation free energy that would depend on coordinates of the system and how buried individual atoms are within an aggregate, and then possibly do the electrostatics next.

So this non-polar solvation is a computational geometry problem but I wrote all the C++ code for this method in JCTC: http://pubs.acs.org/doi/abs/10.1021/ct060025%2B and that method relies heavily on grid-based algorithms for computing relative polymer chain positions and overlap. So I think it is feasible for me to attempt this and my main questions would relate to how I can get my code into LAMMPS.

I would start with just the single CPU version and then put in the MPI and OpenMP support. So perhaps Paul could have some suggestions on how to attempt this implementing code for new fixes and potentials into LAMMPS, because I can handle the computational geometry issues and getting my algorithms working.

Thanks,
Luke

Paul can answer, but I think for a pair style to be
usable with fix gcmc it needs both energy and forces.
Also, I don’t think fix gcmc is the same as what you
may be referring to as hybrid MD/MC where you use
MD to equilibrate before (possibly with a big timestep)
before making the MC accept/reject decision. At
one point we had an initial implementation of a new
command like that, but never released it in LAMMPS.

Steve

4

I will send you the code that we tried out for hybrid MD in a separate

email.

Other comments below.

Steve

5

I do have a question here from looking through the source and documentation. Why does this HMC accept/reject based on kinetic+potential energy?

It would seem like it would work just for accepting on potential energy because we are sampling the configuration partition function Z (and not Q) here and the short MD to propose the move just has to fulfill detailed balance, i.e., reversibility. I know a little bit of the theory here so this sounds right to me, what do you think? Perhaps this was written with a different goal in mind?

Thanks,
Luke

6

Perhaps that’s a bug in this old documentation though, since it says it must have the “fix nve”, so accepting based on kinetic+potential energy has to be a misprint since that should be nearly constant… I will work on it :slight_smile:

7

I no longer remember - too long ago - I believe

that is what Ed Maginn (the expert) thought was
appropriate.

Steve

8

Hi Steve,

I got it working with the current build but there’s a few bugs. I need a “run” command (even “run 0” seems to work!) before the hmc command or it complains that “Energy was not tallied on needed timestep” right away, from compute_pe.cpp line 75… Even though I had “thermo 1” and using the default computes of thermo_temp and thermo_pe.

Also, If the move is too big between accepting and rejecting, for some reason the temperature, PE, KE, and almost everything else besides electrostatic energy explodes to crazy numbers like 1402848593396 and the thing crashes. Something must be missing in the code here but I have no idea what it is.

But smaller steps seem fine and I was able to do 1,000,000 iterations with a small timestep of 2fs (for a coarse-grained model) and accept/rejecting every 10 integrations. The problem seems to be that it’s hard to propose a move with MD that doesn’t increase the PE, because all the velocities are random, so when there are bonds in the system you tend to be stretching them out.

So unfortunately it seems like a terribly inefficient method but perhaps it could be useful for some specific application. It can accept/reject based on PE or KE+PE. Ill send you the code later, I just moved so we didn’t get the Internet set up yet.

Luke

9

Strictly speaking, fix gcmc only needs energies. But practically speaking, you have to define a pair style, and I’ve never tried doing that without defining forces. And if you want to do regular MD moves, you’d certainly need forces.

Paul

10

Ø So perhaps Paul could have some suggestions on how to attempt this implementing code for new fixes and potentials into LAMMPS, because I can handle the computational geometry issues and getting my algorithms working.

We do have a process and suggestions for writing new LAMMPS functions and getting new contributed code into LAMMPS. See:

http://lammps.sandia.gov/contribute.html

http://lammps.sandia.gov/workshops/Aug13/Plimpton/lammps_modify_Aug13.pdf

Paul