Pair_style coul/wolf

LAMMPS LAMMPS Mailing List Mirror
1

Dear all,
I need some clarification in the c++ file of pair_style coul/wolf as I would like to modify this a bit and use a modified potential.

A part from the program is given below,

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

qisq = qtmpqtmp;
e_self = -(e_shift/2.0 + alf/MY_PIS) * qisqqqrd2e;
if (evflag) ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);

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;

if (rsq < cut_coulsq) {
r = sqrt(rsq);
prefactor = qqrd2eqtmpq[j]/r;
erfcc = erfc(alfr);
erfcd = exp(-alfalfrr);
v_sh = (erfcc - e_shiftr) * prefactor;
dvdrr = (erfcc/rsq + 2.0alf/MY_PIS * erfcd/r) + f_shift;
forcecoul = dvdrrrsqprefactor;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul / rsq;

f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}

if (eflag) {
ecoul = v_sh;
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;

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

My questions are

1)What is v_sh calculating? Is that the force term?

  1. And in the next step dvdrr, erf****cc is divided by rsq but the next step is divided by r.

3)What is the importance of factor_coul in this program?

Kindly clarify this. Thanks a lot in advance.

2

Dear all,
             I need some clarification in the c++ file of pair_style
coul/wolf as I would like to modify this a bit and use a modified
potential.

A part from the program is given below,

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

    qisq = qtmp*qtmp;
    e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e;
    if (evflag) ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);

    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;

      if (rsq < cut_coulsq) {
        r = sqrt(rsq);
        prefactor = qqrd2e*qtmp*q[j]/r;
        erfcc = erfc(alf*r);
        erfcd = exp(-alf*alf*r*r);
       * v_sh = (erfcc - e_shift*r) * prefactor;*
* dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift;*
* forcecoul = dvdrr*rsq*prefactor;*
* if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;*
* fpair = forcecoul / rsq;*

* f[i][0] += delx*fpair;*
* f[i][1] += dely*fpair;*
* f[i][2] += delz*fpair;*
* if (newton_pair || j < nlocal) {*
* f[j][0] -= delx*fpair;*
* f[j][1] -= dely*fpair;*
* f[j][2] -= delz*fpair;*
* }*

* if (eflag) {*
* ecoul = v_sh;*
* if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;*
* } else ecoul = 0.0;*

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

My questions are

1)What is *v_sh *calculating? Is that the force term?

​no. v_sh is the potential energy contribution. see the line "ecoul = v_sh"
​

2) And in the next step *dvdrr, **erf**cc* is divided by* rsq* but the
next step is divided by* r*.

i don't follow what you are after here, and there is no question. it seems
you are confusing erfcc and erfcd. ​

3)What is the importance of *factor_coul *in this program?

​it is a scaling factor for 1-2, 1-3, or 1-4 exclusions. it is important,
that this factor is applied to the unscreened coulomb.​
​for more info on exclusions, please see:
http://lammps.sandia.gov/doc/special_bonds.html

​axel.​

3

Thank you so much sir. I have some more doubts.

  1. My function is E= ((q1q2)/r) * erfc(alphar)

where r is the distance between the two atoms.

And this function is screened coulombic function. As it is screened, can I remove the unscreened coulomb part in that program? (that factor_coul portion)

  1. My second question is what dvrr is calculating? rsq-square of the distance between the atoms. r- distance between the atoms. My doubt is the first term is divided by rsq, while the second term is divided by r. Can you explain the motive behind this?

  2. Do we really need to include the self and shifted coulombic energy for the above mentioned equation also or that has been included for that wolf/coul ?

Thanking you.

4

Thank you so much sir. I have some more doubts.

1) My function is E= ((q1*q2)/r) * erfc(alpha*r)

where r is the distance between the two atoms.

And this function is screened coulombic function. As it is screened, can I
remove the unscreened coulomb part in that program? (that factor_coul
portion)

​no, but what you describe here as potential function is pair style
coul/long​, i.e. the real space part of an ewald sum or equivalent
long-range solver.

2) My second question is what dvrr is calculating? rsq-square of the
distance between the atoms. r- distance between the atoms. My doubt is the
first term is divided by rsq, while the second term is divided by r. Can
you explain the motive behind this?

3) Do we really need to include the self and shifted coulombic energy for
the above mentioned equation also or that has been included for that
wolf/coul ?

​sorry, i don't have the time to research this and do what is essentially
your work for you.​

in my personal opinion, nobody should not be modifying pair styles unless
that person understands the physics and math sufficiently well, and *then*
this person would be able to understand what is programmed. the only tricky
thing in the coulomb pair styles is the handling of excluded pairs, and it
actually took quite a while until somebody noticed and figured out how to
do it consistently, i.e. "the right way(tm)".

i don't consider it a good idea to do something first and then hope for
somebody to explain and/or correct it for you. that is not how research
works, IMNSHO.

​axel.​

5

Thank you for your reply sir.