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 = delx*delx + dely*dely + 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_coul*forcecoul * r2inv;

f[i][0] += delx*fpair;
f[i][1] += dely*fpair;

f[i][2] += delz

*fpair;*

if (newton_pair || j < nlocal) {

f[j][0] -= delxfpair;

if (newton_pair || j < nlocal) {

f[j][0] -= delx

f[j][1] -= dely

*fpair;*

f[j][2] -= delzfpair;

f[j][2] -= delz

}

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