Python API fix external virial error

LAMMPS LAMMPS Development
1

I’m using fix external with the virial option enabled in the Python API while running NPT simulations. If I run a simulation for 10 steps, close the lammps object, free the communicator then recreate the simulation and repeat several times, sometimes the simulation pressure will become unstable. I have an example where this happens:

from mpi4py import MPI
from lammps import lammps
import numpy as np

allcomm = MPI.COMM_WORLD
np_global = MPI.COMM_WORLD.Get_size()
rank = allcomm.Get_rank()

# split into one group
color = 0
comm = allcomm.Split(color)

cmds = [
    "units metal",
    "atom_style full",
    "boundary p p p",
    "read_data coord.lmp",
    "change_box all triclinic",
    "include ff.lmp",
    "fix md all nph iso 1 1 1 tchain 5 pchain 5 mtk yes",
    "fix tst all temp/csvr 300 300 0.1 20384",
    "timestep 0.001",
    "velocity all create 300 30094 mom yes dist gaussian",
    "fix com all momentum 100 linear 1 1 1",
    "fix ext all external pf/array 1",
    "fix_modify ext virial yes",
    "run 0 post no",
    "run 10 post no",
]

lmp = lammps(cmdargs=["-screen", "none", "-log", "log.first-" + str(color)], comm=comm)
lmp.commands_list(cmds)
lmp.close()
comm.Free()

for _ in range(200):
    comm = allcomm.Split(color)
    lmp = lammps(
        cmdargs=["-screen", "none", "-log", "log.second-" + str(color)], comm=comm
    )
    lmp.commands_list(cmds)
    lmp.close()
    comm.Free()

lmp.finalize()

ff.lmp (3.2 KB)
coord.lmp (88.8 KB)

After running for long enough I get the following:

ERROR: Non-numeric pressure - simulation unstable (src/fix_nh.cpp:1068)

with the internal pressures returning nan. My LAMMPS version is 20230208. I’ve had a similar problem to this before when the virial array wasn’t being initialised with 0s so there were issues with memory allocation, but that was solved in a different patch. Any help is appreciated.

Cheers,
Blake

2

The problem is that you tell LAMMPS that there will be an external force and virial set, but you never set it. If you change your input as follows, this is addressed:

from mpi4py import MPI
from lammps import lammps
import numpy as np

allcomm = MPI.COMM_WORLD
np_global = MPI.COMM_WORLD.Get_size()
rank = allcomm.Get_rank()

# split into one group
color = 0
comm = allcomm.Split(color)

cmds = [
    "units metal",
    "atom_style full",
    "boundary p p p",
    "read_data coord.lmp",
    "change_box all triclinic",
    "include ff.lmp",
    "fix md all nph iso 1 1 1 tchain 5 pchain 5 mtk yes",
    "fix tst all temp/csvr 300 300 0.1 20384",
    "timestep 0.001",
    "velocity all create 300 30094 mom yes dist gaussian",
    "fix com all momentum 100 linear 1 1 1",
    "fix ext all external pf/array 1",
    "fix_modify ext virial yes",
    "run 0 post no",
]

lmp = lammps(cmdargs=["-screen", "none", "-log", "log.first-" + str(color)], comm=comm)
lmp.commands_list(cmds)
force = lmp.fix_external_get_force("ext");
nlocal = lmp.extract_setting("nlocal");
for i in range(nlocal):
    force[i][0] = 0.0
    force[i][1] = 0.0
    force[i][2] = 0.0
lmp.fix_external_set_energy_global('ext', 0.0)
lmp.fix_external_set_virial_global('ext', [0.0,0.0,0.0,0.0,0.0,0.0])
lmp.command("run 50 post no")
lmp.close()
comm.Free()

for _ in range(200):
    comm = allcomm.Split(color)
    lmp = lammps(
        cmdargs=["-screen", "none", "-log", "log.second-" + str(color)], comm=comm
    )
    lmp.commands_list(cmds)
    force = lmp.fix_external_get_force("ext");
    nlocal = lmp.extract_setting("nlocal");
    for i in range(nlocal):
        force[i][0] = 0.0
        force[i][1] = 0.0
        force[i][2] = 0.0
    lmp.fix_external_set_energy_global('ext', 0.0)
    lmp.fix_external_set_virial_global('ext', [0.0,0.0,0.0,0.0,0.0,0.0])
    lmp.command("run 50 post no")
    lmp.close()
    comm.Free()

lmp.finalize()
3

Hi Axel,

Thank you for your answer, however, even after running your modified code I still get the same problem. At some random point the simulation crashes:

ERROR: Non-numeric pressure - simulation unstable (src/fix_nh.cpp:1068)
Last command: run 0 post no

Printing the pressures out within lammps shows:

pressure: -nan -nan -nan

As it is a run 0 command that is failing upon initialisation, I’m led to believe that it’s still a memory allocation problem.

Cheers,
Blake

4

I disagree. The LAMMPS part is clean. I cannot find anything with a variety of debugging tools. So you would have to provide proof to support your claim.