Implementing a Harmonic Length Restraint for Linear Molecules

LAMMPS LAMMPS Development
1

Hello everyone,

I have a linear molecule, where I have defined, in the molecule file, the atom types of the linear molecule to be in consecutive order, and I have a NVT simulation containing and only containing multiple copies of this linear molecule, where the endpoints of each copy are held stationary.

I want to implement a harmonic force to restrain the total length of the molecule to a certain length, l. Let 1 be the atom type of one endpoint of the molecule, and n be the atom type of the other endpoint of the molecule.

For atom types 1 through n - 1, this restraining force will have a positive contribution of k times (total length of the molecule the given atom belongs to / l - 1) times (position vector of given atom minus position vector of successor atom / distance between the given atom and its successor), where k is a scaling constant.

For atom types 2 through n, this restraining force will have a negative contribution of k times (total length of the molecule the given atom belongs to / l - 1) times (position vector of given atom minus position vector of predecessor atom / distance between the given atom and its predecessor). k here is the same as in the previous paragraph.

Several problems with implementing this immediately arise:

  • Because I am trying to restrain the total length of a molecule, instead of constraining the length of specific bonds, this means that fix shake goes out the window.

  • For fix addforce, implementing this force is simple if only one copy of the molecule exists in the system: I can use compute bond/local to calculate the magnitude and components of the bond lengths of said molecule, sum up the magnitude of the bond lengths to get the total length of the molecule, and then feed this information into an atom-style vector for fix addforce.

  • However, if multiple copies of the molecule exist in the system, because the restraining force on any atom cannot depend on the total lengths of any molecule the atom does not belong to, this would require me to explicitly define a group for every copy of the molecule that exists in the system - which obviously would not scale well.

  • I’ve taken a look at compute chunk/spread/atom to see if I could use it to distribute the total lengths of each molecule to their respective atoms, and it would not actually work because it would require a reduce/chunk of the bond distances, but reduce/chunk requires a atom style variable as an input and compute bond/local, which can give me the bond distances of every bond in the system, does not return an atom style input. This isn’t the only obstacle I’ve encountered so far with trying to implement this using only LAMMPS script commands.

  • The option that my PI wants me to pursue is to write a custom fix and/or bond source file for LAMMPS. However, there are several facts that I do not know how to deal with: the fact that different atoms of the same molecule can be locally owned by different processes, and reassigning forces back to their respective atoms across the different processes (the fact that this is my first time working with MPICH doesn’t help my situation here)

This leads me to my question: Is it possible to implement the harmonic length restraint, as I have described above, entirely in the LAMMPS script (without sacrificing too much parallel performance)? If not, how should I be structuring my source file for this harmonic length restraint?

2

Some general comments:

What you subsequently describe is not a “harmonic force to restrain the total length of the molecule to a certain length”. That harmonic potential is very simply expressed as

U(\mathbf{r}_1, \mathbf{r}_n) = \frac{1}{2} k (|\mathbf{r}_1- \mathbf{r}_n| - r_0)^2

and is easily inserted into your simulation by adding a harmonic bond between the first and last beads of each molecule with the desired harmonic strength. Note how, since the length of the molecule does not (directly) depend on the locations of intermediate beads, a harmonic restraint on the molecule length can’t impose a direct force on those beads either for Hamiltonian consistency.

For more help, you should give the description of the forces in mathematical notation. When you say “positive contribution” and “negative contribution”, for example, that’s mathematically ambiguous (to me at least).

To “spread” the molecular length to an entire molecule, you can try something like this:

# group "ends" contains only terminal particles of molecules

compute mchunk all chunk/atom molecule
compute r2mols ends gyration/chunk mchunk
compute r2mol_peratom all chunk/spread/atom mchunk c_r2mols[4]

# r2mol_peratom now contains (half) the molecular
# end-to-end distance squared as a peratom value
# (possibly additionally scaled depending on particle masses)

Read the documentation carefully for compute gyration/chunk to see why this works (or should work – test before using!).

See this post: Propel polymer bead along bond vector - #4 by srtee for more very hacky advanced ways of using chunks.

This phrase is historically correlated with significant headaches, associated with a PI without sufficient knowledge of what’s “easy” and “hard” in MD, and thus suggesting things without a clear budget for either the programming cost or the physical-significance payout.

Your best bet is to find a local mentor or collaborator to keep you on the right path. A straightforward way to identify such a person is to find the administrator of the cluster you’re intending to run these simulations on, and ask them who to talk to.

3

Positive contribution to what?

Negative contribution to what?

Two points on that: a) there is a difference between a restraint and a constraint. A restraint is an added force that depends on how far the system is from a target value. A constraint is a manipulation that tries to maintain a specific geometry, i.e. a bond constraint manipulates the forces and velocities so that two atoms remain at a given distance with a very small margin. b) I don’t see what fix shake would not be applicable in your case.

Fix addforce is not a good choice here.

That would not be necessary in my opinion, since there are multiple ways to apply a constraint/restraint.

That is only a problem, if the atoms are further apart than the communication cutoff.

This is usually handled outside of the code you need to contribute, so you won’t have to programm that. and the only parts that you would need MPI for would be the restart file handling, and that can be omitted in this case.

That said, I see multiple existing options to address this issue:

  • you use the fix restrain command — LAMMPS documentation. The command line probably would have to be generated based on some script parsing your molecular topology
  • you define an additional “long” bond between the first and the last atom of your molecule at the desired length. If that bond is of “harmonic” style you implement a restraint, but you can also apply fix shake to this one single XXL bond to constrain its length. This requires that the communication cutoff is long enough.
  • you can use pair style list where you add this to your regular pair style using pair style overlay and also add bonds with it to the first/last atoms. This does not require changes to the communication cutoff but makes LAMMPS be subject to minimum image conventions.
4

Hello srtee,

I think I wasn’t being clear enough about what I meant by “length of the molecule”. It seems that by your equation U(r_1,r_n) = \frac{1}{2}k(|r_1-r_n|-r_0)^2, you interpret “length of the molecule” to be the displacement between the endpoints of the molecule. However, that is not what I am trying to restrain - what I am trying to restrain is the arc length of the molecule between the endpoints.

Now that I finally know that LaTeX works here, I can denote what I mean by “positive contributions” and “negative contributions” to the restraining force:

For the sake of simplicitly, all atom indices in the following notation refer to atoms belonging to the same molecule, and the atoms in this molecule are indexed sequentially, where one of the endpoints of this molecule will be indexed 1, and the other will be indexed n.

I will denote d(i,i+1) to be the distance between atom with index i and i+1. This means that the arc length of the molecule will be equal to:

\left( \sum_{i=1}^{n-1}d(i,i+1) \right)

I will denote l_{tot} to be the value that I want to restrain the end-to-end arc length of the molecule to.

I want to implement a quadratic penalty for when the arc length of the molecule is greater or lower than l_{tot}, where the energy associated with this quadratic penalty is:

U = \frac{k}{2} \frac{\left( (\sum_{i=1}^{n-1}d(i,i+1)) - l_{tot} \right)^2}{l_{tot}}

The force, on each particle, that arises from this quadratic penalty is going to be equal to the derivative of the restraint energy with respect to the position of said particle.

For atom 1, I’ve calculated this force to be:

k \left(\frac{\sum_{i=1}^{n-1}d(i,i+1)}{l_{tot}} - 1 \right) \frac{\textbf{x}_1 - \textbf{x}_2}{d(1,2)}

For atom n, I’ve calculated this force to be:

-k \left(\frac{\sum_{i=1}^{N-1}d(i,i+1)}{l_{tot}} - 1 \right) \frac{\textbf{x}_n - \textbf{x}_{n-1}}{d(n,n-1)}

For atoms 2 through n-1, I’ve calculated the force for each of those atoms to be:

k \left(\frac{\sum_{i=1}^{N-1}d(i,i+1)}{l_{tot}} - 1 \right) \left(\frac{\textbf{x}_i - \textbf{x}_{i+1}}{d(i,i+1)} - \frac{\textbf{x}_i - \textbf{x}_{i-1}}{d(i,i-1)} \right)

From these expressions, what I meant when I said

For atom types 1 through n - 1, this restraining force will have a positive contribution of […]

is that for atom types 1 through n - 1, the force from the arc length restraint acting on each of these atoms will have a \frac{\textbf{x}_i - \textbf{x}_{i+1}}{d(i,i+1)} term.

and what I meant when I said

For atom types 2 through n , this restraining force will have a negative contribution of […]

is that for atom types 2 through n , the force from the arc length restraint acting on each of these atoms will have a - \frac{\textbf{x}_i - \textbf{x}_{i-1}}{d(i,i-1)} term.

I hope this clears things up as to what I want to accomplish.

5

Hello Axel,

As I have mentioned in my reply (which should be above this post) to srtee, I am not trying to restrain the displacement between the endpoints of the molecule - what I am trying to restrain is the arc length between the endpoints of the molecule, and that is the only thing about the molecule that I want to restrain. This should explain why fix restrain and fix shake would not accomplish what I want to do - I am not trying to restrain specific bonds/angles/dihedrals within the molecule, but rather the arc length across the entire molecule.

6

Then you should say so. Your original post did the opposite, specifically the following is extremely misleading and in conflict with your second post.

From your revised description, you do want to restrain the length of individual bond (i.e. add a force that pushes those bonds to be longer or shorter), just you want that to be computed dynamically.

Also, you are arguing with atom types, but bonds generally pay no attention at all to atom types.

So the major challenge here would be to obtain the total length as the sum of the individual bond lengths. This is best done with a custom bond style where you would need to add the following kind of code:

  • each molecule should have a unique molecule ID
  • to compute your kind of molecule length do the following
    • determine the largest molecule ID first locally and then via MPI_Allreduce() with MPI_MAX
    • use that value to allocate an array of double and initialize to zero to store a single double value for each molecule ID
    • loop over the list of (local) bonds for the first time to determine the length of each bond and add it to the entry of the molecule ID that bond “belongs” to (i.e. pick the “molecule” property of one of the atoms). This needs special care when newton is set to off for bonds, because then bonds straddling subdomain boundaries are listed twice. So in that case you only add the bond length when the first bond atom is local and the second is a ghost, otherwise add all of the contributions.
    • do an MPI_Allreduce() with MPI_SUM for for your custom array and you have the “length” of each molecule on all MPI processes
  • With this information you should be able to do a second loop over the bonds where you first compute the bond forces as you currently do, i.e. copy from the bond style you are already using. And then compute your additional term that is dependent on the length per molecule.
  • It would be beneficial to create this bond style as a derived class, since then you only need to reimplement a few functions like the compute() function.
7

As a non-LAMMPS comment, this model doesn’t make sense to me, and it would require pretty vigorous defending in review.

After all:

  1. the bond lengths of individual harmonic bonds are Gaussian-distributed;
  2. the arc-length of a simple “harmonic polymer” is just the sum of the individual bond lengths;
  3. the sum of Gaussian-distributed variables is itself approximately Gaussian-distributed
    • (becoming more and more so as the polymer lengthens – the Central Limit Theorem!);
  4. and therefore the polymer arc-length already behaves effectively harmonically.

The model you propose adds some kind of global effect where, as you can tell by your own derivation, each individual bond somehow magically knows what the average stretch of all other bonds is and becomes looser or tighter to suit. Does this change for polymers of different lengths? What if a very long polymer has some confined regions and some free regions? Who knows.

Having said that, notice how your proposed interaction simply takes all existing bond forces and scales them by some factor that is uniform for each molecule. So, to simplify the development process as much as possible, here’s what I’d suggest:

  1. Code up some kind of compute bondlength/atom that, for each atom, runs over all bonds in which the atom participates, adds up the length of each bond, and outputs the summed length as a per-atom variable.
    • You need to be careful about bonds across domains (but you’d have to be careful anyway): if newton_bond is on, some processors need to send their “ghost atom” data back to the owner’s processor. But that’s all the MPI work you have to do, and it’s all handled through LAMMPS’s compute wrapper functions.
  2. Use compute reduce/chunk to sum the total bond lengths across each molecule and create a per-atom variable dividing that number in two. (This saves you the hassle of writing all that MPI code yourself.)
  3. Write a fix bond/scale fix that takes a per-atom variable and, during post_force(), calculates what the bond forces were and subtracts some proportion of it from the x, y and z forces, as supplied by the per-atom variable.

Using this approach you have three separate small code fragments that you can individually test to minimize the chance of bringing serious bugs into your production run.

1 Like