Potential reasons for a "Master list cutoff distance =0"?

LAMMPS LAMMPS Beginners
1

I successfully ran my input scripts (was able to extract viscosity, density, and MSD values out of it) and tried to analyze RDFs, but I found that the RDF output file was empty. I looked at the logs and found that while rdf neighbor searching it hadn’t found any (I triple-checked the atom labels and they all seemed right), but noticed this log while running the RDF script saying that the master list cutoff distance = 0.

After extensively reading the documentation, I’ve deduced that it might have something to do with the pair styles or pair coefficients, but right now I’m at an impasse and would appreciate any help.

I am simulating a system of [EMIM][OAC] ionic liquid with a drude polarized force field.

(sorry, I’d attach the full input files but new users aren’t allowed to do so)
LAMMPS Version 29Sep2021

Prod run pertinent settings:
#Settings
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
pair_style hybrid/overlay lj/cut/coul/long 12.0 12.0 coul/long/cs 12.0 thole 2.600 12.0
pair_modify tail yes
kspace_style pppm 1.0e-5

read_data data-p.lmp extra/special/per/atom 6
include pair.lmp
include pair-drude.lmp

Coefficient File 1:
pair_coeff 1 1 lj/cut/coul/long 0.170000 3.250000 # NA NA
pair_coeff 1 2 lj/cut/coul/long 0.109087 3.396690 # NA CR
[continues…]

Coefficient File 2:
#interactions involving Drude particles
pair_coeff * 14* coul/long/cs
#Thole dipole-dipole damping if more than 1 Drude per molecule
pair_coeff 1 1 thole 1.208
pair_coeff 1 2 thole 1.321
[continues…]

RDF-MSD rerun full script:
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls

read_data data.run.lmp

#mean-squared displacements
group cat molecule 1:200
group ani molecule 201:400
compute MSDcat cat msd
compute MSDani ani msd

#radial distribution functions NA-NA O2-O2 NA-O2 HCR-O2 HCW-O2 HCW-O2
comm_modify cutoff 22.0
neigh_modify one 3000
compute RDF all rdf 200 1 1 6 6 1 6 7 6 4 6 5 6 cutoff 20.0
fix RDF all ave/time 2000 1000 2000000 c_RDF[*] file rdf.lmp mode vector

thermo 2000
thermo_style custom step c_MSDcat[4] c_MSDani[4]

rerun calcDump.lammpstrj dump x y z box yes

2

Exactly. You have no pair style defined at all. Thus there is no pairwise cutoff on record and the largest pairwise cutoff is used to initialize the neighbor list code. But since you specify an explicit cutoff in your rdf command and also provide a communication cutoff (needed for running in parallel), there should be no problem. If you absolutely feel that you need to set this, then you should add a pair_style zero command with a suitable cutoff. That won’t waste any time on doing a pair style force/energy/stress computation but otherwise you have a pair style present and an associated cutoff.

You can always upload files to dropbox, google drive, one drive, or similar and provide a link.

3

Hi,

Thanks for your response. I adeed pair_style and pair_coeff values, but I still get no RDF outputs, and I suspect its due to the following:

I read on the documentation that it has something to do with special_bonds:

If you have a bonded system, then the settings of special_bonds command can remove pairwise interactions between atoms in the same bond, angle, or dihedral. This is the default setting for the special_bonds command, and means those pairwise interactions do not appear in the neighbor list. Because this fix uses a neighbor list, it also means those pairs will not be included in the RDF. This does not apply when using long-range Coulomb interactions (coul/long , coul/msm , coul/wolf or similar. One way to get around this would be to set special_bond scaling factors to very tiny numbers that are not exactly zero (e.g. 1.0e-50). Another workaround is to write a dump file, and use the rerun command to compute the RDF for snapshots in the dump file. The rerun script can use a special_bonds command that includes all pairs in the neighbor list.

hence I put the line:

special_bonds lj/coul 1.0e-50 1.0e-50 0.5 angle yes dihedral yes

to include all pairs in the neighbor list. After doing this, I’m now able to find 1-3 and 1-4 neighbors, but 1-2 ones still aren’t showing up. Any thoughts anyone?

New rerun script:
units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 1.0e-50 1.0e-50 0.5 angle yes dihedral yes

read_data data.run.lmp
pair_style zero 20.0
pair_coeff * *

#mean-squared displacements
(from this point on, same script as posted above)…

Right, apologies this escaped me at the time of writing.

4

I don’t think that this is an issue of LAMMPS. The best way to check is to use a different tool to compute the g(r). Please check out the suggestions here: https://www.lammps.org/prepost.html

If you want somebody to “correct” the LAMMPS implementation, you need to provide some fast/easy to reproduce input example, e.g. based on one of the LAMMPS examples, and confirmation that other tools produce the expected results.