Discrepency in comparison of relative well depth for side-to-side and end-to-end in the gayberne potential

Dear LAMMPS users,

I am using the 8 Feb 2023 version of LAMMPS and comparing the side-to-side well depth with the end-to-end well depth using Gayberne potential.

pair_coeff command allows us to control the side-to-side (\epsilon_a), face-to-face (\epsilon_b), and end-to-end (\epsilon_c) well depths. To check these things, I started writing the coordinate of 2 ellipsoidal particles corresponding to side-to-side and end-to-end orientations. I evaluated the potential energy by giving only 1 run.

I have used the [\epsilon_a \epsilon_b \epsilon_c as 1.0 1.0 0.2]. The side-to-side well depth is expected to be 5 times deeper compared to the end-to-end well depth. But this is what I am not observing; rather, they give the same well depth. By changing \epsilon_c, there was no effect on the side-to-side and end-to-end well depth.

I have tried different combinations as [\epsilon_a \epsilon_b \epsilon_c as 0.2 1.0 1.0]. It gives the side-to-side well depth 5 times deeper compared to end-to-end well depth. But this is not according to the guideline provided along with the Gayberne potential.

I really stuck with this discrepancy, and please help me out with this issue.

The Gay-Berne potential also depends on the shape of the ellipsoids, their contact distance \sigma_c, and the exponents \upsilon,\mu. Therefore, playing with the relative well depths alone is not enough.

I have prepared a small sandbox to play with the parameters, plot the potential energy, and to visualise the position of the particles. Here you go:

Script
# Input variables to create the energy curves.
variable cutoff    equal 5.   # Long cut-off radius == smaller truncation error.
variable efflength equal 4.   # length of the potential curve
variable steps     equal 200  # number of points of the potential curve.
variable box       equal 5.   # box size
variable go        equal 0.4  # Begin of the plotting curve.
variable move_velocity equal "(v_efflength/v_steps)"

# Create the particles.
units lj
boundary f f f
atom_style ellipsoid
atom_modify map yes # needed for accessing a value from an atom vector!
region my_box block 0. ${box} 0. ${box} 0. ${box}
create_box 1 my_box
create_atoms 1 single 0 0 0
create_atoms 1 single 0 0 0
group gb type 1

# Force field.
neighbor        1 bin
pair_style  gayberne 1. 1. 1. ${cutoff}
# oblate
#pair_coeff * * 1. .5 0.2 0.2 1 0.2 0.2 1
#set group gb shape 1 1 .5
# prolate
pair_coeff * * 1. 1. 1.0 1.0 0.2 1.0 1.0 0.2
set group gb shape .5 .5 1
set group gb quat 0.0 0.0 0.0 1.0
group base_grp id 1
group move_grp id 2

# Measure and control the position of the moving particle.
compute pos gb property/atom x y z
variable dist equal "sqrt((c_pos[1][1]-c_pos[2][1])^2 + (c_pos[1][2]-c_pos[2][2])^2 + (c_pos[1][3]-c_pos[2][3])^2)"
variable rx equal "-1*c_pos[2][1]"
variable ry equal "-1*c_pos[2][2]"
variable rz equal "-1*c_pos[2][3]"

# Output.
timestep 1.
thermo_style custom v_dist evdwl
thermo_modify line one format 1 data%11.6f
thermo_modify norm no
thermo_modify flush yes
thermo 1

compute q all property/atom quatw quati quatj quatk
compute diameter all property/atom shapex shapey shapez
dump 1 all custom $(v_steps/10) gb_pot.dump id type x y z c_q[1] c_q[2] c_q[3] c_q[4] c_diameter[1] c_diameter[2] c_diameter[3]
dump_modify 1 colname c_q[1] quatw colname c_q[2] quati colname c_q[3] quatj colname c_q[4] quatk

# Scan along +X
displace_atoms move_grp move ${go} 0 0
fix 1 move_grp move linear ${move_velocity} 0.0 0.0
echo log
log             gb_pot_x.dat
run ${steps}

# Scan along +Y
unfix 1
displace_atoms move_grp move v_rx v_ry v_rz
displace_atoms move_grp move 0 ${go} 0
fix 1 move_grp move linear 0.0 ${move_velocity} 0.0
log             gb_pot_y.dat
run ${steps}

# Scan along +Z
unfix 1
displace_atoms move_grp move v_rx v_ry v_rz
displace_atoms move_grp move 0 0 ${go}
fix 1 move_grp move linear 0.0 0.0 ${move_velocity}
log             gb_pot_z.dat
run ${steps}

And these are the results. Visualise with:
gnuplot> p [:][:2] 'gb_pot_x.dat' u 2:3 w l,'gb_pot_z.dat' u 2:3 w l
Oblate
image image

Prolate
image image

2 Likes

Thanks for your help. Now I am getting the expected results from your script.

I am attaching my own script here. I have taken care of all the parameters, but I didn’t get the expected results and was unable to understand where I was making a mistake. Please help me to figure out the fault.

variable b loop 61
variable r index 1
variable f index sep=1.00.dat sep=1.05.dat sep=1.10.dat sep=1.15.dat sep=1.20.dat sep=1.25.dat sep=1.30.dat sep=1.35.dat sep=1.40.dat sep=1.45.dat sep=1.50.dat sep=1.55.dat sep=1.60.dat sep=1.65.dat sep=1.70.dat sep=1.75.dat sep=1.80.dat sep=1.85.dat sep=1.90.dat sep=1.95.dat sep=2.00.dat sep=2.05.dat sep=2.10.dat sep=2.15.dat sep=2.20.dat sep=2.25.dat sep=2.30.dat sep=2.35.dat sep=2.40.dat sep=2.45.dat sep=2.50.dat sep=2.55.dat sep=2.60.dat sep=2.65.dat sep=2.70.dat sep=2.75.dat sep=2.80.dat sep=2.85.dat sep=2.90.dat sep=2.95.dat sep=3.00.dat sep=3.05.dat sep=3.10.dat sep=3.15.dat sep=3.20.dat sep=3.25.dat sep=3.30.dat sep=3.35.dat sep=3.40.dat sep=3.45.dat sep=3.50.dat sep=3.55.dat sep=3.60.dat sep=3.65.dat sep=3.70.dat sep=3.75.dat sep=3.80.dat sep=3.85.dat sep=3.90.dat sep=3.95.dat sep=4.00.dat

variable g index sep=1.00 sep=1.05 sep=1.10 sep=1.15 sep=1.20 sep=1.25 sep=1.30 sep=1.35 sep=1.40 sep=1.45 sep=1.50 sep=1.55 sep=1.60 sep=1.65 sep=1.70 sep=1.75 sep=1.80 sep=1.85 sep=1.90 sep=1.95 sep=2.00 sep=2.05 sep=2.10 sep=2.15 sep=2.20 sep=2.25 sep=2.30 sep=2.35 sep=2.40 sep=2.45 sep=2.50 sep=2.55 sep=2.60 sep=2.65 sep=2.70 sep=2.75 sep=2.80 sep=2.85 sep=2.90 sep=2.95 sep=3.00 sep=3.05 sep=3.10 sep=3.15 sep=3.20 sep=3.25 sep=3.30 sep=3.35 sep=3.40 sep=3.45 sep=3.50 sep=3.55 sep=3.60 sep=3.65 sep=3.70 sep=3.75 sep=3.80 sep=3.85 sep=3.90 sep=3.95 sep=4.00

units lj
newton on
atom_style ellipsoid
dimension 3
boundary p p p
read_data config/$f
variable T index 1

pair_style gayberne 1.0 3.0 1.0 5.0
pair_coeff * * 1.0 1.0 1.0 1.0 0.2 1.0 1.0 0.2 5.0

comm_modify mode single vel yes
neighbor 1.9 bin #multi
neigh_modify every 1 delay 0 check yes page 5000000 one 500000

velocity all create $T 91259 loop geom
fix 2 all nve/asphere

compute q all property/atom quatw quati quatj quatk
compute sh all property/atom shapex shapey shapez

thermo_style custom step temp evdwl pe ke epair density

thermo_modify norm no
thermo_modify flush yes
thermo 1
log log/$g=.dat

dump 20npbc all custom 1 dump/dump$g.patch id type x y z c_q[1] c_q[2] c_q[3] c_q[4] c_sh[1] c_sh[2] c_sh[3]
dump_modify 20npbc colname c_q[1] quatw colname c_q[2] quati colname c_q[3] quatj colname c_q[4] quatk
dump_modify 20npbc sort id

timestep 0.001
reset_timestep 0
run 1

clear
next b
next g
next f
#next T
jump in.patch

config folder having files for particle separations as one file sep=1.00.dat

2 atoms

2 ellipsoids

1 atom types

0.0 10.0 xlo xhi
0.0 10.0 ylo yhi
0.0 10.0 zlo zhi

Atoms

1 1 1 0.4 0.0 0.0 0.0
2 1 1 0.4 1.00 0.0 0.0

Velocities

1 0.0 0.0 0.0 0.0 0.0 0.0
2 0.0 0.0 0.0 0.0 0.0 0.0

Ellipsoids

1 3.0 1.0 1.0 0.0 0.0 0.0 1.0
2 3.0 1.0 1.0 0.0 0.0 0.0 1.0

First of all, you should make the effort of formatting your input script.
Secondly, I don’t quite understand what are you trying to achieve. Please explain why you want to use data files with atoms in different positions. It’s a cumbersome approach, plus I have took the time to prepare a script that uses the fix move to achieve the same in a more elegant way.

Thanks for your reply. I will spend more time editing the script further.
I was using the read data file to simulate the rigid body of three atoms.

I was trying to make the correspondence relation in parameters (T in lammps and \epsilon_0 in vmmc that perform MC simulations) to get similar energies. For that, I started with two particles and only the Gayberne potential. So, I also use the data file reading for GB particles too.

Actually, I was not aware of the “fix move” command. Sorry for that, I will read all available commands.

Thanks, Hothello, for providing me with help!
Sorry to ask a small query again.

Previously, I was making a mistake in understanding the role of anisotropy parameters; if we take the orientation of particles along x or y, then the role of \epsilon_a, \epsilon_b, and \epsilon_c got changed. The orientation direction for ellipsoidal particles always represents the end-to-end interaction. e.g.,

For the shape coordinate, 3 1 1, the \epsilon_a represents the end-to-end interaction.
For the shape coordinate, 1 3 1, the \epsilon_b represents the end-to-end interaction.
For the shape coordinate, 1 1 3, the \epsilon_c represents the end-to-end interaction.

Is that a correct understanding?

I am still struggling to understand one more point. For the x or y orientation of particles for k’=5, the ratio of side-to-side and end-to-end energies are 5, but the relative well depth of both potentials are different compared to the orientation of particles along the z direction.

For clear visualization, I am attaching the figures for energy comparison when particle
orientations are taken along x, y, and z directions.



The relative difference is of factor 3 for shape 3 1 1
the factor of 5 for shape 4 1 1
the factor of 1.65 for shape 2 1 1
for completely matching the well-depth for any orientation of particles.

Please help me to understand the reason beyond this.