# [lammps-users] Problem with Coulomb interactions in deforming simulation boxes

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.

1. 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

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?

1. 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

The K-space grid is not adapted during either
a minimization or NPT run. The assumption
is that you are not changing the volume so dramatically
that it needs to be adapted. If that's not true,
then you should restart the simulation (e.g. 2 runs
in a row), in which case the K-space grid will
get recomputed (with the current volume). That
will introduce a (small) discontinuity in the energy function,
so you shouldn't expect to be able to minimize
consistently across that change. But you should have
no expectation that you can do a single minimize for
a large volume change anyway - that isn't a consistent
mathematical problem.

Steve

Thanks for the clarification.When I use fix deform I only change the dimension of the box by 1%, which I presume counts as a small amount, so that keeping the k-space grid constant should be okay.

However, here’s a word of warning for the record: there are some problems for which this isn’t accurate enough. In my case I’m calculating thermal expansion coefficients, which requires determining the lattice parameter to 4, preferably 5 significant figures. In the example I gave in my first email the lattice constant corresponding to the minimum energy (at T=0) is 5.123 when I deform from lz_lo to lz_hi, but it’s 5.124 when I deform from lz_hi to lz_lo, and this isn’t accurate enough. Clearly a further simulation is required, in which the deformation is restricted to a much smaller range to find the true minimum.

Philip

-------- Original-Nachricht --------

I'll just add that the minimize doc page lists some issues/settings that
can damage your ability to converge a minimization to hi accuracy. One
of them is lo tolerance in PPPM or Ewald, which may be an issue for you.

Steve