Hello lammps users,
lammps version 11 august 2017.
I have two questions
-
pair style lj/cut/coul/wolf doest not work. I get the following error unknown pair style. Can someone please verify this.
-
I am trying to verify the energies form pair style lj/cut/coul/dsf with the calculation that I have done offline.
Consider the following input script.
Two LJ particles with charge +1 and -1 at a distance of 1A.
units real
atom_style full
boundary f f f
region 1 block 0 3 0 3 0 3
create_box 1 1
create_atoms 1 single 1 0 0
create_atoms 1 single 2 0 0
mass 1 1.0
set atom 11 charge 1
set atom 22 charge -1
pair_style lj/cut/coul/dsf 0.2 2
pair_coeff 1 1 0 0
thermo 1
thermo_style custom pe evdwl ecoul elong
run 0
write_data data.test
The output I get from lammps is -412.4329 kcal/mole.
PotEng E_vdwl E_coul E_long
-412.4329 0 -412.4329 0
However when I calculate it offline using the expression in the doc pages I get -83.8224 kcal/mole
To further verify my offline calculation I tried with pair style lj/cut/cou/cut 2 2 I get
PotEng E_vdwl E_coul E_long
-332.06371 0 -332.06371 0
and my calculations give Ecoul = -332.06 which matches with lammps.
I would appreciate if anybody could clarify how pair style lj/cut/coul/dsf calculations are done or point out if there is any obvious mistake that I am making.
The following octave script might be useful if anyone is trying to verify.
clc;clear all;
q= [-1 1] 18.2223;
alpha = 0.2;
r=1;
rc = 2;
Ecoul = (q(1)q(2))/r
t1=erfc(alpharc)/(rcrc);
t2=(2alphaexp(-alphaalpharcrc))/(sqrt(pi)rc);
Esf = q(1)q(2)(erfc(alphar)/r - erfc(alpharc)/rc + (r-rc)*(t1 + t2))
save as dsf.m
and run as
$octave dsf.m
Thank you
Arun