[lammps-users] NaN values while generating thermo output values

Hello,

I am using the LAMMPS version 29Oct2020. I have a model where I have a graphene layer placed in between 8 layers of copper on either side of it. I made it using Materials Studio using a cvff forcefield, and I used msi2lmp to generate a data file for it. I have a small LAMMPS code which does tensile loading using displacement control on it. But the simulation gives my ‘NaN’ values for the thermo_style output values. I know that this occurs due to the overlapping of atoms (one of the reasons). I tried reducing the time-step, but it did not work. Does anyone have any suggestions as to how I can overcome this issue? I am attaching the LAMMPS code and a picture of my model along with this. I have the data file, and if required, I will send it. Thank you!

Regards,
Rajesh

MD simulation Cu-Graphene tensile loading

units metal
dimension 3
boundary p p p

atom_style full
neighbor 2.5 bin
neigh_modify delay 5

pair_style lj/cut/coul/long 10.0
bond_style harmonic
angle_style harmonic
dihedral_style harmonic 3.0
improper_style cvff

read_data SquareModel.data

kspace_style ewald 1.0e-4
kspace_modify gewald 0.2428

Assign original velocities to atoms

compute new all temp
velocity all create 298.1 487639 temp new

Set up ensemble

fix 1 all nvt temp 298.1 298.1 100.0
fix 3 all temp/rescale 10 298.1 298.1 0.01 1.0

timestep 0.001

Output temperature, energy, pressure of the system

thermo_style custom step temp etotal pxx pyy pzz
thermo_modify flush yes
thermo 100
dump 1 all xyz 20000 dump.xyz
#dump_modify 1 element Si
dump_modify 1 element C Cu

compute str all stress/atom NULL

dump 2 all custom 50000 dump.*.stress id type x y z c_str[1] c_str[2] c_str[3] c_str[4] c_str[5] c_str[6]

delete_atoms overlap 0.1 all all
run 400000
quit

Fix rigid boundary atoms

compute new2 mobile temp
fix 2 boundary setforce 0.0 0.0 0.0
unfix 3
fix 4 mobile temp/rescale 100 298.1 298.1 0.01 1.0
fix_modify 4 temp new2

Apply displacement control loading

velocity up set 0.0 0.0 0.02 units box
velocity low set 0.0 0.0 -0.02 units box
velocity mobile ramp vz -0.02 0.02 z -70.0 70.0 sum yes

thermo_style custom step temp etotal pxx pyy pzz pxy pxz pyz
thermo_modify flush yes
thermo 100

run 600000

there is some other possible reason. you have this in your input:

fix 1 all nvt temp 298.1 298.1 100.0
fix 3 all temp/rescale 10 298.1 298.1 0.01 1.0

this is very bad. also those temp/rescale settings are such that it is almost useless (too frequent and to small a tolerance until rescaling happens).
fix temp/rescale is a very bad choice for a thermostat under almost all circumstances. but even then, having two competing thermostats on the same atoms is even worse.

You have just eliminated one problem, if you also have overlapping atom coordinates or inconsistent box dimensions, then you have to correct for that as well.
There is no simple workaround for that. MD follows the GI-GO principle, so you will only have well behaving and useful simulations if you eliminate all possible problems.

Hello,

Thank you for the response, Dr. Kohlmeyer. I have eliminated the fix temp/rescale step and run the code again. I am still getting ‘NaN’ values in my thermo output. I am stumped now.

Regards,
Rajesh