[lammps-users] nvt integration for partially rigid molecules and minimization

I have a system composed by a drople bisferrocene bound to an alkyl group.
The cp rings have a dummy atom that binds to the fc. the two groups composed by the fc ring and the dummy atom are fixed with fix-rigid.
the bonds/angles/dihedrals to the central iron are left free to move.
I know I can’t use energy minimization properly because it doesn’t work well with fix rigid.
What I’m not understanding is:
1)How do I perform time integration with fix nvt without integrating twice but at the same time while considering the rigid body in the time integration
2)How should I initialize my system for simulated annealing? I have to perform some steps of energy minimization before dynamics, but ramping down with fix nvt from a temperature of 300K remains stuck at the first step of time integration with very high initial energies (maybe because some molecules end up too close) and so the system blows up in the first steps without minimization occurring

thanks in advance
Mattia Siviero

I have a system composed by a drople bisferrocene bound to an alkyl group.
The cp rings have a dummy atom that binds to the fc. the two groups composed by the fc ring and the dummy atom are fixed with fix-rigid.
the bonds/angles/dihedrals to the central iron are left free to move.
I know I can’t use energy minimization properly because it doesn’t work well with fix rigid.
What I’m not understanding is:
1)How do I perform time integration with fix nvt without integrating twice but at the same time while considering the rigid body in the time integration

you define two groups. one group has all atoms that are part of rigid bodies, the other group has all other atoms.
fix rigid (or its variants) will include time integration for the atoms that are in the rigid bodies, so you only need explicit time integration for the remaining atoms.

2)How should I initialize my system for simulated annealing? I have to perform some steps of energy minimization before dynamics, but ramping down with fix nvt from a temperature of 300K remains stuck at the first step of time integration with very high initial energies (maybe because some molecules end up too close) and so the system blows up in the first steps without minimization occurring

while you cannot do minimization for the rigid bodies, you can do minimization for all other atoms. just use fix setforce 0.0 0.0 0.0 on the rigid body atoms and they will remain in place while all other atoms can be minimized. if after this minimization you still have problems relaxing the system, you can do time integration with fix nve/limit on the non-rigid atom group and have no fix on the rigid bodies. that will (again) only move the atoms not included in rigid bodies. after that you can do time integration with fix rigid only and no time integration for anything else to only relax the rigid bodies somewhat. at that point you should have a reasonably relaxed system for starting dynamics for all. and then you should run with fix nve and fix rigid (without any thermostat) and figure out which length of timestep you need to have sufficient conservation of energy. rigid body integrators can be quite sensitive and require short time steps (for normal atomic/molecular systems they can be in the range 0.1- 0.5fs how large depends on the shape of the bodies and their moments of inertia. the more spherical the larger).

once this is working well, you can turn on thermostats and do whatever cooling you want. you can set an initial kinetic energy with the velocity command and then (re)create the rigid body fix and define your temperature ramp.

Axel.

Thanks for the answer. what I noticed was happening running damped dynamics with fix nvt on all the system was that the non rigid sidechains stretched well beyond what should happen in a harmonic bond where the potential profile is quadratic. So I ended up with heads that were relatively stationary and sidechains that stretched enormously. I find this especially weird when after extending, the alkyl sidechain don’t seem to relax to their equilibrium angle and bond values.

I’ll try to relax the sidechains before relaxing the whole system as you told me. I’m wondering about the role of fix nve in relaxing the system: if I’m imposing conservation of energy how am I relaxing the system?
thank you,
Mattia Siviero

I’m wondering about the role of fix nve in relaxing the system: if I’m imposing conservation of energy how am I relaxing the system?

fix nve implements plain velocity verlet timestepping with no other changes (e.g. thermostatting). it does not add or remove energy so you will get an NVE ensemble simulation, if multiple conditions are met:

  • there are no other system manipulations (changes to forces or velocities or positions)
  • you have a fully periodic system
  • the system is in equilibrium
  • the choice of simulation parameters is such that energy is conserved

in short, fix nve does not “enforce” energy conservation (or an NVE ensemble) but is a prerequisite for being able to achieve it.

axel.