Go-like Models in HOOMD

[Joshua]

I'm about to embark on doing a bunch of simulations on a Go-like
polymer model and am trying to pick a simulation package and was
considering using HOOMD. Go-like models are essentially a bead-spring
polymer model except that non-bonded terms between two beads are
either (1) LJ-like if the two beads make contact in a reference
conformation, or (2) purely repulsive (e.g. eps*(sig/r)^12 - I think
you can get this from the current potentials by shifting the standard
LJ potential and applying the appropriate cutoff) otherwise. I think
this can be pulled off if every bead has a unique identifier/type.
These types of models are generally used to look at protein folding to
some native state (which, by definition in this model has the lowest
energy).

We are trying to gather extensive sampling/statistics for a folding
event in a polymer that is on the order of a couple of hundred beads.
My understanding from looking at the mailing list and looking at the
benchmarks, is that the GPU performance can be negatively impacted by
using a small number of particles (correct me if I'm mistaken about
this). I was wondering if it was possible to play a little trick to
effectively increase the system size in order to boost performance.
Since there are no long range interactions in this system, I was
thinking you could build a lattice of the polymers, harmonically
restraining a single atom in each replica to a point in space such
that each realization of the polymer was separated from the others so
that they don't interact. You could then scale the number of particles
in the system up to the sweet spot for the GPU.

Any suggestions, feedback, advice pertaining to simulating this sort
of system with HOOMD would be appreciated. It looks like a great
package and it would be great to be able to use it for our work.

[Axel]

hi joshua,

On Nov 25, 4:02 pm, Joshua <joshua.adel...@gmail.com> wrote:
> I'm about to embark on doing a bunch of simulations on a Go-like
> polymer model and am trying to pick a simulation package and was
> considering using HOOMD. Go-like models are essentially a bead-spring
> polymer model except that non-bonded terms between two beads are
> either (1) LJ-like if the two beads make contact in a reference
> conformation, or (2) purely repulsive (e.g. eps*(sig/r)^12 - I think
> you can get this from the current potentials by shifting the standard
> LJ potential and applying the appropriate cutoff) otherwise. I think

in HOOMD the latter is even easier. the 12-6 LJ potential has
an "alpha" parameter, setting it to 0.0 will retain only the repulsive
part. HOOMD in its current form only supports one global cutoff.

> this can be pulled off if every bead has a unique identifier/type.
> These types of models are generally used to look at protein folding to
> some native state (which, by definition in this model has the lowest
> energy).
>
> We are trying to gather extensive sampling/statistics for a folding
> event in a polymer that is on the order of a couple of hundred beads.
> My understanding from looking at the mailing list and looking at the
> benchmarks, is that the GPU performance can be negatively impacted by
> using a small number of particles (correct me if I'm mistaken about

yes. if i recall correctly, a high-end nvidia GPU has 30 multi-
processors,
each with 8 cores and each of them will be scheduled to execute
4 threads simultaneously (similar to hyperthreading) which means that
about a thousand threads need to be executed at the same time to
keep it busy. if you also consider the overhead of launching a kernel,
that you need about 10x as many threads to make that overhead
become less important. thus a decent performance of HOOMD requires
of the order of 10,000 particles to reach its full speed.

> this). I was wondering if it was possible to play a little trick to
> effectively increase the system size in order to boost performance.
> Since there are no long range interactions in this system, I was
> thinking you could build a lattice of the polymers, harmonically
> restraining a single atom in each replica to a point in space such
> that each realization of the polymer was separated from the others so
> that they don't interact. You could then scale the number of particles
> in the system up to the sweet spot for the GPU.

in principle it should be possible to do something like that, but
there
are currently no position restraints. it would be better to have one
of
those on the center of mass of a group of particles to avoid
artifacts.

HTH,
axel.

[Joshua Anderson]

On Nov 25, 2009, at 1:01 PM, Axel wrote:

> hi joshua,
>
> On Nov 25, 4:02 pm, Joshua <joshua.adel...@gmail.com> wrote:
>> I'm about to embark on doing a bunch of simulations on a Go-like
>> polymer model and am trying to pick a simulation package and was
>> considering using HOOMD. Go-like models are essentially a bead-spring
>> polymer model except that non-bonded terms between two beads are
>> either (1) LJ-like if the two beads make contact in a reference
>> conformation, or (2) purely repulsive (e.g. eps*(sig/r)^12 - I think
>> you can get this from the current potentials by shifting the standard
>> LJ potential and applying the appropriate cutoff) otherwise. I think
>
> in HOOMD the latter is even easier. the 12-6 LJ potential has
> an "alpha" parameter, setting it to 0.0 will retain only the repulsive
> part. HOOMD in its current form only supports one global cutoff.

If a true WCA potential is required, you will not need to wait long. I
have that functionality working in the template-pair-potentials branch
now.

>> this can be pulled off if every bead has a unique identifier/type.
>> These types of models are generally used to look at protein folding
>> to
>> some native state (which, by definition in this model has the lowest
>> energy).

There is a limitation of around 40 particle types in hoomd. It stores
all N^2 different type parameters and the limit is imposed by the
amount of fast cache memory. Surely you should be able to implement
the LJ/WCA combination with only 2 particle types?

>>
>> We are trying to gather extensive sampling/statistics for a folding
>> event in a polymer that is on the order of a couple of hundred beads.
>> My understanding from looking at the mailing list and looking at the
>> benchmarks, is that the GPU performance can be negatively impacted by
>> using a small number of particles (correct me if I'm mistaken about
>
> yes. if i recall correctly, a high-end nvidia GPU has 30 multi-
> processors,
> each with 8 cores and each of them will be scheduled to execute
> 4 threads simultaneously (similar to hyperthreading) which means that
> about a thousand threads need to be executed at the same time to
> keep it busy. if you also consider the overhead of launching a kernel,
> that you need about 10x as many threads to make that overhead
> become less important. thus a decent performance of HOOMD requires
> of the order of 10,000 particles to reach its full speed.

Axel's quick analysis is indeed correct. I've done some speedup
benchmarks vs. system size: look at slide 24 here http://gladiator.ncsa.uiuc.edu/PDFs/accelerators/day3/breakouts/chem/anderson.pdf
The speedup is about 20x at 5k particles, 40x at 10k and 50x at 20k.
So while hoomd doesn't run at optimal efficiency with small numbers of
particles, 5k is still large enough to get decent speedups.

>> this). I was wondering if it was possible to play a little trick to
>> effectively increase the system size in order to boost performance.
>> Since there are no long range interactions in this system, I was
>> thinking you could build a lattice of the polymers, harmonically
>> restraining a single atom in each replica to a point in space such
>> that each realization of the polymer was separated from the others so
>> that they don't interact. You could then scale the number of
>> particles
>> in the system up to the sweet spot for the GPU.
>
> in principle it should be possible to do something like that, but
> there
> are currently no position restraints. it would be better to have one
> of
> those on the center of mass of a group of particles to avoid
> artifacts.
>
> HTH,
> axel.
>
>> Any suggestions, feedback, advice pertaining to simulating this sort
>> of system with HOOMD would be appreciated. It looks like a great
>> package and it would be great to be able to use it for our work.
>

> --
>
> 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
> .
>
>
>
>

[Joshua Adelman]

Hi Josh,

Thanks for your response. I think the limitation on the number of particle types may pose a problem for me. In a Go-type model two beads only interact through a LJ potential if they are within some cut-off distance (7-8 angstroms) in the reference structure, which is usually the folded structure determined via x-ray crystallography. The parameters are chosen in such a way to match the distance between the pair of beads to the distance in the reference structure. All other interactions are purely repulsive. I believe this means that you would need to define a unique interaction between each pair of atoms in the system. In my case that would translate to several hundred particle types.

I'll try to think if there is a way to get around the particle type limitation on the model side of things. Hopefully there is, because I was quite excited about HOOMD for this project.

Also, do you have any intention of adding positional restraints to HOOMD's functionality soon? I'd be interested in something similar to the following functionality found in LAMMPS:
http://lammps.sandia.gov/doc/fix_spring.html

Best wishes,
Josh

[Joshua Anderson]

I understand the issue with the particle types now. With just a little
bit of c++ hacking, an easier method should be possible. Simply apply
the repulsion term to all particle pairs via the standard
neighborlist / pair force methods in hoomd. Then construct a fixed
"neighbor list" embedding the connectivity information from the
reference structure and apply the attractive part of the potential
with that list.

If your connectivity has less than 12 neighbors per particle, it can
be done a little simpler. Apply the WCA potential between all
particles (single type again), but add the connected particles to the
neighbor list exclusions. Then apply the standard LJ potential with
the fixed "neighbor list" as above. This is less work because you
don't have to create a new pair potential that is just the attractive
part, and you can use PotentialPairGPU<PairEvaluatorLJ> for both. The
only additional c++ code required is to derive from NeighborList a
simple class that reads in your connectivity information from a file
(for example).

The functionality of fix/spring (and more) is planned. There is no
current ETA
https://codeblue.umich.edu/hoomd-blue/trac/ticket/286

--------
Joshua A. Anderson, Ph.D.
Chemical Engineering Department, University of Michigan

[Joshua Adelman]

I'm not sure that the hack solution will work in this case. Since the
force-field we are trying to construct is meant to create a global
energy minimum corresponding to the reference conformation of the
polymer, the LJ term between a particular pair of beads needs to have
a specific \sigma_{i,j}, rather than a generic \sigma for all "native"
contacts (e.g. two different pairs of atoms have a desired minimum at
two distinct distances). This also impacts the parameterization of the
WCA potential in the hack (in the implementation the purely repulsive
does not need to be WCA, though and could just be the 12 term in the
standard LJ potential with a single \sigma for all "non-native" pairs).

[Joshua Anderson]

There is a way to remove the limit on the number of particle types. It
would be relatively easy to modify the lj potential to read the pair
coefficients from the texture cache. You would still have O(Ntypes^2)
storage requirements, but that would only limit you to potentially
thousands of types. Performance would suffer: probably 25-30% overall.
Maybe Fermi's larger caches would help mitigate the damage, maybe not.

--------
Joshua A. Anderson, Ph.D.
Chemical Engineering Department, University of Michigan

[Joshua Adelman]

We might be able to live with that sort of performance hit since the
algorithm that we are using is trivially parallel, similar to the
Replica Exchange Methods that Axel coded up.

Once the particle type limit has been removed, would there be a
secondary problem if we want to define many unique sets of bond terms,
angle terms, etc for pairs and groups of those particle types? For
example we have a linear polymer with beads sequential typed as T1,
T2, T3 ... TN. We want to define a unique (i.e unique parameters for
that specific interaction) bond between each adjacent pair of atoms in
the chain, a unique angle spring between every 3 adjacent beads, and a
unique dihedral between every 4 adjacent beads.

Thanks again for all of your input on this.

>>>>>> To post to this group, send email to hoomd-
>>>>>> us...@googlegroups.com.

[Joshua Anderson]

On Nov 30, 2009, at 11:52 AM, Joshua Adelman wrote:

> We might be able to live with that sort of performance hit since the
> algorithm that we are using is trivially parallel, similar to the
> Replica Exchange Methods that Axel coded up.
>
> Once the particle type limit has been removed, would there be a
> secondary problem if we want to define many unique sets of bond terms,
> angle terms, etc for pairs and groups of those particle types? For
> example we have a linear polymer with beads sequential typed as T1,
> T2, T3 ... TN. We want to define a unique (i.e unique parameters for
> that specific interaction) bond between each adjacent pair of atoms in
> the chain, a unique angle spring between every 3 adjacent beads, and a
> unique dihedral between every 4 adjacent beads.

Bond/angle/dihedral parameters are already read from textures, so
there should be no limitation.

> Thanks again for all of your input on this.

No problem.