Dear lammp users,
I’m currently working on a cholesterol model using martini force field. For molecules with ring structures like cholesterol, it requires constrained bonds to keep the rigidity of molecule during simulation. While Lammps has ’ fix shake’ command which uses ‘shake’ algorithm as ‘constraint’ function does in Gromacs, the two commands have different input values. ‘Constraint’ command fixes the distances between two atoms while ‘fix shake’ command inputs atom, mass, bond and angle types. I’m thinking using ‘bond’ command as an alternative. Since the bond type is harmonic, the alternative could be bond_style harmonic, bond_coeff ro K, where the force constant K is set to be a very large number (10^7 for instance). In this way the ‘spring’ is very rigid, and bond distance will not fluctuate. I’m not sure whether this is a good way.
Thank you,
Sean
Dear lammp users,
I'm currently working on a cholesterol model using martini force field. For
molecules with ring structures like cholesterol, it requires constrained
bonds to keep the rigidity of molecule during simulation. While Lammps has '
fix shake' command which uses 'shake' algorithm as 'constraint' function
does in Gromacs, the two commands have different input values. 'Constraint'
command fixes the distances between two atoms while 'fix shake' command
inputs atom, mass, bond and angle types. I'm thinking using 'bond' command
as an alternative. Since the bond type is harmonic, the alternative could be
bond_style harmonic, bond_coeff ro K, where the force constant K is set to
be a very large number (10^7 for instance). In this way the 'spring' is very
rigid, and bond distance will not fluctuate. I'm not sure whether this is a
good way.
the shake solver in LAMMPS does not attempt to handle clusters of
connected bonds, as that would negatively impact parallel efficiency.
to keep larger bodies rigid, you need more than just bonds.
using stiff bonds has the problem that the largest possible time step
is restricted, thus you'd have to use a multi-timestepping scheme like
r-RESPA to achieve suitable efficiency. on the other hand it is fairly
straightforward to implement and easy to monitor how effective it is.
an alternate approach is to use a rigid body integrator on the rigid
body atoms and a "normal" integrator on the remaining atoms.
axel.
Dear Axel,
Thank you for your advice. I tried to use run_style respa command to solve the timestep issue related with stiff bonds. It gives warning message during simulation:
run_style respa 6 2 2 2 2 2
Respa levels:
1 = bond angle dihedral improper
2 =
3 =
4 =
5 =
6 = pair kspace
WARNING: One or more respa levels compute no forces (…/respa.cpp:210)
I’m using default setting of run respa, and my system only contains bonds, angles and impropers. Is it because there is no dihedral that it gives warning?
Sean
Dear Axel,
Thank you for your advice. I tried to use run_style respa command to solve
the timestep issue related with stiff bonds. It gives warning message during
simulation:
run_style respa 6 2 2 2 2 2
Respa levels:
1 = bond angle dihedral improper
2 =
3 =
4 =
5 =
6 = pair kspace
WARNING: One or more respa levels compute no forces (../respa.cpp:210)
I'm using default setting of run respa, and my system only contains bonds,
angles and impropers. Is it because there is no dihedral that it gives
warning?
not it is because if the reason that the warning already tells you.
you have 6 respa levels, but you are computing forces on only two of
those levels. that is pointless. also, i doubt that the factor 32
difference in time step between the inner and outer timestep is what
you want. two respa levels should be sufficient for this. chances are,
you are misunderstanding some of the documentation of the respa
runstyle.
axel.
Axel,
Thank you for your reply. My system has a very small toleration on timestep. When timestep is larger than 1 fs, it gives bond atom missing error. Therefore I have to make sure the inner timestep where the bond force is calculated is smaller than 1 fs. The timestep corresponding with respa command is 5 fs. But since I have multi- timestep commands with the largest timestep=30 fs. I decide to use 6 levels. I also don’t want to change the timestep, because I want the pair potential to be calculated at same timestep (highest level). I’am pasting the relating commands.
timestep 0.1
run 5000
timestep 0.5
run 5000
timestep 1.0
run 10000
timestep 5.0
run_style respa 6 2 2 2 2 2
run 25000
timestep 10.0
run 50000
timestep 20.0
run 100000
timestep 30.0
run 500000
Axel,
Thank you for your reply. My system has a very small toleration on timestep.
When timestep is larger than 1 fs, it gives bond atom missing error.
Therefore I have to make sure the inner timestep where the bond force is
calculated is smaller than 1 fs. The timestep corresponding with respa
command is 5 fs. But since I have multi- timestep commands with the largest
timestep=30 fs. I decide to use 6 levels. I also don't want to change the
timestep, because I want the pair potential to be calculated at same
timestep (highest level). I'am pasting the relating commands.
as i already mentioned, you are misunderstanding the respa run style.
what you do:
run_style respa 6 2 2 2 2 2
is equivalent to:
run_style respa 2 32
axel.