Why Protein molecule not moving after giving velocity

Hi…

I am writing Lammps script for measuring friction coefficient between two protein molecules. For that, I am giving velocity to the top protein molecule and not giving velocity to the bottom protein molecule. But to my surprise, the top protein molecule is not moving. It is only vibrating at its own plave. Can someone help me out with this?

I am posting here my lammps input script.

#------Initialize Simulation------
units metal
dimension 3
boundary p p p

atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic

#------Define Interatomic Potential------
pair_style lj/charmm/coul/long 8 10
pair_modify mix arithmetic
kspace_style pppm 1.0e-4
neighbor 0.3 bin
neigh_modify delay 2 every 1

#------Reading LAMMPS data_file------
read_data data.assembly
read_data data.assembly add append shift 0.0 12.0 0.0

pair_coeff * * 100.0 2.0
special_bonds charmm

group collagen id <= 1275
group lubricin id > 1275
group mobile union collagen lubricin
group rest subtract all mobile

#------Energy Minimization------
reset_timestep 0
timestep 0.001
thermo 100
thermo_style custom step etotal emol press temp lx ly lz
min_style cg
minimize 1e-25 1e-25 100 500
run 200

#------Equilibration------ Physical propeties become consant after #equilibration for the rest of the simulation

#velocity all create 300.0 1231 mom yes rot yes
#fix 1 collagen nve/limit 1.0

fix 3 lubricin rigid/nve single force * on on on torque * off off off langevin 300.0 300.0 10.0 904297

#fix 2 lubricin langevin 300.0 300.0 10.0 904297
#fix 4 collagen setforce 0.0 0.0 0.0

velocity lubricin set 0.0 0.0 1.0

#velocity collagen set 0.0 0.0 0.0

reset_timestep 0
thermo 500
timestep 0.001
thermo_style custom step temp press
dump 1 all custom 50 trajectory.lammpstrj id x y z ix iy iz

run 3000

SIMULATION DONE

print “All done”

Hi..

I am writing Lammps script for measuring friction coefficient between two
protein molecules. For that, I am giving velocity to the top protein

not quite clear to me, what you mean by "friction" in this context.
also, the potential parameters of

pair_coeff * * 100.0 2.0

seem *very* unusual for a protein. this is a despite using CHARMM
force field styles throughout, you are not using CHARMM parameters,
right?
it is also very unusual to have all particles with the same
parameters. this is a very deep and thus quite narrow potential shape.

molecule and not giving velocity to the bottom protein molecule. But to my
surprise, the top protein molecule is not moving. It is only vibrating at
its own plave. Can someone help me out with this?

difficult to say without knowing more details about your system, e.g.
whether there is a solvent or how many particles of what kind are
around beyond the protein.

one thing that stands out, besides the unusual parameters is the use
of a langevin thermostat for translation. that one is *supposed* to
dissipate translation.
hence the friction that you would actually see, would be primarily
that of the langevin thermostat (or rather the viscosity of the
implicit solvent it would represent), and when reducing the whole
protein to just 3 translational degrees of freedom, it does not
surprise me, that you get unusual results.

axel.

Thanks axel.

“not quite clear to me, what you mean by “friction” in this context.”

as you can see I have used read_data commands twice. so, they are creating two beds one above the other as I have used shift command after read_data.

read_data data.assembly
read_data data.assembly add append shift 0.0 12.0 0.0

now, i am giving velocity to the top bed which is group lubricin containing atoms with id greater than 1275 in the combined data file. actually, in data.assembly file there are 1275 atoms. after that, I am finding the shear force and normal force on the atoms in the lubricin group. the (friction coefficient = shear force/normal force).
i have not shown the complete input script for calculating the coefficient of friction above because the unexpected result( top protein bed not moving) is coming well before I am finding the forces.

“difficult to say without knowing more details about your system, e.g.
whether there is a solvent or how many particles of what kind are
around beyond the protein.”

there are no other solvents or other particles in my system beyond the two sets of proteins.

" pair_coeff * * 100.0 2.0

seem very unusual for a protein. this is a despite using CHARMM
force field styles throughout, you are not using CHARMM parameters,
right?
it is also very unusual to have all particles with the same
parameters. this is a very deep and thus quite narrow potential shape."

even if, I don’t use pair_coeff still there is no change in the result. what parameters are you suggesting me to change for different particles?

“one thing that stands out, besides the unusual parameters is the use
of a langevin thermostat for translation. that one is supposed to
dissipate translation.”

instead of langevin thermostat, if i use temp/rescale thermostat still the situation not changes.

I am attaching photo of my lammps trajectory file. I can’t attach the dynamic view in which top protein bed is not moving translationally but only vibrating at its place.
PL help me out with this.

protein-bed.PNG

i cannot "fix" your input. many details sound as if you have very
little experience with MD. the same goes for your input parameters and
settings. they all don't make much sense to me. this all looks as if
you have just assembled it with cut-n-paste from various pieces
without exactly knowing what you are doing. there is not much to
infer from your visualization apart from that it looks like you
*should* have a solvent here. the fact that you don't worry about
wiping out force field settings is even more unsettling.
finally, i also still don't get what you are trying to compute here.

how about trying a simpler model system first? something where you
have reference data from a publication that you can reproduce and that
demonstrates that the method you want to use is actually working. and
where you can prove to yourself, that you are using it correctly?

axel.

thanks for your reply Axel. my script has worked and I got the desired result.

you are right , I am new to LAMMPS and this was my first script.

I have few more doubts to ask. I got the following error “out of range atoms” during the execution of my program. When I searched for the possible explanation for this, i found following:

Error: “Out of range atoms - cannot compute PPPM”

Explanation: One or more atoms are attempting to map their charge to a PPPM grid point that is not owned by a processor. This is likely for one of two reasons, both of them bad. First, it may mean that an atom near the boundary of a processor’s sub-domain has moved more than 1/2 the neighbor skin distance without neighbor lists being rebuilt and atoms being migrated to new processors.This also means you may be missing pairwise interactions that need to be computed. The solution is to change the re-neighboring criteria via the neigh_modify command. The safest settings are “delay 0 every 1 check yes”. Second, it may mean that an atom has moved far outside a processor’s sub-domain or even the entire simulation box. This indicates bad physics, e.g. due to highly overlapping atoms, too large a timestep, etc.

The explanation is not quite clear to me.

  1. What does first line mean “One or more atoms are attempting to map their charge to a PPPM grid point that is not owned by a processor”?

  2. What does it mean “you may be missing pairwise interactions”?

Can you pl elaborate more on this?

thanks for your reply Axel. my script has worked and I got the desired
result.
you are right , I am new to LAMMPS and this was my first script.

I have few more doubts to ask. I got the following error "out of range
atoms" during the execution of my program. When I searched for the possible
explanation for this, i found following:

Error: "Out of range atoms - cannot compute PPPM"

Explanation: One or more atoms are attempting to map their charge to a PPPM
grid point that is not owned by a processor. This is likely for one of two
reasons, both of them bad. First, it may mean that an atom near the boundary
of a processor's sub-domain has moved more than 1/2 the neighbor skin
distance without neighbor lists being rebuilt and atoms being migrated to
new processors.This also means you may be missing pairwise interactions that
need to be computed. The solution is to change the re-neighboring criteria
via the neigh_modify command. The safest settings are "delay 0 every 1 check
yes". Second, it may mean that an atom has moved far outside a processor's
sub-domain or even the entire simulation box. This indicates bad physics,
e.g. due to highly overlapping atoms, too large a timestep, etc.

The explanation is not quite clear to me.

1. What does first line mean "One or more atoms are attempting to map their
charge to a PPPM grid point that is not owned by a processor"?

LAMMPS does domain decomposition. this means, each processor "owns"
some part of the total system, but it also has a "skin" of atoms that
are copies from neighboring domains. the former atoms are referred to
a "local" atoms, the latter as "ghost" atoms.

LAMMPS - like most MD codes - uses neighbor lists. this means, rather
than looping over all pairs of atoms (which is an O(N**2) process), it
maintains lists of possible neighbors (which will result in O(N)
performance). this results in a significant speedup, especially for
systems with many atoms. but if your atoms move around, these lists
needs to be regenerated and there commands in LAMMPS to tweak them.

when you get this error it almost always means, that your input is
bad. either your neighbor list settings are bad, but more likely your
force field parameters or initial geometry are bad and thus atoms are
"catapulted" through your system.

2. What does it mean "you may be missing pairwise interactions"?

see above. when your neighbor lists are out-of-date, e.g. due to fast
moving atoms, then there are atoms within the cutoff range that are
not considered, which in turn means, your simulation is bogus.

Can you pl elaborate more on this?

these topics have been discussed on this mailing list many, many
times. your are encouraged to revisit old discussions in the mailing
list archives. but more importantly, you need to consult with a local
MD expert, that can give you sufficient guidance (these are not really
LAMMPS problems, but issues pertaining to all parallel MD codes). this
mailing list is not a classroom.

axel.

Thanks a lot Axel for such elaborative explanation.