Proposal for information flow for backmapping from a DPD simulation in HOOMD Blue

[mrshirts]

• Overall problem:

  • We start with a complete MuPT representation with both AA and CG (repeat unit) Primitives, which includes topological and hierarchical connections for these Primitives, initially without coordinates. We want to end up with CG and AA coordinates, both consistent with a dense melt.

• Detailed description of starting information:

  • We have complete representations with both CG and AA definitions, lacking only coordinates.
  • We have 1 or more AA sets of relative atomic coordinates, generated by small molecule conformer generators, for each AA, call these "repeat unit templates". These can be rotated and translated as rigid bodies as needed.
  • CG Primitives are optionally assigned an enveloping rigid body Shape, such as a Sphere, Ellipsoid, in addition to a specification of center of mass coordinates.

• Proposed procedure:

  • Determine three quantities from the AA templates:
    • D, The distance between COM of the CG descriptions
      • calculated by center of mass or distance between connectors for the AA units)
    • S, The size of the bounding sphere
      • calculated by the distance from the COM of the AA template to the furthest atom out.
    • B, Distance between connectors
      • Connectors are the midpoints on the bonds between two linked beads
      • Ddistance between the two anchor atoms in the AA template.
    • Write a GMSO CG XML file with the following parameters
      • DPD particle size parameter is S
      • DPD bond equilibrium length is D
      • DPD equilibrium angle theta for the 3 bead distances CG-CG-CG is 2 arcsin (B/D). This is because if the connectors are set at the midpoints of the bonds connecting each particle, then sin(theta/2) = (B/2)/(D/2). Force constant is somewhat, but not too weak, so equilibrium distance between connectors remains within 10% or so of B. Size S does not come into play here for this algorithm.
  • Generate a CG non-self avoiding angle-constrained random walk using D and theta. Size is not directly used here (S could be for a self-avoiding random walk, but we don't need to for the current DPD procedure).
  • This CG representation gets written to GMSO, and parameterized with the exported XML file.
    • A periodic box is defined that simply bounds all the particles - it will be rapidly shrunk in the DPD process, I think? Or do we just impose a box of the desired size, and not worry how the particle overlap? @StephMcCallum, can you answer this?
  • the DPD process runs in HOOMD, to obtain a dense melt in the desired periodic box. @StephMcCallum, is this an NVT or NPT process?
  • the CG positions (particle centers) are read/handed back into MuPT to inject into the CG representation.
    • Congratulations! There is now CG melt ready to write out to any supported engine (not just HooMD Blue).
  • the AA positions are are assigned by MuPT by rotating and translating the AA template(s) so that the AA connectors line up with the CG connectors. The DPD angle guarantees that the CG connectors remain at about B during the simulation. If the bonds are a little long, that is OK, the won't be TOO long.
  • the AA level of the representation is exported to GMSO or Interchange, and then to the desired simulation engine.
  • In the simulation engine, the AA configuration is parameterized and energy minimized.

• At this point, congratulations! There is now an AA melt ready to write out to any supported engine, which presumably will not crash.

[StephMcCallum]

This CG representation gets written to GMSO (but no coordinates yet - if we can define a GMSO w/o coordinates), and parameterized with the exported XML file.

Positions are required for GMSO. So we could do this step after assigning coordinates with a random walk and before running the DPD simulation. This seems like an extra step to me. Why would we want to export and parameterize in GMSO and then back to mupt?

[mrshirts]

Positions are required for GMSO

Well, these could be zero, but we can certainly do it inside MuPT (just move the code inside)- I was just trying to simplify things.

Why would we want to export and parameterize in GMSO and then back to mupt?

I don't think we are? The two GMSO steps are CG (which, if one wants to do CG, is the only output) and the output for the AA at the end. Let me know where it's confusing, and I will adjust.

[timbernat]

We want to end up with CG and AA coordinates, both consistent with a dense melt.

True, though more precisely we also want to store those positions in the representation tree self-consistently (i.e. atom and bead positions "agree")

• Detailed description of starting information:

I might fill in starting info in a little bit more detail, as:

• Detailed description of starting information:
  • We have MuPT representations storing topological and hierarchical connections for AA and CG (repeat unit) Primitives, initially without coordinates.
  • We can readily use small mol conformer generators to obtain correct relative atomic positions for each chemically-distinct CG Primitive, call these "repeat unit templates"
  • CG Primitives are optionally assigned an enveloping rigid body Shape, such as a Sphere, Ellipsoid, or just center of mass coord

then precede "Proposed Procedure" with:

• Loosely-packed initial CG Primitive positions are obtained by an angle-constrained random walk, which is cheap to evaluate (since no history needs to be tracked) and only requires an effective diameter for each unit
  • Primitives support rigid transformations, so the "reference" position for each CG Primitive is rigidly moved onto its corresponding random walk step, and the underlying Atom Primitives, with their correct relative positions, are moved in lockstep to their absolute positions

Determine three quantities from the AA templates

Strictly 4, since the center of mass "C" also needs to be (and can easily be) determined from the atomic positions

Connectors are the midpoints on the bonds between two linked beads

Only when the separation between adjacent CG beads is precisely 0. Connectors consist of an anchor point, which is coincident with the anchor atom's position, and a linker point, which should be coincident with the adjacent atom (and anchor) in the "next" bead.

In the most general case, one would have a trapezoid:

     A1 --- A2
    /        \
  C1 -------- C2

where C1, C2 are the centers of each bead and A1 is the anchor of the first connector AND the linker of the second, and A2 is defined similarly for the second bead, mutadis mutandis. However, the rigid transforms in the random walk (already in) roll each next bead to force these collinear, as:

C1 --- A1 ------- A2 ---- C2

Incidentally, getting the A1-A2 connectors to line up is also not given, and is what the bond alignment procedures (also already in) handle.

calculated by the distance from the COM of the AA template to the furthest atom out.

To align better with DPD, the diameter and center of the sphere should really be "B" and the midpoint of the two anchors on one bead, respectively.

DPD particle size parameter is S
DPD bond equilibrium length is D

These are determined by the parameters defined above, namely:

• S = B/2, and in general won't be the same for all beads. Since the LJ particles are points, this will get absorbed into the separation distance instead as:
• D = 1/2B_1 + b + 1/2B_2, where b = |A1 - A2|, the underlying atomic separation or equivalently the distance between the surfaces of the two beads along the center-to-center line.

The DPD angle guarantees that the CG connectors remain at about B during the simulation. If the bonds are a little long, that is OK

Exactly!

the AA representation is exported to GMSO or InterMol, and then to the desired simulation engine

I want to be precise that the AA parts of the representation are exported; we're not creating a new tree, but using only the coordinates and elements in the "bottommost" parts of the one we've just updated from the CG simulation. Also, are we using Intermol or Interchange? I though the former was deprecated

In the simulation engine, the AA configuration is energy minimized.

Do we then do an analogue of what happen above after the CG simulation, where we read back in and update the atomic positions in the repr?

[mrshirts]

True, though more precisely we also want to store those positions in the representation tree self-consistently (i.e. atom and bead positions "agree")

Yes, in the more general case! I'm just working on what we are currently doing, where we specifically want the CG to be a dense melt, and then generate an AA configuration that agrees. "agreement" is not entirely defined, as backmapping is not uniquely defined, and we can use other definitions of Coarse-graining that are not center of mass or center of geometry averaging. But yes, we want them to agree with whatever the the user or future developers implement.

[mrshirts]

Also, are we using Intermol or Interchange? I though the former was deprecated

Yes, I will edit!

[mrshirts]

I want to be precise that the AA parts of the representation are exported; we're not creating a new tree, but using only the coordinates and elements in the "bottommost" parts of the one we've just updated from the CG simulation.

Yes, I will clarify.

[mrshirts]

Primitives support rigid transformations, so the "reference" position for each CG Primitive is rigidly moved onto its corresponding random walk step, and the underlying Atom Primitives, with their correct relative positions, are moved in lockstep to their absolute positions

So, this is done at the end, right, after the simulation? When you have updated CG primitive coordinates, then the Atomic primitives can have their relative positions translated to be consistent with the CG Primitive, correct? But these relative positions are not carried around in the simulation.

I guess it could be done during the non-self-avoiding random walk, but not sure if it needs to be, since we are not going to do anything with the AA coordinates at this step (I don't think).

The relative orientations need to be set as well before you can actually do anything with the AA coordinates, as well, or there will be really long bonds that will cause problems.

[mrshirts]

S = B/2, and in general won't be the same for all beads. Since the LJ particles are points, this will get absorbed into the separation distance instead as:

S could be different from B/2. For example, the LJ sphere radius is frequently larger than B (bonded atoms are usually well-within the repulsive core for QM reason. For water, Sigma_O = 3.15/2 Angstroms, and B is 0.95 Angstroms. But in THIS algorithm I don't think we will need S.

[mrshirts]

Strictly 4, since the center of mass "C" also needs to be (and can easily be) determined from the atomic positions

that doesn't need to be defined FOR the XML, though? It's a 3d vector, not a parameter.

[mrshirts]

Connectors are the midpoints on the bonds between two linked beads

Only when the separation between adjacent CG beads is precisely 0.

No, I think this is always the case. S can be bigger, and I don't think is involved.

[mrshirts]

I think we are a little differing on the geometry. I see D as the distance between particle centers. B will be shorter than D, and will be equal to D only if the connectors are opposite each other. If the connectors are not opposite each other, then B will be less than D, and we impose the angle so that the bonds can be the right distance even if the centers of each particle are further apart. How does that picture work?

[timbernat]

When you have updated CG primitive coordinates, then the Atomic primitives can have their relative positions translated to be consistent with the CG Primitive, correct? But these relative positions are not carried around in the simulation.

Right, each bead gets a start-to-end of simulation transformation which also updates the underlying atom positions when passed back to the representation

So, this is done at the end, right, after the simulation?

The second time it happens, yes

I guess it could be done during the non-self-avoiding random walk, but not sure if it needs to be, since we are not going to do anything with the AA coordinates at this step (I don't think).

This is the first time where it happens, when beads are moved into place along the RW such that the anchors and centers are all in-line. You can't really avoid this step (only substitute it), since you need a starting configuration for the DPD simulation. The RW gives "reference" positions for each repeat unit an absolute position WRT neighbors.

The fact that an angle-constrained RW gets the relative orientations of neighbors basically right means DPD should converge more quickly, and really acts just to "densify" the loose random walk chains with one another.

that doesn't need to be defined FOR the XML, though? It's a 3d vector, not a parameter.

Good point, that's needed for the structure output to mbuild, but not directly in the XML

I think we are a little differing on the geometry ... How does that picture work?

I've included a napkin sketch with some pictures. There are at least two different ways to pick a bounding sphere; with diameter from COM to anchors (1) or with diameter between anchors directly (2); I'm describing the latter, and I think you're describing the former. There are at least a couple reasons I prefer (2), namely:

• The P3HT example we looked at, where long sidechains skew the COM and force repeat units with the correct atomic inter-monomer bond distances to overlap as CG spheres
• As depicted in (4), the COM will in general NOT be equidistant to both anchor sites; you'd need to project it onto their bisector plane to get both anchors to lie on the surface of the sphere, or else always have one anchor floating inside or outside of the sphere

As I look at it though, I realize the drawback of (2) is that the center and anchors are ALWAYS collinear, whereas we'd want each 2-functional bead to have a characteristic angle. That would let you get the "trapezoid" alignment I described, drawn in (3).

Come to think of it, this ambiguity about which distance to use might be why the bonds I got from the draft DPD placer were so long. Good to make this more precise!