15
$\begingroup$

I have modeled a complex structure of a protein and a Carbon nano-cone (used as an indenter to the protein) using PACKMOL. The LAMMPS data file of the complex structure consists of 4 atom types (C O N S). But I constantly get an error saying 'Atoms lost' (I have tried with different potentials like ReaxFF, LJ, etc...). I know this means a few atoms have been blown off so far from the simulation box that they have been lost. LAMMPS mailing list archives all suggest for a bad geometry, but there is nothing I can do to modify the data file. Hence I really need some help with correcting the code(seen below). The main motive of the code is for the Carbon nano cone to indent the protein. The code below is the code I have written for LJ potential.

# -------------- Initialisation -------------------#

units           metal
boundary        p p p
atom_style      atomic
timestep        0.001

read_data       SE_atomic.dat
mass    1 14.0067   #Nitrogen
mass    2 12.0107   #Carbon
mass    3 15.9940   #Oxygen
mass    4 32.0650   #Sulphur


#Substrate

region           sars block INF INF INF INF INF 300 units box
group            sars region sars

#Indenter

region           top block INF INF INF INF 300 INF units box
group            top region top

#Potentials


pair_style       lj/cut 2.5

pair_coeff  1  1  0.0163859696   2.9565
pair_coeff  1  2  0.07146236771  3.1365
pair_coeff  1  3  0.07168744659  2.9165
pair_coeff  1  4  0.02856991424  2.8232
pair_coeff  2  2  0.12219124354  3.2615
pair_coeff  2  3  0.12257609881  3.0415
pair_coeff  2  4  0.04885079323  2.9482
pair_coeff  3  3  0.07146236771  3.1365
pair_coeff  3  4  0.07168744659  2.9165
pair_coeff  4  4  0.09181993247  3.0815


compute         peratom all pe/atom
neighbor        2.0 bin
neigh_modify delay 20 every 1 check yes page 500000 one 50000

#Initialisation

velocity all create 353 511124 rot yes mom yes
variable        t_step equal "step"
variable        t_temp equal "temp"
thermo          100
thermo_style    custom step press temp pe

#Relaxation

fix rigid sars addforce 0.0 0.0 0.0
fix relax top nvt temp 250 250 0.01
variable t equal 0.005
timestep ${t}
dump 1 all custom 500 equil.*.dump id type x y z fx fy fz

run 3000


#Indentation

fix 1 top addforce 0.0 -0.05 0.0
fix fix1 all nvt temp 250 250 0.01
dump 2 all custom 100 sars-Indent.*.dump id type x y z fx fy fz

run 1000


print "All done"

This code is written for ReaxFF potential

# -------------- Initialisation -------------------#

units           metal
boundary        p p p
atom_style      charge
timestep        0.001

read_data       SE_charge.dat
mass    1 14.0067   #Nitrogen
mass    2 12.0107   #Carbon
mass    3 15.9940   #Oxygen
mass    4 32.0650   #Sulphur


#Substrate

region           sars block INF INF INF INF INF 300 units box
group            sars region sars

#Indenter

region           top block INF INF INF INF 300 INF units box
group            top region top

#Potentials

pair_style      reax/c NULL
pair_coeff      * * SARS.ff N C O S
compute         peratom all pe/atom
neighbor        0.5 bin
neigh_modify delay 20 every 1 check yes page 500000 one 50000

#Initialisation

velocity all create 353 511124 rot yes mom yes
variable        t_step equal "step"
variable        t_temp equal "temp"
thermo          100
thermo_style    custom step press temp pe

#Relaxation

fix rigid sars addforce 0.0 0.0 0.0
fix relax top nvt temp 250 250 0.01
fix charge all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
variable t equal 0.005
timestep ${t}
dump 1 all custom 500 equil.*.dump id type x y z fx fy fz

run 2000


#Indentation

fix 1 top addforce 0.0 -0.05 0.0
fix fix1 all nvt temp 250 250 0.01
dump 2 all custom 100 sars-Indent.*.dump id type x y z fx fy fz

run 1000


print "All done"

If anyone in the community can help regarding the same, it would be really helpful for me! Cheers!

I have uploaded my LAMMPS datafile for the LJ case and my datafile for the ReaxFF case, and all data files, input files, and log files are located in this folder.

$\endgroup$
10
  • $\begingroup$ I wanted to share the data files too, but I am new to stack exchange and have no idea how upload the data file. Can you help me with that? $\endgroup$ Jun 25, 2020 at 2:37
  • $\begingroup$ Can you check how many lines it is? $\endgroup$ Jun 25, 2020 at 2:48
  • 2
    $\begingroup$ On many SE sites like TeX.SE, it is considered fine to include such data in the question if it's "fewer than a few hundred lines". Since you have 30,000, I will create a GitHub repository that you can upload the file to. Other SE sites have resorted to things like this as well. $\endgroup$ Jun 25, 2020 at 2:54
  • 2
    $\begingroup$ I created the repository, but while doing so and thinking about how it would work, it I remembered that I already had a repository for this type of thing, since long before this site existed. Please create a folder called Issue018 here: github.com/HPQC-LABS/QCBugs Then upload the data files, and possibly also the input files and output files, then create an "Issue" and perhaps just post a link to this question there, instead of explaining the bug over there. $\endgroup$ Jun 25, 2020 at 3:01
  • 1
    $\begingroup$ You should change the units to lj for the LJ input or to real for the reaxFF simulation. In LAMMPS the units you should use is given by the type of force field. $\endgroup$
    – Fabian
    Jun 25, 2020 at 7:01

1 Answer 1

14
$\begingroup$

The "Lost atoms" error typically happens when huge forces blow up an MD calculation. This is also the problem here:

The main issue are the units. The input parameters for reax force fields are given in "real" units by default. Unless they are converted real units have to be used in the lammps input file. Setting the units appropriately eliminates the "Lost atoms" error because the forces now have reasonable values.

Two more (not directly related) issue are standing out: The "fix nvt" and the "region" commands. Unless explicitly unset a fix persists until the end of the simulation. Therefore,

fix relax top nvt temp 250 250 0.01

and

fix fix1 all nvt temp 250 250 0.01

each apply a time integrator to all atoms in "top". This will certainly lead to undesired results.

The simulation box defined in "SE_charge.dat" only extends to 216, so the region top is always empty.

region           top block INF INF INF INF 300 INF units box

Therefore, the entire Relaxation block does nothing because the nvt time integrator is applies to 0 atoms.

TL;DR With a few tweaks to the input file the simulation with the reaxFF is running. All I changed is

# The units for the parameters in reaxFF input decks are real
units           real                 # time is given in fs

# Your simulation box only extends to 216.506012, setting 300 in the commands below leaves the region top empty. Set this value to something usefull
region           sars block INF INF INF INF INF 200 units box
region           top block INF INF INF INF 200 INF units box

# adjust the time step to the new units
fix relax top nvt temp 250 250 100
variable t equal 1

# use unfix to remove the nvt fix, otherwise you have two time integrators active!
unfix rigid
unfix relax

# again adjust the time step
fix fix1 all nvt temp 250 250 100

EDIT: The entire input file with changes is here.

$\endgroup$
6
  • $\begingroup$ Hey thanks a lot for taking time to help me with this! would it be okay if you could mail me the input file ( [email protected]) ? $\endgroup$ Jun 25, 2020 at 8:36
  • 1
    $\begingroup$ Can you instead push it to github so that we can close the issue there? Let me know your username. But Sai, doesn't this answer already say what exactly to change? $\endgroup$ Jun 25, 2020 at 13:34
  • $\begingroup$ @NikeDattani I have added a link to the edited input file $\endgroup$
    – Fabian
    Jun 25, 2020 at 15:44
  • $\begingroup$ The code runs smoothly, but as I observed the dumpfile animations in OVITO the nano cone does not indent the protein at all. The tip of the nano cone just splits for some reason. I have been trying to get it indented from past 9 hrs (morning till evening here in India) but no useful result. $\endgroup$ Jun 26, 2020 at 3:33
  • $\begingroup$ Sai, is that a new question? Or part of the question of the input file not working? $\endgroup$ Jun 26, 2020 at 13:09

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .