Dear LAMMPS community,

I am using LAMMPS to simulate clay systems using ellipsoidal particles and the Gay-Berne potential. All of my simulations are run with periodic boundary conditions. After an equlibration phase, virtual samples are subjected to compression.

I initially performed isotropic compression tests controlling the external applied pressure using the command fix npt/asphere temp Tstart Tstop Tdamp iso Pstart Pstop Pdamp. My temperature slightly oscillated around the target value and the pressure smoothly reached the prescribed value. I printed the stress tensor of my system and, as expected, that was isotropic with Pxx=Pyy=Pzz. On the side, to better understand LAMMPS and the commands I was using, I developed a post-processing MATLAB code which, starting from LAMMPS output on particles positions and orientation, identifies all interacting pairs, computes the GB potential for each interacting pair and the components of the force (according to documentation of GB and source code). With this script I was able to compute the stress tensor of the system which was exactly like that computed by LAMMPS.

I moved then to 1D compression tests (with the aim of replicating some lab tests). I performed isotropic compression up to a point where particles were contacting and then I applied a uniaxial compression only in the z direction, maintaining fixed the x and y directions. I employed the command fix npt/asphere temp Tstart Tstop Tdamp z Pstart Pstop Pdamp (with and without nreset keyword). I monitored the length of the box in the three directions and I saw that, correctly, z dimension smoothly decreased and x and y dimensions remained fixed throughout the simulation. The temperature and the total energy (potential + translational and rotational kinetic) were monitored as well and no unexpected behaviour was observed. By visualising my system in OVITO I could not spot anything wrong as the configuration of the particles is as I expected it to be. However, when checking the pressure tensor given by LAMMPS I noticed that it was still isotropic with Pxx=Pyy=Pzz, while I was expecting something anisotropic, given I was compressing only in one direction maintaining the others "constrained". When checking with my MATLAB script, I was able to compute the same potential energy as given by LAMMPS at the end of the run. However, when computing the stress tensor (similarly to what I did for 3D compression), I obtained, as expected, an anisotropic stress state with Pzz>Pxx=Pyy (specifically Pzz~2*Pstop and Pxx=Pzz~0.5*Pstop so that (Pxx+Pyy+Pzz)/3=Pstop).

Has someone experienced this issue before with non-spherical particles and Gay-Berne potential when performing 1D constrained compression?

Thanks very much in advance.

Kind regards,

Sara