[lammps-users] Issue with 1D compression and Gay-Berne potential

LAMMPS LAMMPS Mailing List Mirror
1

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

2

Thanks for a careful explanation of the issue. I don’t know of any
particular reason why LAMMPS vs your code would disagree when
a system is compressed in one dimension. Also if the system is liquid-like
it may not sustain large pressure differences in x,y,z.

I think you are saying you get an exact match (to high precision?)
between your MATLAB post-processing code and LAMMPS with no compression. I assume
by Pxx = Pyy = Pzz you mean approximately equal in an instantaneous sense.

To debug this I suggest the following. For as small a system as possible,
write out a data file (and dump file to see if they are consistent) for an early
time when the 2 codes agree, and a late time (compressed) when they
do not agree. Verify that if LAMMPS reads the data file it produces the
same thermo output (pressure components, etc) on the new step 0 as at the end of
the previous run. Then you can just do a 0-step run with the new data file.

The pressure (virial) for both codes is just a summation
of pairwise terms. If the 2 codes disagree, you should be able to identify
a pairwise term that is different between the 2 codes. Since the virial
for a pair of atoms is just the outer product of the delr and f 3-vectors,
you can print out those 6 values from your code and from the inner loop
of the PairGayBerne::compute() method, near the botttom before
the call to ev_tally_xyz(). The 6 components of the delr and f 3-vectors
are arguments to that method.

Steve

3

Dear Steve,

Thanks so much for your reply. I’ll try to debug the code as you suggested.

Kind regards,

Sara