error with conserved quantity when using fix external pf/callback and fix NPT or MSST

Hi All,

I am using pf/callback to run MD simulations with DFTB+ (written in fortran). I’ve had success with running NVE and NVT simulations, where the conserved quantity oscillates about a constant value given a reasonable time step, etc. However, for flexible cell simulations such as NPT or MSST, I notice discontinuities in the conserved quantity after a small number of time steps (< 100). These discontinuities are not present when running the same system (and virtually the same input file) with an earlier version of the DFTB+/LAMMPS hybrid, where I used a dftb pair style that writes out DFTB+ input files, uses a system call to run DFTB+, and then reads in the total energy, atomic forces, and stress tensor from file. I should also point out that in the DFTB+ callback function, the box dimensions are set by using call lammps_extract_global(boxxlo, lmp, 'boxxlo’), etc. Atomic coordinates in DFTB+ are indexed by the LAMMPS id array that is passed in to the callback function. The potential energy from DFTB+ is set using a wrapper function to FixExternal::set_energy, i.e., call lammps_set_user_energy (lmp, etotal).

Below is an example input file for this hybrid code. Note that the timestep and total number of MD steps are set in DFTB+ itself using calls to the lammps_command function.

Any thoughts from anyone?

Thanks,

Nir

units real

atom_style charge

atom_modify map array

atom_modify sort 0 0.0

read_data data.lammps

restart 100 h2o2_9kms.restart

neighbor 1.0 bin

neigh_modify once yes

velocity all create 300.0 9999

fix 2 all external pf/callback 1 1

fix_modify 2 energy yes

fix msst all msst x 0.09 q 500 mu 0.3 e0 -292964.884961449 p0 5.0e3 tscale 0.01

fix_modify msst energy yes

variable dhug equal f_msst[1]

variable dray equal f_msst[2]

variable lgr_vel equal f_msst[3]

variable lgr_pos equal f_msst[4]

thermo_style custom step temp ke pe lx ly lz press pxx pyy pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos f_msst

thermo 1

dump 1 all custom 5 pos.xyz id type q x y z

dump 2 all custom 5 vel.xyz id type q vx vy vz

If you’re saying the only differnence in a right

answer vs a wrong answer (no energy discontinuity

vs yes), is the method of hooking LAMMPS to

DFTB+, then it should be relatively easy

to debug. Just run both versions and dump

out everything you can every timestep from

both LAMMPS and DFTB+. You should

be able to spot when something is different,

either the values DFTB+ is using as input,

or the values LAMMPS is using to perform

NPT or MSST and do its timestep.

I also suggest you write out all your debug
info (and the file version input values to DFTB+)

at high precision to avoid round-off drift (when

the file values are read back in).

Steve

Hi Nir,

From your description, I can’t tell whether the atom positions are being rescaled along with the box dimensions. If not, this would cause the barostat to misbehave.

Aidan

Hi Nir,

From your description, I can't tell whether the atom positions are being
rescaled along with the box dimensions. If not, this would cause the
barostat to misbehave.

this reminds me of something.

aidan, didn't you some time ago change the default reference point for
the atom position rescaling from the origin to the center of the box?
and make it configurable from the input?
since there are no specifics, it is impossible to say, if this
happened before of after the "earlier version" that was mentioned. it
has been a while, too.

axel.

We added the option of making the fixed point configurable in 2012, but the default behavior is the same as what it has been for a long time, the center of the box remains fixed. Under normal conditions, using a different fixed point should not affect the energy at all.

Thank you everyone for your input. FWIW, the culprit seems to have been a
subtle issue with how DFTB+ performs updates for the reciprocal part of
the Ewald sum. As a quick plug, I believe I have a working version of
DFTB+ that accesses LAMMPS functionality through pf/callback, if any of
the LAMMPS usership is interested.

Best,

Nir