Using replicate command applies forces on optimised structures

LAMMPS LAMMPS Beginners
1

Dear LAMMPS Users,

I hope this message finds you well.

I am a beginner with LAMMPS and currently attempting to replicate the results of a study that performed structure optimization and phonon dispersion calculations for naphthalene (J. Chem. Theory Comput. 2020, 16, 2716−2735).

For the structure optimization, I am using the primitive unit cell (1×1×1) of the experimental structure as the simulation box, applying periodic boundary conditions, and optimizing the lattice constants with the box/relax tri method. While the results do not match the values reported in the paper exactly, likely due to differences in force fields and limited information provided in the study, I have obtained very similar lattice constant values.

To calculate the phonon dispersion using the frozen phonon method, I created a supercell from the optimized primitive unit cell using the replicate command. However, I noticed that when expanding the cell in any direction, some atoms experience nonzero forces. These forces range from what I consider sufficiently converged values to about 0.1 kcal/Å.

From my understanding of physics, if the primitive unit cell is properly optimized, the resulting expanded supercell should maintain uniform structural optimization. While I expect minor errors during the expansion process, I find the observed magnitudes of force to be unexpectedly high.

I would like to seek advice from the community regarding the following questions:

  1. Is it appropriate to perform structure optimization for the primitive unit cell of molecular crystals using LAMMPS?
  2. If the answer to (1) is yes, why might forces appear on atoms when the optimized structure is expanded into a supercell? Are there any ways to address this issue?

Thank you very much for your time and insight. I greatly appreciate your support.

Best regards,
naphalene.input (1.4 KB)
naphtalene.lmp (20.1 KB)
log.lammps (8.1 KB)
opt_force_primitive_dump_final.dat (3.4 KB)
no_opt_force_no_pbc_10_10_10_dump_final.dat (6.4 KB)
S.H.

2

Using a single unit cell enforces symmetries that may be broken in a supercell. Whether that happens depends strongly on the force field in use.

Since LAMMPS uses floating point math, all operations are subject to floating-point rounding and precision issues, this is even more true for non-orthogonal cells, since many operations require conversion to fractional coordinates and back.

Here are some websites with discussions of the floating-point issues. Those are good to know about in general when doing simulations, not just for the topic at hand.