Syntax error while transcribing simulation trajectories

LAMMPS LAMMPS Beginners
1

Dear LAMMPS users,

I am attempting to post-process some GCMC simulations generated with the DL_MONTE simulation software using LAMMPS’ rerun feature. To do this, I use a python script to transcribe DL_MONTE output files into LAMMPS dump files.

The simulation is of gas adsorption in a MOF, and my only pairwise interactions of interest are intermolecular gas-gas and gas-MOF interactions. To this end I have attempted to switch off any unwanted interaction pairs with the neigh_modify exclude command:

  • neigh_modify exclude type I J for unwanted MOF-MOF interactions, and
  • neigh_modify exclude molecule/intra gas for intramolecular gas-gas interactions.

While the former exclude statements appear to be working, my simulated interaction energies are not changing when the latter statements are turned on/off. I believe this is an issue with how I’ve transcribed different molecule types into the dump file. Please could you take a look at the following MWE and suggest where in the dump or input files my molecule definition is going wrong?

Thanks,
Joe Manning
University of Manchester, UK

Minimum working example

Input file:

# 1) Initialization
units real
dimension 3
atom_style full
pair_style lj/cut 12.5
boundary p p p

# 2) System definition
region simulation_box block 0.0 25.832 0.0 25.832 0.0 25.832
create_box 3 simulation_box

# 3) Simulation settings
# Acetonitrile atom type definitions: 1 = C, 2 = Me, 3 = N
mass 1 12
mass 2 17
mass 3 14

pair_coeff 1 1 0.1192 3.55
pair_coeff 1 2 0.1192 3.25
pair_coeff 1 3 0.1524 3.65
pair_coeff 2 2 0.1192 2.95
pair_coeff 2 3 0.1524 3.35
pair_coeff 3 3 0.1947 3.75


# 4) grouping atoms for processing
group gas type 1 2 3
neigh_modify exclude molecule/intra sorb

# 5) outputting dummy data
thermo_style   custom step temp etotal
thermo 1

rerun acn_mwe.dump every 1 dump x y z purge yes add yes replace no format native

Dump file:

ITEM: TIMESTEP 
0
ITEM: NUMBER OF ATOMS 
6
ITEM: BOX BOUNDS pp pp pp
0.0 25.832
0.0 25.832
0.0 25.832
ITEM: ATOMS id type x y z q mol
N  3 15.9568715 23.164098  10.9439578 0.269  1
C  1 14.9471131 23.5570332 10.5382007 0.129  1
Me 2 13.6030958 24.0800411 9.9981265  -0.398 1
N  3 10.1803461 23.3906558 9.1216064  0.269  2
C  1 9.3700299  22.6295461 8.8010619  0.129  2
Me 2 8.291476   21.6164873 8.3744081  -0.398 2

Expected versus actual behaviour

From the thermo command in the input file, I expect the following to be printed to stdout:

Step    Temp    TotEng    
0       0       -0.13326

N.B. This energy value was calculated manually in python

Instead, I get this:

Step    Temp    TotEng    
0          0          1190783.9

System info:

OS: Windows 10
CPU: Intel i7-10750H @ 2.60GHz
LAMMPS version: LAMMPS 64 bit 2 Aug 2023
No MPI version

2

By default, LAMMPS does not evaluate pair interactions for 1-2 and 1-3 neighbours. For documentation and how this default can be changed, see special_bonds.

If your acetonitrile molecule is a three-bead linear chain united atom model and you haven’t changed the defaults, then there is no intramolecular pair contribution to exclude in the first place, and your results make sense.

1 Like
3

Dear Dr. Tee,

Thank you for your quick reply. I believe you’ve correctly identified my issue - bonding information has not been correctly transcribed along with the trajectories (hence the clearly unphysical potential energy of ca. +1.2 x107 kCal/mol from my earlier example).

To resolve this, I’m now attempting to assign bonds on the fly to the trajectory as read in by rerun, as the documentation indicates that read_dump cannot read bonding information from the trajectory files. However, while attempting to insert the following code block into my runfile I am running into a second issue: although the bonds are now defined, they are not being applied to newly-inserted atoms as they’re imported. (N.B. the simulation still runs to completion with identical outputs to my original example).

bond_style zero nocoeff
bond_coeff  1 
create_bonds many all all 1 0.5 1.8

I guess this leaves me with 3 options:

  • reformat the dump file so that rerun can read in bonding information
  • modify the runfile somehow so create_bonds is affecting the imported atoms
  • return to plan A and exclude all LJ/Coulombic interactions below a certain distance threshold e.g. with neigh_modify collection/interval

Please could you advise to let me know if any of the above 3 options sounds feasible?

Thanks,
Joe Manning