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.