LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Pair_style coul/wolf
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Pair_style coul/wolf

From: Kousika A <kousika12@...24...>
Date: Sat, 12 Aug 2017 20:18:00 +0530

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)

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 ?

Thanking you.

On Fri, Aug 11, 2017 at 5:42 PM, Axel Kohlmeyer <akohlmey@...24...> wrote:

On Fri, Aug 11, 2017 at 5:34 AM, Kousika A <kousika12@...24...> wrote:
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,

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, erfcc 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:



Kindly clarify this. Thanks a lot in advance.

Best regards,
Research Scholar
Department of Metallurgical and Materials Engineering

Check out the vibrant tech community on one of the world's most
engaging tech sites,!
lammps-users mailing list

Dr. Axel Kohlmeyer  akohlmey@...24...
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

Best regards,
Research Scholar
Department of Metallurgical and Materials Engineering