Looking for development recommendations to speed up PLUMED calculation

[John Haberstroh]

Hello PLUMED users!

Not sure if this is the developers list as well as the users list, but I thought I'd give it a shot.

I love PLUMED and I would like to use it to replace a current GROMACS hack that I have going. I think it is a wonderful piece of software for the MD community, and I found modifying it delightfully simple.

My calculation at hand is a biased ensemble calculation of an electrostatic coupling of pseudo-charges within a small ligand to the rest of the system. When I say pseudo-charges, I mean charges that are only intended to be relevant to this calculation, and that should *not* be treated as real charges within the dynamics, especially since these pseudo-charges do not interact with one another. My system has about 100k "environment" atoms and 40 pseudo-charged atoms within the ligand.

I modified the DHENERGY class to accept an argument for PSEDUOCHARGE_A, a list of pseudocharges for GROUP_A. I then compute all of the distances and interactions with these charges to the system. I have a few problems, though:

1. The calculation is SLOW. It slows down gromacs by a factor of 10 for my hacked code, and a factor of 30 for the PLUMED version. As it should be: all of the other long-range coulomb couplings are taken care of with a nice FFT. However, due to the structure of this calculation, it does not parallelize for PLUMED nicely. I have three ideas for solutions.

1A. I could run this calculation in parallel, by communicating the pseudocharges and positions to every node, and then do the computation in parallel. I don't know how to do this within the PLUMED framework, but it seems like one of the easier fixes.

1B. I could implement biasing via Hybrid Monte Carlo. This would be excellent, and I was looking into hacking this into gromacs directly, but then I wondered if the PLUMED team has thought about implementing this. It could be a nice addition to the current functionality of plumed for order parameters that are too costly to compute every MD step.

1C. I could implement the coupling via some sort of PME modification for coupling two disjoint sets of atoms, which really sounds like a pain...

2. The PLUMED calculation is *wrong*. My system is in a non-rectilinear (it takes the compact representation of a rhombic dodecahedron), and I think that DHENERGY (and possibly all of PLUMED) is computing PBC with a rectilinear nearest image convention. I was doing the same with my home-written code and I was getting wildly unstable answers like I am currently getting with PLUMED. I was able to compare to the built-in GROMACS dist() function, and I tracked down the error. When you are using a non-rectilinear PBC, you need to convert to a compact representation, which requires first using the rectilinear NIC and then checking all neighbor cells for a closer image. Can anyone verify whether PLUMED is doing distances correctly for triclinic boxes?

Thanks!

John

[Giovanni Bussi]

Hi,

in principle there is the plumed2-git list. Anyway, that was mostly to distinguish plumed2 related questions historically, perhaps it is not that necessary now.

1. DHENERGY inherits from CoordinationBase and as such is already parallel: atoms are scattered on all processors (this is true for all atoms requested by PLUMED) and the sum is parallelized with both MPI and Openmp (if you set PLUMED_NUM_THREADS properly).

1A For your particular case, it might be worthwhile scattering only the 40 atoms and leave the 100k parallel. This should be supported someway in master branch, though probably not yet well documented. This require some more effort to dig in the code.

1B Try multiple time stepping: http://plumed.github.io/doc-master/user-doc/html/_m_t_s.html . If you have a functional form similar for DHENERGY you should be able to tolerate STRIDE=2 or STRIDE=3 (see http://pubs.acs.org/doi/abs/10.1021/ct5007086). This is straightforward and will give a boost x2 or x3. Remember to check effective energy conservation.

1C This is a pain as you guessed.

2 PLUMED should do minimal image correctly for skewed boxes. There is also a paranoia test for this (regtest/based/rt-make-1). You can also play with WRAPAROUND on water to see that images are done right. In case you find a discrepancy, I would first assume that gromacs dist() does it wrong :-). Actually, if you can find an example where PLUMED it is not computing minimal images correctly this should be considered as bug, please report it and we will be very happy to fix it.

Giovanni

 

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/ca88c5ce-1ab8-4862-ac5d-ff5f9912e142%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[John Haberstroh]

1A. I am looking into the source to see if I can scatter just the 40 atoms, thanks for the tip!

1B. Do you have any experience with MTS compared to Hybrid MC? Hybrid MC seems good to me because it is exact and tunable -- and the worst that can happen is that it slows down the biased sampling. MTS seems a little dangerous, though a stride of 2 doesn't feel too scary and seems like a nice little free boost :)

1C. Noted!

2. I am quite confident that gromacs dist() is implemented correctly for this situation -- I did an A/B test of my gromacs implementation against a separate python implementation, and I found that I had a very problematic conceptual misunderstanding of minimum image convention in skewed boxes, and when I fixed it, my results matched the gromacs results to a T. The PLUMED results look much more like my original python code. I'll see if I can come up with an exemplary test case that demonstrates the issue (if one does exist!).

--
You received this message because you are subscribed to a topic in the Google Groups "PLUMED users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/plumed-users/K2NlFi7Qwgc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to plumed-users...@googlegroups.com.

To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAPLm8ZJ9GWD9CW%3D0J-ydCzqF-zvV8kHZVi4QE5zQXTpTM04-2Q%40mail.gmail.com.

[Giovanni Bussi]

On Tue, Jun 28, 2016 at 1:19 AM, John Haberstroh <jhabe...@berkeley.edu> wrote:

1A. I am looking into the source to see if I can scatter just the 40 atoms, thanks for the tip!

1B. Do you have any experience with MTS compared to Hybrid MC? Hybrid MC seems good to me because it is exact and tunable -- and the worst that can happen is that it slows down the biased sampling. MTS seems a little dangerous, though a stride of 2 doesn't feel too scary and seems like a nice little free boost :)

HMC will not work for a large system. The acceptance will depends on fluctuations of the total energy and will rapidly go to zero.

If you look in the paper about MTS in plumed it is explained how to monitor the change in total energy arising by the MTS on the PLUMED bias alone. Anyway, also this number is likely too fluctuating to be usable in a HMC procedure (if CV involved all the solvent).

Finally, implementing a HMC is a bit complicated technically because of the way PLUMED is interacting with GROMACS, so even doing a quick test would require some non trivial coding.

 

1C. Noted!

2. I am quite confident that gromacs dist() is implemented correctly for this situation -- I did an A/B test of my gromacs implementation against a separate python implementation, and I found that I had a very problematic conceptual misunderstanding of minimum image convention in skewed boxes, and when I fixed it, my results matched the gromacs results to a T. The PLUMED results look much more like my original python code. I'll see if I can come up with an exemplary test case that demonstrates the issue (if one does exist!).

Try WRAPAROUND ATOMS=water-atoms GROPBY=3 AROUND=center-of-your-solute

You will see that distances are correctly computed in the rhombic dodecahedral cell

Giovanni

 

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAP3B6WVP7SqStVxPUqYhHPCnshPqZkn%3DWoeKMOH0844E0YYNNA%40mail.gmail.com.