box/relax in long-range Coulomb systems

Dear all.

We do box/relax relaxations on a system with interactions described by a hybrid potential:
pair_style hybrid/overlay tersoff coul/streitz 14.0 ewald
kspace_style ewald 1e-6
kspace_modify gewald 0.3
Whenever there are nonzero charges on atoms (zero overall charge neutrality), pressure never reaches zero. CG minimization always ends with linesearch alpha is zero (this does not happen when all charges are set to zero). I found this old post (https://sourceforge.net/p/lammps/mailman/message/31580237/) which seems to talk about a similar problem. Of course, pressure is partly calculated from temperature and partly from virial. LAMMPS manual says that virial is everything without KE. Does it mean that the term /(3*V) arising from the Ewald summation is added as explained by Hummer et al., J. Chem. Phys. 109 (1998) 2791 in Eq. (15)? It is not included in the virial. I looked for the addition of this term in fix_box_relax.cpp and ewald.cpp but could not find it; perhaps it is done some other way. Sorry, if this is a silly question.

three comments:

  1. virial contributions from ewald summation is included: https://github.com/lammps/lammps/blob/master/src/KSPACE/ewald.cpp#L460

  2. there is no guarantee that a minimization will reach a zero pressure.

  3. shouldn’t coul/streitz be used in combination with fix qeq/slater?

If you do minimization, you are essentially doing a 0K calculation and thus temperature should be set to 0K so that you don’t get a bogus contribution to the pressure which could result in inconsistent results.

Axel.

Dear Axel.

Many thanks for your comments:

  1. Ok, I get it now. I read your paper and can see that the long-range electrostatic contribution is there. There seems to be an ambiguity about what is called virial. In your case, it is everything without the KE contribution. In the paper of Hummer, there is a virial contribution to pressure (KE part + short-range Coulomb part, i.e. ecoul) plus elong/(3*V).

  2. We are actually fitting the Tersoff part of the potential so that the sum of the Tersoff + long-range electrostatics reproduces experimental data. For each trial set of Tersoff parameters we check to what extent the cohesive energy, lattice parameters and elastic moduli agree with experiments and update these parameters (Nelder-Mead method). Pressure is never zero, although we were able to make it smaller by setting gewald to 0.3. Without the electrostatic part pressure is always zero. Should we also include the condition of zero pressure in the objective function?

  3. Yes, QEq is used: fix fq all qeq/slater 1 14.0 1e-8 1000 coul/streitz

We use molecular statics (minimize) and the temperature is thus automatically 0 K. Or, it has to be given somewhere explicitly?

Roman

Dear Axel.

Many thanks for your comments:

  1. Ok, I get it now. I read your paper and can see that the long-range electrostatic contribution is there. There seems to be an ambiguity about what is called virial. In your case, it is everything without the KE contribution. In the paper of Hummer, there is a virial contribution to pressure (KE part + short-range Coulomb part, i.e. ecoul) plus elong/(3*V).

The ambiguity is mostly due to you mistaking the keyword “virial” with the meaning of “what is a virial?”. The documentation clearly speaks of “kinetic virial contributions” but then consistently refers to the resulting property as “stress tensor” (which is then converted to pressure by dividing through the relevant volume).

  1. We are actually fitting the Tersoff part of the potential so that the sum of the Tersoff + long-range electrostatics reproduces experimental data. For each trial set of Tersoff parameters we check to what extent the cohesive energy, lattice parameters and elastic moduli agree with experiments and update these parameters (Nelder-Mead method). Pressure is never zero, although we were able to make it smaller by setting gewald to 0.3. Without the electrostatic part pressure is always zero. Should we also include the condition of zero pressure in the objective function?

As mentioned before, there is no guarantee that a minimization will lead to a zero pressure. A minimization will always lead to a local minimum, which may have some residual stress. Also a minimization will stop if it cannot find a pathway to a lower energy state. This can be different due to the minimization algorithm and/or the ruggedness/noisyness of the potential hypersurface. Minimization algorithms assume a smooth hypersurface for faster convergence. Also, sometimes you can get caught in a saddle point. The latter is a reason why people sometimes do a “stability test”, i.e. randomize positions by a suitably random margin and minimize again. The resulting symmetry breaking and displacement across (small) barriers could open a path to a lower energy state.

When you include charges, due to the long-range nature of their interactions and particularly in combination with charge equilibration, the risk of getting “trapped” or being “stuck” in minimization is increased.

  1. Yes, QEq is used: fix fq all qeq/slater 1 14.0 1e-8 1000 coul/streitz

We use molecular statics (minimize) and the temperature is thus automatically 0 K. Or, it has to be given somewhere explicitly?

Minimization itself will disregard kinetic energy contribution, but the pressure computation does not (by default), since the same code path is used for MD as well. The atom velocities will be set to 0.0 if not set explicitly, but if you read a restart or data file from an MD run, you will also import the velocities and to eliminate the the kinetic energy contribution to the stress, you need to either explicitly set the pressure compute for fix box/relax or explicitly set velocities to zero.

Axel.