# LAMMPS code that I wrote not behaving as expected. Some issue with the group & region commands

I am actually working on SARS-CoV-2 proteins. Specifically I am trying to tribologically disengage the Spike Glycoprotein from the Membrane protein using LAMMPS. The method is to use amorphous carbon as an abrasive to rub (mechanically wear) the contact point of both the above mentioned proteins.

The detailed explanation of my code (because I might have erred anwhere)

1. How I generated the datafile :
-I downloaded the PDB (Protein Database) of both the proteins spike and membrane.
-Used atomsk to convert the PDB files to lammps data files. Used the command atomsk protein_name.pdb lammps.
-This resulted in 2 LAMMPS datafiles. But for the Membrane protein, instead of one molecule of the protein, 4 molecules packed into a Parallelogram like simulation box was generated
I don't know why that happened but I was pretty happy that it did because I had to do build a layer of the Membrane protein anyway.
This is an image of the same Membrane protein before putting through Atomsk.

-You can observe at the top right corner area of each image to see that Bonds checkbox has disappeared from the OVITO visualization after the pdb file was run through Atomsk and LAMMPS datafile was generated. Again I have no idea about why that happenned!
-Next, I "placed" the spike glycoprotein on the parallelogram membrane layer using PACKMOL. I am actually worried about the junction where the spike protein meets the membrane protein because of manual placement of the proteins.
The code for Placing the spike on the membrane is as follows :

tolerance 3.0
filetype pdb
output E_S.pdb

structure M.pdb
number 1
inside box -44. 0. 0. 156. 50. 130.
center
fixed 0. 0. 0. 0. 0. 0.
end structure

structure S.pdb
number 1
inside box -44. 50. 0. 100. 200. 150.
center
fixed 0. 116. 0. 4.71238898 0. 0.
end structure


-After all this steps the model looked something like this:

-When visualized in Pymol the spike glycoprotein on the membrane is very clearly visible.

2. LAMMPS CODE

# Tribological detachment of SARS-CoV-2 Spike Glycoprotein.

# definition
units           real                                                        # types of units used
dimension       3                                                           # Defines a 3D simulation
processors      * * *                                                       # Command for optimum usage of processors
boundary        p p p                                                       # Defines periodic boundary conditions
atom_style      charge                                                      # Defines atom type to be charge

# SARS Input
mass            1 12.0107                                                   #Carbon
mass            2 14.0067                                                   #Nitrogen
mass            3 15.9940                                                   #Oxygen
mass            4 32.0650                                                   #Sulphur

# Group the abrasive atoms
region          carbonatoms block -96 -79 23 34 -8 11                                        # Defines a 3D block region called "abrasive" that is made up of the unit cell in the x-, y-, and z-direction for the given dimensions
group           carbonatoms region carbonatoms                                          # Assigns the name carbonatoms to atom type 1.

# Group the  SARS data file
region          sars block -100 100 -26 185 -65 65 units box                # Create a region for the datafile
group           sars region sars                                            # group the input file with the name "sars"

# Interatomic potentials
pair_style      reax/c NULL                                                 # Pair potential style ReaxFF
pair_coeff      * * ffield.reax.FC C N O S                                  # Assign Respective atoms

# Settings
compute         peratom all pe/atom                                         # Compute potential energy per atoms
neighbor        2.0 bin                                                     # NEVER KNEW WHAT THIS IS !!!!!!!!!
neigh_modify    delay 20 every 1 check yes page 500000 one 50000            # Helps with the lost atoms error!!!

# Initialization
velocity        all create 350 123456                                       # Setting temperature to 350 K
variable        t_step equal "step"                                         # Assigning a variable for step
variable        t_temp equal "temp"                                         # Assigning a variable for temperature
thermo          100                                                         # Show [#] for every 100 steps
thermo_style    custom step press temp pe                                   # Show [temp and step & PE ]

# Relaxation
#fix             rigid sars addforce 0.0 0.0 0.0                            # Make the sars group immobile by reducing force to 'zero' in all direction
#fix             relax carbonatoms npt temp 250 250 0.5 iso 0 0 0.5 drag 1  # The start and end temperatures of abrasive is 250K and start and end pressure of abrasive is 0 and the abrasive is damped.
fix             charge all qeq/reax 1 0.0 10.0 1.0e-6 reax/c                # Fix needed for reaxFF. DONNO WHAT THIS DOES !!!!!!
variable        t equal 1                                                   # Assigning timestep
timestep        ${t} # Equating timestep dump 1 all custom 200 equil.*.dump id type x y z fx fy fz # Dump 'atom id' 'atom type' 'x,y,z coordinates' 'force along x,y,z coordinates' run 1 # Run for 2000 time steps # always remove fixes that are no longer needed #unfix rigid #unfix relax #unfix charge #Scratching fix 2 carbonatoms move linear 5 0 0 units box # Apply force in x direction for wear @ 5 Angstroms/femtosecond dump 3 all custom 100 sars-Scratch.*.dump id type x y z fx fy fz # Dump 'atom id' 'atom type' 'x,y,z coordinates' 'force along x,y,z coordinates' run 1000 # Run for 2000 time steps # End simulation print "All done" # HOPE TO SEE THIS AT THE END OF SCREEN  -Since I usually have a lot of issues with LAMMPS (I'm still in the learning curve) in the code, I have mentioned what I think that command does in my code right next to the command itself. Please correct me if I'm wrong.\ 2.1 Regarding LAMMPS data file (Lipid.dat in the code) -When visualized it looks something like this. The individual atom pointed out is one of the many amorphous carbon atoms likely to be used as the abrasive for tribological disengagement. I have used just one atom just to see if the code works right.The coordinates of the pointed out atom is written is inside the data file itself (Refer Atom number 51227 in line number 51236 of the data file(lipid.dat)) Also the Bonds checkbox is missing for some reason when visualized 2.2 Results I'm getting:) OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:94) using 1 OpenMP thread(s) per MPI task # Tribological detachment of SARS-CoV-2 Spike Glycoprotein. # definition units real # types of units used dimension 3 # Defines a 3D simulation processors * * * # Command for optimum usage of processors boundary p p p # Defines periodic boundary conditions atom_style charge # Defines atom type to be charge # SARS Input read_data lipid.dat # Reads the datafile orthogonal box = (-99.923 -24.73 -64.058) to (99.923 184.496 64.058) 1 by 1 by 1 MPI processor grid reading atoms ... 51227 atoms read_data CPU = 0.091445 secs mass 1 12.0107 #Carbon mass 2 14.0067 #Nitrogen mass 3 15.9940 #Oxygen mass 4 32.0650 #Sulphur # create the abrasive atoms #lattice diamond 3.57 # Defines a diamond lattice with unit length 3.57A region carbonatoms block -96 -79 23 34 -8 11 # Defines a 3D block region called "abrasive" that is made up of the unit cell in the x-, y-, and z-direction for the given dimensions #create_box 1 box # Creates a simulation box for the abrasive region. #create_atoms 1 box # Creates atoms within the simulation box. #mass 1 12.0107 # Assign the mass of carbon. group carbonatoms region carbonatoms # Assigns the name carbonatoms to atom type 1. 3 atoms in group carbonatoms # Group the SARS data file region sars block -100 100 -26 185 -65 65 units box # Create a region for the datafile group sars region sars # group the input file with the name "sars" 51227 atoms in group sars # Interatomic potentials pair_style reax/c NULL # Pair potential style ReaxFF pair_coeff * * ffield.reax.FC C N O S # Assign Respective atoms Reading potential file ffield.reax.FC with DATE: 2013-06-28 WARNING: Changed valency_val to valency_boc for X (../reaxc_ffield.cpp:315) # Settings compute peratom all pe/atom # Compute potential energy per atoms neighbor 2.0 bin # NEVER KNEW WHAT THIS IS !!!!!!!!! #neigh_modify delay 5 #------------""----------""-------- neigh_modify delay 20 every 1 check yes page 500000 one 50000 # Initialization velocity all create 350 123456 # Setting temperature to 350 K variable t_step equal "step" # Assigning a variable for step variable t_temp equal "temp" # Assigning a variable for temperature thermo 100 # Show [#] for every 100 steps thermo_style custom step press temp pe # Show [temp and step & PE ] # Relaxation #fix rigid sars addforce 0.0 0.0 0.0 # Make the sars group immobile by reducing force to 'zero' in all direction #fix relax carbonatoms npt temp 250 250 0.5 iso 0 0 0.5 drag 1 # The start and end temperatures of abrasive is 250K and start and end pressure of abrasive is 0 and the abrasive is damped. fix charge all qeq/reax 1 0.0 10.0 1.0e-6 reax/c # Fix needed for reaxFF. DONNO WHAT THIS DOES !!!!!! variable t equal 1 # Assigning timestep timestep${t}                                                        # Equating timestep
timestep        1
dump            1 all custom 200 equil.*.dump id type x y z fx fy fz        # Dump 'atom id' 'atom type' 'x,y,z coordinates' 'force along x,y,z coordinates'
run             1                                                        # Run for 2000 time steps
Neighbor list info ...
update every 1 steps, delay 20 steps, check yes
max neighbors/atom: 50000, page size: 500000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 34 35 22
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair reax/c, perpetual
attributes: half, newton off, ghost
pair build: half/bin/newtoff/ghost
stencil: half/ghost/bin/3d/newtoff
bin: standard
(2) fix qeq/reax, perpetual, copy from (1)
attributes: half, newton off, ghost
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 1299 | 1299 | 1299 Mbytes
Step Press Temp PotEng
0   -10766.276          350   -3905465.6
1    -10766.28          350   -3905465.6
Loop time of 5.4404 on 1 procs for 1 steps with 51227 atoms

Performance: 0.016 ns/day, 1511.221 hours/ns, 0.184 timesteps/s

Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 3.7494     | 3.7494     | 3.7494     |   0.0 | 68.92
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 0.0001034  | 0.0001034  | 0.0001034  |   0.0 |  0.00
Output  | 0.0027405  | 0.0027405  | 0.0027405  |   0.0 |  0.05
Modify  | 1.688      | 1.688      | 1.688      |   0.0 | 31.03
Other   |            | 0.0001042  |            |       |  0.00

Nlocal:    51227 ave 51227 max 51227 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    13650 ave 13650 max 13650 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    6.23875e+06 ave 6.23875e+06 max 6.23875e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 6238751
Ave neighs/atom = 121.786
Neighbor list builds = 0
Dangerous builds = 0

# always remove fixes that are no longer needed
#unfix rigid
#unfix relax
#unfix charge
#Scratching

fix  2 carbonatoms move linear -5 0 0 units box                              # Apply force in x direction for wear @ 5 Angstroms/femtosecond
dump 3 all custom 100 sars-Scratch.*.dump id type x y z fx fy fz            # Dump 'atom id' 'atom type' 'x,y,z coordinates' 'force along x,y,z coordinates'
run 1000                                                                   # Run for 2000 time steps
Per MPI rank memory allocation (min/avg/max) = 1300 | 1300 | 1300 Mbytes
Step Press Temp PotEng
1   -10766.277          350   -3905465.6
100   -562.05515    8180.9385   -3905397.2
200   -562.05562    8180.9385   -3905397.2
300   -562.05498    8180.9385   -3905397.2
400   -562.05579    8180.9385   -3905397.2
500   -562.05559    8180.9385   -3905397.2
600    -562.0559    8180.9385   -3905397.2
700    -562.0552    8180.9385   -3905397.2
800   -562.05578    8180.9385   -3905397.2
900   -562.05555    8180.9385   -3905397.2
1000   -562.05588    8180.9385   -3905397.2
1001   -562.19269    8180.9385     -3905428
Loop time of 7549.05 on 1 procs for 1000 steps with 51227 atoms

Performance: 0.011 ns/day, 2096.958 hours/ns, 0.132 timesteps/s

Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 4092.2     | 4092.2     | 4092.2     |   0.0 | 54.21
Neigh   | 10.099     | 10.099     | 10.099     |   0.0 |  0.13
Comm    | 0.1589     | 0.1589     | 0.1589     |   0.0 |  0.00
Output  | 2721.7     | 2721.7     | 2721.7     |   0.0 | 36.05
Modify  | 724.69     | 724.69     | 724.69     |   0.0 |  9.60
Other   |            | 0.1211     |            |       |  0.00

Nlocal:    51227 ave 51227 max 51227 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    13650 ave 13650 max 13650 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    6.23872e+06 ave 6.23872e+06 max 6.23872e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 6238720
Ave neighs/atom = 121.786
Neighbor list builds = 50
Dangerous builds = 50

# End simulation
print "All done"                                                            # HOPE TO SEE THIS AT THE END OF SCREEN
All done

Please see the log.cite file for references relevant to this simulation

Total wall time: 2:06:06


-The above is the result after the simulation is done.
-The single atom(abrasive) which is supposed to move in x direction according to code does not move at all. Instead some other 2-3 atoms start moving outside the simulation box.
I don't understand what is happening with this.

3.Appendix
You can find all the necessary file for this issue over here

Long story short
-I expected the single abrasive atom to to move along the x-direction and do something to the S protein and M protein junction. It didn't. So if anyone can help me out regarding the same it would be amazing!

• I am so sorry for the extremely long question but it needed to be properly described in order to get a feasible solution! Jul 23 '20 at 4:28
• Welcome aboard! It may take a little longer to get an answer, since it will be more specific, but nonetheless it is most helpful to have the background info. :D Jul 24 '20 at 7:35
• @Etienne Palos, Haha I expected so. Jul 24 '20 at 7:36

Not a full answer, and it's too big for as a comment, but I can explain the linefix charge all qeq/reax 1 0.0 10.0 1.0e-6 reax/c.

This function uses the charge equilibration method (QEq) from Rappe and Goddard (https://doi.org/10.1021/j100161a070) to calculate the partial charges on each atom in the simulation. It's not moving anything in the simulation and in this case I don't think the partial charges are even feeding into anything so I reckon you can remove that line without impacting the model. If you're using a forcefield that needs partial charges then I normally use something like the following code to print the the partial charges:

group       type1 type 1    # For atomtype 1
compute     charge1 type1 property/atom q
compute     q1 type1 reduce ave c_charge1


Followed later by

thermo_style    custom step pe c_q1 ...    # Where ... is any other computes etc.


Hope that helps explain at least some of the code!

• +1 Beautiful. There's a reason multiple answers (from various users) are aloud, so no worries about the partial answer. Jul 24 '20 at 10:56
• Thanks that cleared a lot of doubts but if i remove the fix charge all qeq/reax 1 0.0 10.0 1.0e-6 reax/c command there will be an error. Jul 25 '20 at 3:50
• It seems if i use pair_style reax/c it is mandatory to use fix charge all qeq/reax 1 0.0 10.0 1.0e-6 reax/c althoight he numbers can vary~ Jul 25 '20 at 3:51
• @Matt any idea on how/why the individual atom is not moving? Jul 25 '20 at 3:52
• @SaiSmaranSBPES1201701189PES That's interesting to know, I didn't know the fix actually applied the charge to the atoms as it calculates them. I'm not sure why atom is not moving I'm afraid. There is no clear error so I would keep playing with the function till it works or exits on error.
– Matt
Jul 25 '20 at 13:11