Proteins+MARTINI+HOOMD

784 views
Skip to first unread message

Martin

unread,
Jun 27, 2012, 10:24:54 AM6/27/12
to hoomd...@googlegroups.com
Hi Hoomd-Blue community,

I've been hard at work implementing and testing the MARTINI scheme in HOOMD. It works very well for lipids and surfactants simulations. Lipid bilayer membranes (LBMs) of tens of nanometers in side length can be simulated over microseconds in just of a few days on Tesla 2070s and GTX580s. Typical example: for a system of size Lx = Ly = Lz ~= 20nm over ~ 1 microsec it takes around 24h. Speed could be further improved with an anisotropic version of the NPT barostat/thermostat since at the moment the NPH integrator is used to maintain pressure while non-conservative DPD interactions are used to maintain temperature which slows down the simulation. An anisotropic barostat is necessary in membrane simulations.

Irving-Kirkwood stress calculation for slices in the system has been implemented in HOOMD and should be generalized for volume elements shortly. It is done on the CPU at the moment as it is an infrequent calculation. A GPU version could be written but I'm not there yet.

Now we arrive to the main topic: simulating proteins within the MARTINI scheme. The problem here is that they use constraints (LINCS algorithm) to maintain some bond distances fixed which is easily feasible with rigid bodies but these constraints can be used within the anisotropic NPT integration scheme in GROMACS which is not the case at the moment for rigid bodies in HOOMD. I was wondering if you guys would know of a hackish way around this or if the NPH thermostat really needs to be adapted to rigid bodies. Could one simulate the rest of the system in NPH and the rigid bodies in NVE without ending up with some kind of explosion !? Some things would of course be violated, but the effect could potentially be mitigated no ? Just asking so that my work efforts can be focused in the appropriate direction.

My implementation of the MARTINI scheme is not the cleanest, but I would definitely be ready to share it with the community if anyone is interested.

Martin

Carolyn Phillips

unread,
Jun 27, 2012, 11:55:36 AM6/27/12
to hoomd...@googlegroups.com
Sorry Martin,
I got a little confused as to what you are precisely asking about.
Can you clarify?

You said,
>The problem here is that [The MARTINI scheme] use[s] constraints (LINCS algorithm) to
> maintain some bond distances fixed which is easily feasible with rigid
> bodies but these constraints can be used within the anisotropic NPT
> integration scheme in GROMACS which is not the case at the moment for rigid
> bodies in HOOMD.

What is not the case in HOOMD? An anisotropic NPT integration scheme
that works with rigid bodies?

I was wondering if you guys would know of a hackish way
> around this or if the NPH thermostat really needs to be adapted to rigid
> bodies.

Again, confused. HOOMD has npt and npt_rigid. Are you asking if
the HOOMD *did* have a anistropic NPT integrator, how would you then
hack it to work with rigid bodies?

Could one simulate the rest of the system in NPH and the rigid
> bodies in NVE without ending up with some kind of explosion !? Some things
> would of course be violated, but the effect could potentially be mitigated
> no ? Just asking so that my work efforts can be focused in the appropriate
> direction.


> --
> You received this message because you are subscribed to the Google Groups
> "hoomd-users" group.
> To view this discussion on the web visit
> https://groups.google.com/d/msg/hoomd-users/-/pIzVVNQHHroJ.
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to
> hoomd-users...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/hoomd-users?hl=en.

David LeBard

unread,
Jun 27, 2012, 12:40:34 PM6/27/12
to hoomd...@googlegroups.com
I think he means the LINCS constraint algorithm, similar to SHAKE, RATTLE, SETTLE, etc.  I would be very happy to see this implemented as well, because it would mean we could also run atomistic biomolecular force fields in HOOMD.  

However, I am not sure of a "hackish" way around constraints other than raising the bond force constant and dramatically lowering the simulation time step.  One unfortunate side effect to this hack could be the slowing down of the simulation throughput.  That said, we have done tests with the Shinoda-Devane-Klein CG model (CGCMM) in HOOMD where we lowered the timestep by a factor of 2, and found the simulation speed remained constant!  The reason was due to the reduced number of neighbor list rebuilds with a shorter time step.  So, it is possible you may be able to get away with this hack for the Martini protein model as well.   But of course constraints would be much cleaner and more useful in the long run...

Dave

Martin

unread,
Jun 28, 2012, 8:47:51 AM6/28/12
to hoomd...@googlegroups.com
Sorry if my post created confusion. Let us try again.

Proteins can be simulated within the MARTINI scheme using GROMACS. MARTINI proteins require that some bond distances be fixed which is done using the LINCS algorithm in GROMACS. If one wants to achieve this using HOOMD, I think that rigid bodies could definitely do the trick. So far so good.

Using GROMACS, one can simulate a MARTINI protein embedded in a lipid bilayer membrane. The membrane tension and the system pressure are kept constant using a semi-isotropic version of the NPT integrator that support LINCS. In HOOMD, a membrane alone can be simulated quite well using a combination of the NPH integrator and DPD interactions. However, once rigid constraints are introduced, one would need a version of the NPH integrator that supports rigid bodies which is not the case now. Hence the following two questions:

1) How hard would it be to write a rigid body compatible version of the NPH integrator. I have extensively modified HOOMD for my purposes, but comparing TwoStepNPT with TwoStepNPT_rigid makes me dizzy. A little guidance would be much appreciated !

2) Could one integrate regular particles in the NPH ensemble and the rigid bodies in the NVE ensemble. Assuming the time step is small enough and that the number of rigid body constraints is relatively this could work no ?

Thanks for helping yet again !

Martin


On Wednesday, June 27, 2012 11:55:36 AM UTC-4, Carolyn wrote:
Sorry Martin,
I got a little confused as to what you are precisely asking about.
Can you clarify?

You said,
>The problem here is that [The MARTINI scheme] use[s] constraints (LINCS algorithm) to
> maintain some bond distances fixed which is easily feasible with rigid
> bodies but these constraints can be used within the anisotropic NPT
> integration scheme in GROMACS which is not the case at the moment for rigid
> bodies in HOOMD.

What is not the case in HOOMD?   An anisotropic NPT integration scheme
that works with rigid bodies?

 I was wondering if you guys would know of a hackish way
> around this or if the NPH thermostat really needs to be adapted to rigid
> bodies.

Again, confused.   HOOMD has npt and npt_rigid.   Are you asking if
the HOOMD *did* have a anistropic NPT integrator, how would you then
hack it to work with rigid bodies?

Could one simulate the rest of the system in NPH and the rigid
> bodies in NVE without ending up with some kind of explosion !? Some things
> would of course be violated, but the effect could potentially be mitigated
> no ? Just asking so that my work efforts can be focused in the appropriate
> direction.


> hoomd-users+unsubscribe@googlegroups.com.

Carolyn Phillips

unread,
Jun 28, 2012, 10:37:07 AM6/28/12
to hoomd...@googlegroups.com, ndt...@umich.edu
Hi Martin,

I see what you are after now!

The expert on implementing new integrators for rigid bodies is Trung
Nguyen, who has not weighed in here yet.

I am going to make the quick observation that yes I expect it can be
done because in Lammps, there is a fix that applies NPH applied to
extended spherical, and the documentation notes that the position,
velocity, AND angular velocity is updated in this case. What that
indicates to me is that the technical questions of "how do you apply a
NPH integrator to a non-point body" has been worked out and is
probably not that difficult.

As for mixing the two systems. I *think* that if the number of
things that were not being integrated by the NPH integrator were
small, that they would simply make the rate at which the system as a
whole responded to changes in the parameters slower (e.g. the coupling
would effectively be weaker than you set it). However, I am assuming
that the changing box size would not somehow cause non-NPH integrated
particles to end of "out of the box".

Alternatively, you can puzzle out how to make an NPH_rigid yourself
(and I am sympathetic to the difficulties of that!) or hope that a
developer familiar with the rigid parts of the code can find some time
to put it in.

-Carolyn
>> > hoomd-users...@googlegroups.com.
>> > For more options, visit this group at
>> > http://groups.google.com/group/hoomd-users?hl=en.
>
> --
> You received this message because you are subscribed to the Google Groups
> "hoomd-users" group.
> To view this discussion on the web visit
> https://groups.google.com/d/msg/hoomd-users/-/WrE8ckbbtZUJ.
>
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to
> hoomd-users...@googlegroups.com.

Trung Nguyen

unread,
Jun 28, 2012, 10:57:59 AM6/28/12
to hoomd...@googlegroups.com
Hi Martin,

I think implementing an NPH integrator for rigid bodies is doable but
it seems to be much more messy for anisotropic pressure tensors than
for the isotropic case. And I don't find any time to work on it in the
near future. Because the NPH rigid integrator is supposed to resize
the box (like NPT rigid), using it along with another NPH integrator
would invalidate the other (and vice versa)- which I suppose you don't
want. Also, if you apply the NPH rigid integrator for a small number
of atoms (in the rigid body group), the fluctuation in pressure would
be large and the simulation is more likely to crash.

I suppose using NVE rigid coupled with an NPH integrator for a sea of
non-rigid particles would be fine, the NVE rigid does nothing more
than moving the bodies according to the Newton's second law (despite
the fact that V is there). Similar to BoxResizeUpdater, the non-rigid
NPH integrator would rescale the bodies' center of mass (if any
present) after resizing the box dims. As long as you can sample a
Boltzmann distribution of V from the mixed integrators, I think the
system is reasonably in an NPH ensemble.

Cheers,
-Trung

Joshua A. Anderson

unread,
Jul 2, 2012, 8:17:42 AM7/2/12
to hoomd...@googlegroups.com
- Constraint dynamics: We have plans to implement SHAKE and/or LINCS and/or (insert bond constraint algorithm here) in hoomd. I do not have an ETA on completion at this time.

- Anisotropic NPT (for point particles): This has been requested so many times…. The author of the NPT code has offered to update it. I also do not have an ETA on this.
- Anisotropic NPT (for rigid bodies): Is there even a known derived method to do this?

Going only slightly hackish, a Berendesen pressure updater would be relatively easy to implement. LAMMPS has an anisotropic one, so that appears to be solved. If I understand how it works, it should function with any combination of particle and/or rigid body integrators. And maybe it is not that hackish - doesn't MARTINI specify that a berendsen npt scheme is to be used?

The only barrier to implementing a Berendsen pressure updater is to implement the box scale routine on the GPU. update.box_resize currently works on the CPU. I guess this code is also replicated in the various npt schemes - perhaps we should unify them into one utility class that can be used by all.
--------
Joshua A. Anderson, Ph.D.
Chemical Engineering Department, University of Michigan

Rastko Sknepnek

unread,
Jul 2, 2012, 9:00:32 AM7/2/12
to hoomd...@googlegroups.com
I am moving to a new position this fall, which means that I need to finish a number of projects none of which are MD related. In short, unfortunately I won't have time to look into the NPT code for at least next five to six months.

I have an extensive set of notes on the anisotropic (still orthogonal box, thou) NPT integrator. If anyone is interested in implementing it, I could share my notes and discuss what needs to be done. The implementation should be straightforward.

Rastko

-----------------------------------------------
Rastko Sknepnek
Research associate
Department of Materials Science
Northwestern University
(847) 467 2130 (office)
(515) 520 1598 (cell)

Jens Glaser

unread,
Jul 2, 2012, 9:51:51 AM7/2/12
to hoomd...@googlegroups.com
Hi all,

actually I started work on the anisotropic NPT integrator. Thus far, I have confirmed that the equations of motion of

Yu, T.-Q., Alejandre, J., López-Rendón, R., Martyna, G. J., & Tuckerman, M. E. (2010). Measure-preserving integrators for molecular dynamics in the isothermal–isobaric ensemble derived from the Liouville operator. Chemical Physics, 370(1-3), 294-305. Elsevier B.V. doi:10.1016/j.chemphys.2010.02.014

give the right ensemble, when specialized to orthorhombic and tetragonal integration modes. I plan on implementing the integrator this week (finally),
and let you know of any progress. I will also analyze HOOMD's existing NPT integrator which you wrote. Should I have any questions, I will contact you.

Best,
Jens.

Martin

unread,
Jul 3, 2012, 10:06:57 AM7/3/12
to hoomd...@googlegroups.com
It's true that the foundation MARTINI papers refer to a Berendsen NPT integration scheme. It might be the best solution to the problem initially posed in this post given that it shouldn't be too hard to implement. I will take a look at it and see if I can copy some of the work done in the LAMMPS code.

M
>>>>> For more options, visit this group at
>>>>> http://groups.google.com/group/hoomd-users?hl=en.
>>>
>>> --
>>> You received this message because you are subscribed to the Google Groups
>>> "hoomd-users" group.
>>> To view this discussion on the web visit
>>> https://groups.google.com/d/msg/hoomd-users/-/WrE8ckbbtZUJ.
>>>
>>> To post to this group, send email to hoomd...@googlegroups.com.
>>> To unsubscribe from this group, send email to
>>> For more options, visit this group at
>>> http://groups.google.com/group/hoomd-users?hl=en.
>>
>> --
>> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
>> To post to this group, send email to hoomd...@googlegroups.com.
>> To unsubscribe from this group, send email to hoomd-users+unsubscribe@googlegroups.com.
>> For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.
>>
>
> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to hoomd-users+unsubscribe@googlegroups.com.

Joshua A. Anderson

unread,
Jul 3, 2012, 10:11:11 AM7/3/12
to hoomd...@googlegroups.com
FYI - I can't accept any patches that copy even a single line of code from LAMMPS. LAMMPS is licensed under the GPL.
--------
Joshua A. Anderson, Ph.D.
Chemical Engineering Department, University of Michigan

> >>>>> hoomd-users...@googlegroups.com.
> >>>>> For more options, visit this group at
> >>>>> http://groups.google.com/group/hoomd-users?hl=en.
> >>>
> >>> --
> >>> You received this message because you are subscribed to the Google Groups
> >>> "hoomd-users" group.
> >>> To view this discussion on the web visit
> >>> https://groups.google.com/d/msg/hoomd-users/-/WrE8ckbbtZUJ.
> >>>
> >>> To post to this group, send email to hoomd...@googlegroups.com.
> >>> To unsubscribe from this group, send email to
> >>> hoomd-users...@googlegroups.com.
> >>> For more options, visit this group at
> >>> http://groups.google.com/group/hoomd-users?hl=en.
> >>
> >> --
> >> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> >> To post to this group, send email to hoomd...@googlegroups.com.
> >> To unsubscribe from this group, send email to hoomd-users...@googlegroups.com.
> >> For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.
> >>
> >
> > --
> > You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> > To post to this group, send email to hoomd...@googlegroups.com.
> > To unsubscribe from this group, send email to hoomd-users...@googlegroups.com.
> > For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.
> >
> >
> >
>
>
> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To view this discussion on the web visit https://groups.google.com/d/msg/hoomd-users/-/oPktTkKKb34J.
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to hoomd-users...@googlegroups.com.

Martin

unread,
Jul 3, 2012, 10:41:08 AM7/3/12
to hoomd...@googlegroups.com
No problem, no code will be directly copied. In fact the Berendsen scheme is so simple to implement that code need not be copied.

M
> >>>>> For more options, visit this group at
> >>>>> http://groups.google.com/group/hoomd-users?hl=en.
> >>>
> >>> --
> >>> You received this message because you are subscribed to the Google Groups
> >>> "hoomd-users" group.
> >>> To view this discussion on the web visit
> >>> https://groups.google.com/d/msg/hoomd-users/-/WrE8ckbbtZUJ.
> >>>
> >>> To post to this group, send email to hoomd...@googlegroups.com.
> >>> To unsubscribe from this group, send email to
> >>> For more options, visit this group at
> >>> http://groups.google.com/group/hoomd-users?hl=en.
> >>
> >> --
> >> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> >> To post to this group, send email to hoomd...@googlegroups.com.
> >> To unsubscribe from this group, send email to hoomd-users+unsubscribe@googlegroups.com.
> >> For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.
> >>
> >
> > --
> > You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> > To post to this group, send email to hoomd...@googlegroups.com.
> > To unsubscribe from this group, send email to hoomd-users+unsubscribe@googlegroups.com.
> > For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.
> >
> >
> >
>
>
> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To view this discussion on the web visit https://groups.google.com/d/msg/hoomd-users/-/oPktTkKKb34J.
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to hoomd-users+unsubscribe@googlegroups.com.

Martin

unread,
Jul 5, 2012, 2:05:37 PM7/5/12
to hoomd...@googlegroups.com
An update on the Berendsen thermostat and anisotropic barostat supporting rigid bodies. 

After a bit of reading I decided that for the moment, the best I could do was to implement simple updaters for both and integrate the system forward using nve and nve_rigid. I therefore wrote a plugin that contains both the thermostat and an anistropic version of the barostat that supports three modes much like the NPH integrator written by Jens: isotropic, anisotropic, semi-isotropic.

All tests done so far show that it works quite well. To test the semi-isotropic version of the barostat I simulated a MARTINI DPPC bilayer and obtained results equal to those previously obtained. I also simulated a square colloid (rigid body) suspended in a LJ fluid and found that both temperature and pressure were maintained.

The updaters are in a plugin at the moment, and a little rough on the edges. They could probably be combined in one to optimize the rigid body parts especially the calls to set_RV. This will be done at a later moment as some serious testing is about to be undergone. However, if any of you is interested in using what I have written (with no guarantees whatsoever !) I will definitely share.

M

On Monday, July 2, 2012 8:17:42 AM UTC-4, Joshua A. Anderson wrote:
>>>>> For more options, visit this group at
>>>>> http://groups.google.com/group/hoomd-users?hl=en.
>>>
>>> --
>>> You received this message because you are subscribed to the Google Groups
>>> "hoomd-users" group.
>>> To view this discussion on the web visit
>>> https://groups.google.com/d/msg/hoomd-users/-/WrE8ckbbtZUJ.
>>>
>>> To post to this group, send email to hoomd...@googlegroups.com.
>>> To unsubscribe from this group, send email to
>>> For more options, visit this group at
>>> http://groups.google.com/group/hoomd-users?hl=en.
>>
>> --
>> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
>> To post to this group, send email to hoomd...@googlegroups.com.
>> To unsubscribe from this group, send email to hoomd-users+unsubscribe@googlegroups.com.
>> For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.
>>
>
> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to hoomd-users+unsubscribe@googlegroups.com.

Jens Glaser

unread,
Jul 5, 2012, 2:16:56 PM7/5/12
to hoomd...@googlegroups.com
Martin,

did you check that your equations of motion give the correct ensemble? Is the integrator symmetric and measure-preserving? Is the conserved quantity really conserved, i.e. fluctuating around its initial value,
or is there a drift?

The change from cubic to orthorhombic or tetragonal is in principle not hard, but it is necessary to ascertain oneself of the above points before writing or modifying an integrator.
This said, I am working on the anisotropic NPT, too, but we can discuss in order not to duplicate efforts.

Jens

To view this discussion on the web visit https://groups.google.com/d/msg/hoomd-users/-/pgrOL-CeRT4J.

To post to this group, send email to hoomd...@googlegroups.com.
To unsubscribe from this group, send email to hoomd-users...@googlegroups.com.

Martin

unread,
Jul 5, 2012, 2:40:00 PM7/5/12
to hoomd...@googlegroups.com
Hi Jens,

As I said, pretty much all properties of the systems simulated using the NPH integrator were preserved within uncertainty margins. Pressure and temperature appear to be truly conserved, that is they fluctuate around the target values. I still haven't tested the anisotropic or semi-isotropic barostat with rigid bodies. This will be done shortly.

Again, I did not write a new integrator. I proceeded in much the same way LAMMPS proceeds, that is the Berendsen thermostat and barostat are updaters activated alongside the nve and nve_rigid integrators.

It is still debated in the literature whether the Berendsen scheme is physical or not. But it sure appears that in many cases, it is more than sufficient. As I said, I will share the code with anybody that is interested, and I'd be pleased to collaborate with you to reduce the overall workload.

M

Jens

unread,
Jul 12, 2012, 5:26:22 PM7/12/12
to hoomd...@googlegroups.com
Hi all,

I implemented an anisotropic NPT scheme that is designed to (rigorously!) sample the NPT ensemble. It supports cubic, orthorhombic and tetragonal integration modes,
and is based on the Martyna-Tobias-Klein barostat/thermostat.

Those who are interested can try out the patch here
(only for the brave!), or wait until it is hopefully implemented in HOOMD-blue.

Best,
Jens.

GPU-Greybeard

unread,
Nov 20, 2012, 4:53:35 AM11/20/12
to hoomd...@googlegroups.com
On Wednesday, June 27, 2012 4:24:54 PM UTC+2, Martin wrote:
Hi Hoomd-Blue community,

I've been hard at work implementing and testing the MARTINI scheme in HOOMD. It works very well for lipids and surfactants simulations.
...

Hi Martin,

how did you solve the G96Angle-problem in HOOMD? As you know, MARTINI uses a GROMACS's cosine-based angle term (aka #2) in contrast to the quadratic potential available in HOOMD.

Thanks & regards

GB

Martin Bertrand

unread,
Nov 20, 2012, 2:44:20 PM11/20/12
to hoomd...@googlegroups.com
Hi GB,

I wrote a martini plugin which contains this COSINE2 potential that I just posted online. Here's the message I sent to another inquirer:

***************************************************************************************
So I finally have posted my MARTINI plugin online on GitHUB !


 At the moment it will only work with HOOMD-Blue v0.10 NOT v0.11. The small issue is with BoxDim which I should correct soon enough.

All you need to do is:
- make sure v0.10 is the active version
- download plugin
- make a directory to build the plugin and cd to it
- cmake <plugin_directory>
- in your simulation script "from hoomd_plugins import martini_plugin asmartini"
- to initialize your system provided the proper xml file you run the following command "system,pairlj,paircoul,harm,angle = martini.utils.quick_setup('system.xml','martini_v2.0.itp')" where 'system.xml' is the name of your xml file !

That's it, once installed, one line and the utility takes care of everything ... or almost. Dihedrals are not set up at the moment.
***************************************************************************************

You can download the code and check the cosine2 potential in there. As I say in the above message, it only works with v0.10 of HOOMD but the changes for v0.11 are minimal and should be done at a later time.

If you just want to use the potential you can run the following commands:
from hoomd_plugins import martini_plugin as martini
cos2 = martini.angle.cos2()
cos2.set_coeff(<type>,k=<spring constant>,t0=<eq. angle>)

Any other questions, feel free to ask !

M


--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To view this discussion on the web visit https://groups.google.com/d/msg/hoomd-users/-/zf3IKIDSigcJ.

To post to this group, send email to hoomd...@googlegroups.com.
To unsubscribe from this group, send email to hoomd-users...@googlegroups.com.

For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.



--
Martin Bertrand, Ph.D.
Department of Physics
University of Ottawa


Reply all
Reply to author
Forward
0 new messages