COMB3 energy drain

LAMMPS LAMMPS Mailing List Mirror
1

Hey folks,

I’ve been trying to simulate an interfacial system with the COMB3 potential, but have had some problems that I can’t seem to diagnose. For simplicity I recreated the problem with a small bulk alpha-alumina (Al2O3) system, with the initial data file and two input scripts attached.

The input in.COMB3.corundum_newCOMB3 is an input file I am using for a newly compiled version of LAMMPS (29 Mar 2019) which includes an updated version of the COMB3 forcefield files and fix qeq/comb source code from the group of Susan Sinnott https://research.matse.psu.edu/sinnott/software.

in.COMB3.corundum was used with LAMMPS(22 Aug 2018).

Regardless of whether I use the newly compiled version or the 22 Aug 2018 version I have the same issue.

Running an NVE simulation with COMB3 and fix qeq/comb the total energy decays extremely rapidly (~0 K in < 20 ps). While I know that due to the charge dynamics, the energy isn’t conserved, this seems to be something beyond a small drain on the energy due to the charge dynamics.

If I instead use

fix 2 all nvt temp 1043.0 1043.0 0.01

I instead see something different. For an NVT simulation, there is a buildup in the center-of-mass momentum, while the thermal motions go to zero. The thermostat keeps the temperature at 1043 K, but all of the kinetic energy is in center of mass translation.

If I use the same simulation setup for bulk aluminum instead of alumina, I don’t see this behavior. For an NVE simulation of bulk Al, there is a much slower decay of the total energy, more consistent with the loss of energy conservation due to the charge dynamics. In an NVT simulation, There doesn’t appear to be the buildup of total momentum, even over 200 ps or more.

As a final note, if I track the global scalar for the new version of fix qeq/comb from the Sinnott group, it fluctuates around zero.

An example thermodynamic output file is attached as comb3-3728228.

Any ideas for what might be going on here would be greatly appreciated.

Thanks,

Seth Martin

PhD Candidate

Department of Chemistry

The University of Kansas

data.small_corundum (10 KB)

in.COMB3.29Mar19.corundum (781 Bytes)

in.COMB3.22Aug18.corundum (785 Bytes)

comb3-3728228 (947 KB)

2

seth,

what you describe is rather well known and often referred to as “the flying ice cube syndrome”.
certain geometries are more sensitive to this, as they introduce some small asymmetry in the forces (sometimes due to numerics, sometimes due to weaknesses in the algorithm), but sometimes this is also induced by the equilibration procedure, which can initialize the system to have a small net momentum, which will then invariably grow and grow massively, as soon as (part of) the system starts to freeze out.

there is no general solution that solves all of these issues, but one has to determine the optimal approach for each system.
sometimes adding a very weak tether to some larger part of the system can remove this effect. sometimes it is sufficient to regularly remove/redistribute any net momentum. sometimes it is sufficient to remove it at initialization of the velocities. sometimes it can help to exclude a small number of atoms (<10) spread out through the system from time integration, so they remain immobile and all others interacting with them as well. if you search through the mailing list archive for “flying icecube”, you should find more discussions on the topic and probably more suggestions for workarounds and report of how well they work.

HTH,
axel.

3

You could try a different thermostat, like Langevin. I believe
it is less susceptible to the “flying ice cube” issue.

Steve