LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Pair_style coul/wolf

# Re: [lammps-users] Pair_style coul/wolf

 From: Kousika A 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 wrote:

On Fri, Aug 11, 2017 at 5:34 AM, Kousika A 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];

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, 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.​

​axel.​

Kindly clarify this. Thanks a lot in advance.

--
Best regards,
A.KOUSIKA
Research Scholar
Department of Metallurgical and Materials Engineering

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
_______________________________________________
lammps-users mailing list
lammps-users@...655....net
https://lists.sourceforge.net/lists/listinfo/lammps-users

--
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

--
Best regards,
A.KOUSIKA
Research Scholar
Department of Metallurgical and Materials Engineering