[lammps-users] Issue With Electrostaic Force Calculation

Dear LAMMPS Users,

I am simulating a system composed of 125 ionic liquid ion-pairs (the cation has atom types 1 to 3 and the anion atom type is 4) and single water molecule (atom type 5 for O and 6 for H) . The aim is to compute the electrostatic energy/force on the O atom due to all IL atoms from a pre-equlibrated trajectory by invoking the rerun command. To make the test scenario simple, the pair style is set to

pair_style coul/cut 16.0.

The group ID for IL atoms (based on types 1 2 3 4) is IP and water oxygen atom is OW (type 5).

group IP type 1 2 3 4

group OW type 5

The pairwise electrostatic interaction energy on OW are computed in two ways:

compute 1 OW group/group IP
compute 3 OW pe/atom pair
compute 4 OW reduce sum c_3

with the output command

thermo_style custom c_1 c_4
thermo 1

I am not understanding why c_1 and c_4 values are different in the log file. Also, electrostaic force on OW is computed and dumped by the following commands:

compute 2 OW property/atom type fx fy fz

dump 1 OW custom 1 allForce.txt c_2[*]

dump_modify 1 sort id

The force components computed in this way is different from the force components from compute ID 1 (c_1[*]).

Please correct me on the protocols I am using. The script is attached as a reference.

Sourav

Penn State

inp.txt (1 KB)

Dear LAMMPS Users,

I am simulating a system composed of 125 ionic liquid ion-pairs (the cation has atom types 1 to 3 and the anion atom type is 4) and single water molecule (atom type 5 for O and 6 for H) . The aim is to compute the electrostatic energy/force on the O atom due to all IL atoms from a pre-equlibrated trajectory by invoking the rerun command. To make the test scenario simple, the pair style is set to

pair_style coul/cut 16.0.

The group ID for IL atoms (based on types 1 2 3 4) is IP and water oxygen atom is OW (type 5).

group IP type 1 2 3 4

group OW type 5

The pairwise electrostatic interaction energy on OW are computed in two ways:

compute 1 OW group/group IP
compute 3 OW pe/atom pair
compute 4 OW reduce sum c_3

with the output command

thermo_style custom c_1 c_4
thermo 1

I am not understanding why c_1 and c_4 values are different in the log file. Also, electrostaic force on OW is computed and dumped by the following commands:

these two computations are not computing the same thing. with pe/atom all interacting pairs are computed and then for each pair half of the energy is added to each of the two participating atoms. while for compute group/group the sum is for the whole interaction for each pair. this difference would become more significant, if you had more atoms included that are neither in group OW or group IP. it is difficult to make more specific comments without knowing more about the specifics of your system.

compute 2 OW property/atom type fx fy fz

dump 1 OW custom 1 allForce.txt c_2[*]

dump_modify 1 sort id

The force components computed in this way is different from the force components from compute ID 1 (c_1[*]).

this is the same issue.

axel.

Dear Axel,

Thanks for the initial clarifications. I am providing a bit more details of the simulated system for further reference:

The IL is a 4-site model where the cation is represented by rigid 3 (type 1, 2 and 3) beads and anion is a single bead (type 4).
Total no of IL ion-pairs is 125. A single water molecule (type 5 for O and type 6 for H) is embeded in the IL box. The water is treated with rigid SPC/E model . The trajectory was generated with pair_style lj/cut/coul/long 16.0. The aim is to decompose the total electrostatic force/energy arising due to all IL atoms on the water O into cationic and anionic contributions
(Fx=Fx_cat + Fx_anion, Fy=Fy_cat + Fy_anion and Fz=Fz_cat + Fz_anion and E=E_cat+E_anion). The initial test rerun is performed on simpler pair_style coul/cut 16.0 to check the consistency of the decompositions. So, I created groups for cation (type 1, 2 and 3) , anion (type 4) and water O (type 5) as follows:

group IP type 1 2 3 4 # grouping of IL atom types

group cat type 1 2 3 # grouping of IL cation atom types

group anion type 4 # grouping of IL anion atom types

group OW type 5 #grouping of H2O O

Then the following computes are invoked:

compute 1 OW group/group IP # pairwise interaction among all IL atoms with O
compute 2 OW property/atom type fx fy fz # force arising due to pairwise interactions on O
compute 3 OW pe/atom pair #PE arising due to pairwise interactions on O

compute 4 OW reduce sum c_3 # sum to dump with thermo command (will not have any effect due single O atom in the system)
compute 5 OW group/group cat # pairwise interaction among all IL cation atoms with O
compute 6 OW group/group anion # pairwise interaction among all IL anion atoms with O

I expect that

c_1=c_5+c_6 (PE, which I found consistent)
c_1[]=c_5[]+c_6[] (Force components, which I found consistent)
c_1=c_3 (Not the case)
c_1[
]=c_2[] (Not the case)
c_2[
]=c_5[]+c_6[] (Not the case).

So, I am confused which is the protocol I should take to solve the problem.

Thanks

Sourav

Penn State

Well, I have given you the necessary information to tell the two approaches apart. With that you should read the documentation again and then make an informed decision. That is something you have to do yourself. You know best what exactly you want.
Axel.