We aim to simulate a two-dimensional (2D) system of ellipses interacting via the Gay-Berne (GB) potential using the LAMMPS package. Each ellipse in the system has an aspect ratio of 3:1. We also impose an interaction strength ratio of 5 between the side-by-side and end-to-end configurations. In terms of the standard GB parameters, we intend to simulate a (3, 5, 2, 1) GB system.
Monte Carlo Implementation
Before implementing this in LAMMPS, we developed a simple Monte Carlo (MC) code to simulate the system. We used two different representations of the GB pair potential:
- Vector Approach (VA) [1]
- Scalar Approach (SA) [2]
From an implementation standpoint, we found the SA significantly easier to code. However, the VA appeared more elegant and structured once completed.
After fixing the GB model parameters to (3, 5, 2, 1), the following matrices play a key role in the VA formulation:
sig_x = 3.0; sig_y = 1.0; sig_z = 1.0; // Shape parameters for (3,1,1) GB ellipsoids
eps_x = 0.2; eps_y = 1.0; eps_z = 1.0; // Energy parameters for (3,1,1) GB ellipsoids
// Shape anisotropy matrix (Eq. 64 in [1])
S[3][3] = {
{sig_x, 0, 0},
{0, sig_y, 0},
{0, 0, sig_z},
};
// Energy anisotropy matrix (Eq. 65 in [1])
E[3][3] = {
{1.0 / sqrt(eps_x), 0, 0},
{0, 1.0 / sqrt(eps_y), 0},
{0, 0, 1.0 / sqrt(eps_z)},
};
Matching Energy Between VA and SA
In implementing the VA, we found that the original expression for the energy strength coefficient epsA from Eq. (60) in [1]:
epsA = (sig_x * sig_y + sig_z * sig_z) * sqrt(2 * sig_x * sig_y / detA);
did not yield consistent pair interaction energies compared to the SA in 2D. However, when we modified the expression slightly to:
epsA = (sig_x * sig_x + sig_y * sig_z) * sqrt(2 * sig_y * sig_z / detA);
the energy matched between both approaches.
LAMMPS Configuration
For the LAMMPS simulation, we used the following commands in the input script:
units lj
atom_style ellipsoid
dimension 2
set type 1 shape 3 1 1
pair_style gayberne 1 1 2 4
pair_coeff * * 1 1 0.2 1 1 0 0 0 4
However, to ensure consistency with the corrected expression for epsA , we made changes in the LAMMPS source code:
- For CPU:
pair_gayberne.cppin/src/ASPHERE/ - For GPU:
pair_gayberne_gpu.cppin/src/GPU/
Suggested Modification to LAMMPS Source Code
This message is akin to a bug submission report concerning the current implementation of the Gay-Berne potential in the LAMMPS ASPHERE package.
In 2D simulations, the GB particle should ideally exhibit shape anisotropy as either (3,1,1) or (1,3,1) — not (1,1,3). However, we observed that LAMMPS currently appears to be implemented only for the (1,1,3) shape-anisotropy case when calculating the strength anisotropy term \eta_{12}, Eq. 11 in [3] (or equivalently Eq. 60 in [1]).
As such, we propose the following correction to the code section that computes the lshape parameter, to ensure correct computation of \eta_{12} for 2D simulations. This correction is based on Eq. 12 in [3], which corresponds to the appropriate 2D formulation of the strength anisotropy term.
The current line in the source code:
lshape[i] = (shape1[i][0]*shape1[i][1] + shape1[i][2]*shape1[i][2]) *
sqrt(shape1[i][0] * shape1[i][1]);
should be modified to (for (3,1,1) GB ellipses):
lshape[i] = (shape1[i][0]*shape1[i][0] + shape1[i][1]*shape1[i][2]) *
sqrt(shape1[i][1] * shape1[i][2]);
We kindly request the LAMMPS development community to review this suggestion and provide their feedback. We believe this modification will help ensure that 2D simulations of GB ellipses (with (3,1,1) shape anisotropy) using the LAMMPS ASPHERE package yield physically consistent results.
References
- Ricci, M., Roscioni, O. M., Querciagrossa, L. & Zannoni, C. MOLC. A reversible coarse grained approach using anisotropic beads for the modelling of organic functional materials. Phys. Chem. Chem. Phys. 21, 26195–26211 (2019).
- Gay, J. G. & Berne, B. J. Modification of the overlap potential to mimic a linear site–site potential. J. Chem. Phys. 74, 3316–3319 (1981).
- Everaers, R. & Ejtehadi, M. R. Interaction potentials for soft and hard ellipsoids. Phys. Rev. E 67, 041710 (2003).