Hi,

I’m trying to determine the ground-state lattice parameters of a tetragonal crystal (ax = ay < az). I’m modelling the interactions with a Buckingham potential, i.e. exp(-c*r) + A/r^6 + Coulomb, and I’m using kspace_style pppm 1e-6 to calculate the Coulomb interactions with periodic boundary conditions.

I find that the potential energy depends on the initial lattice parameter. Since I didn’t have this problem with other non-Coulombic potentials, it makes me wonder whether the k-space vectors are being properly adapted as the size of the simulation box changes. I’m using LAMMPS from 1 May 2010, which is after the bug fix on 27 April that claimed to solve this problem for minimize. I’ve tried two approaches, but they both exhibit the same problem.

- minimization with variable box size

Here are the relevant parts of the input file:

variable ax equal 5.082

variable ay equal 5.082

variable az equal 5.126

variable rel equal {az}/{ax}

lattice custom {ax} a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 {rel} origin …

…

pair_style buck/coul/long ${cutoff_dist}

kspace_style pppm 1.0e-6

…

fix RELAX all box/relax x 0.0 y 0.0 z 0.0 couple xy

min_style sd

min_modify line quadratic dmax 0.0001

minimize 0.0 0.0 1000 10000

unfix RELAX

This converges to the lattice parameters ax = ay = 5.081999, az = 5.088998. However, if I use az=5.128 as the starting value the same script doesn’t fully converge, but oscillates between ax=ay=5.082151, az=5.088490 and ax=ay=5.0823169, az=5.0886607. In all cases the final total potential energy is different. Surely ‘minimize’ should converge to (exactly) the same energy and lattice parameter?

- fix deform for the z-direction

I keep ax=ay fixed and vary the box-size in the z-direction. The variation of the potential energy shows a nice minimum. However, when I start at az_lo and increase to az_hi I get a different curve for when I start at az_hi and decrease to az_lo. The two potential energy curves are shifted sideways relative to each other, so the minimum energy occurs at a different value of the lattice parameter. Here are the relevant parts of the input file:

variable ax equal 5.082

variable ay equal 5.082

variable az equal 5.126

variable z_final equal 5.070

variable rel equal {az}/{ax}

lattice custom {ax} a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 {rel} origin …

…

pair_style buck/coul/long ${cutoff_dist}

kspace_style pppm 1.0e-6

…

variable lz_final equal {z_final}*{z_size}

fix DEFORM all deform 1 z final 0 ${lz_final} units box

run ${run_deform}

unfix DEFORM

I interpret the results as implying that the kspace-grid was set at the beginning of the deformation and not subsequently changed. However, if it is being properly adapted, why don’t I get the same energy for a given intermediate lattice parameter irrespective of the direction of deformation? How confident can I be that the k-space grid will be properly adapted when I move on to NPT simulations?

Thanks for any help,

Philip