Harmonic and Table Dihedral functions

[Nathan Horst]

Hello HOOMD-users,

A colleague and I have been working on implementing dihedrals into our system, but have come across a few problems when verifying that they are working properly.

Our dihedral potential is implemented as is shown in the attached file "potential.png" for a system that consists of a chain of particles 32 beads long. There is an energy well centered about dihedral angles of 0 corresponding to the trans configuration of the polymer chain, as well as wells corresponding to the gauche+ and gauche- configurations at angles of +/- pi/3, respectively. 

With this information in mind, we would expect that we would see a transition from an all-trans configuration at low temperatures to a configuration with equal probability for trans and gauche at high temperatures. We have tried implementing this potential using both dihedral.harmonic() and dihedral.table() by running simulations interpolating from 20K to 500K and we do not see what we expect at all. I have attached a figure entitled "dihedral_pot_comparison.png" in which you can see that the percentage of dihedrals in the trans configuration drops immediately even at low temperatures for both functions considered.

We ran further tests, where specifically we looked at cases where we ran simulations with dihedral.table() where we used our own potential, as well as a simple harmonic potential as given in the HOOMD example scripts. Running these systems at a constant temperature of 10K once again gave us results contrary to what we expected. I have attached a file "dihedral_table_10Kconstant.png" in which we show both versions of the dihedral table potential failing where a simulation without any dihedral angles stays in the trans configuration. This seems to suggest that no matter how we implement dihedral angles into the system, the system will be in error.

We looked into further debugging this problem, by running the HOOMD unit tests for the dihedral potentials. However, upon further inspection, these unit tests simply seem to check whether or not the potential is passed through the simulation without failing and do not incorporate any testing of how accurately the potentials are incorporated.

Does anyone have any example they could share with us where the dihedrals have been working as expected? Alternatively, does anyone know where we might be in error? I have attached a runfile (simtestchains.py) and xml file (phi1.xml) so anyone can try this same system out on their own machines.

Thanks,

Nathan Horst

[Joshua Anderson]

Thanks for the detailed report. I don't have time to look into this right now, perhaps someone else on the mailing list can help.

I just wanted to say that the dihedrals are unit tested with C++ tests:
<https://bitbucket.org/glotzer/hoomd-blue/src/86926b80e67991b7fe86d40b4e26c1770edbe762/hoomd/md/test/test_harmonic_dihedral_force.cc?at=master&fileviewer=file-view-default>
<https://bitbucket.org/glotzer/hoomd-blue/src/86926b80e67991b7fe86d40b4e26c1770edbe762/hoomd/md/test/test_table_dihedral_force.cc?at=master&fileviewer=file-view-default>

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
Phone: 734-647-8244
http://www-personal.umich.edu/~joaander/

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
To post to this group, send email to hoomd...@googlegroups.com.
Visit this group at https://groups.google.com/group/hoomd-users.
For more options, visit https://groups.google.com/d/optout.

[Michael Howard]

Hi Nathan -- I don't see anything obviously wrong in your script, but it's a little tricky to debug the full system that you have because there is a lot of physics going on. I might suggest you go to a simpler test case -- a 4 bead chain with just bonds, angles, and dihedrals where all beads are the same type -- and then sample the distribution of the dihedral angle. For this test, it would be good to also use units so that temperature, energies, etc. are of order 1 (yours are very small -- your max temperature is order 1e-2, for example) just to make sure there are no precision errors at fault (are you running single or double precision?). Can you try doing such a test and let us know what happens?

Regards... Mike

[Nathan Horst]

Hi Mike,

I have attached the script we used for a chain of 4 particles(run_four_chain.py), as well as the initial xml file (scaled_four_chain.xml). We left multiple names for the particles for the purpose of our analysis scripts, but you will notice that the interactions are all set to be exactly equal. We noticed some interesting things when looking at the angle distribution from the simulation results, namely that it seems the dihedral potential wants to be flipped opposite to what we expect. I have included two figures for two opposite potentials which show the angle distribution throughout the simulation, as well as the regions that are expected to be the lowest energy regions according to the potential we defined(Potential_AngleDistribution.png and OppositePotential_AngleDistribution.png). Is this something as simple as a sign error?

Best,

Nathan Horst

[Curt Waltmann]

Hi MIke,

I'm Nathan's colleague who he referred to in the original post. In addition to testing the four-bead system we have now tested our original 32 bead system with bonds, angles, and dihedral angles using the newly released HOOMD 2.0.  Again we had no issues until we added the dihedral angles. Running HCexample_script.py on our transchain.xml using the same potential (potential.png from the original post) we saw a much lower percentage of dihedral angles in the trans conformation than we expected just like we did when using earlier versions of HOOMD. We graphed this in normal_dihedral_potential_20-500k.png. We also checked to see if we would get the same behavior as in our second post when we switched the sign of our dihedral potential. In reverse_dihedral_potential_20-500k.png, the dihedral angles behave like they did in our second post by going completely into the trans configuration despite the fact that our potential should now be dictating that they don't enter the trans configuration. It is also worth noting that all of our energy terms have been scaled up 30 times since the original post. Again we are not sure if the problem is as simple as a sign error that got transferred over to HOOMD 2.0 or something else in our files, but we would appreciate any ideas you or anyone else could give us.

Best,

Curt Waltmann and Nathan Horst

[Michael Howard]

Hi Curt, Nathan -- Thanks for running the additional tests. I am out of office on vacation until later this week, but I will take a look as soon as I can.

--Mike

[Michael Howard]

Hi Curt, Nathan -- I created a simple test case for the harmonic dihedral (which probably just replicates one of the unit tests... attached anyway), where I fix the position of three particles, and then log the energy of a harmonic dihedral as I rotate the fourth particle around a unit circle. When I do this, I get perfect agreement between the expected potential and the measured one. Of course, one could also do a proper run to also test to check if the forces on the particles are correct, but I imagine they would be (and that the unit tests are validating this). I know that we checked that for the OPLS dihedral at any rate.

This might be a stupid question, but are you sure that whatever code you're using to analyze the dihedrals isn't shifted by 180*?

Regards... Mike

[Nathan Horst]

Thanks for the response!

I think that our analysis code works fine with the test we ran, though we will certainly check it again to make sure. 

After looking at how you did this test, I ran it as you did and got the same results, that the energy seems to match the potential given. 

I then went ahead and made a few modifications to your scripts, and once again did not see what I would have expected, as I will try to explain below.

The potential that we want to replicate in our work is of the Ryckaert-Bellamans form where it is listed as the sum from i=0 to i=5 of a(i)*(cos(x)^i) where a(i) are constants.

We noted that there were not enough terms in the HOOMD OPLS potential to include this potential, without considerable error, so instead we modified the potential to be 5 separate potentials of the form specified in dihedral.harmonic(), where the only difference is a constant term for a shift in energy. We plotted the addition of these five potentials for use in HOOMD against the original form, and as expected saw only a constant shift in the potential (dihedral_potentials.png). Our approach, then, became defining five different dihedral types for each dihedral in our system assuming that their incorporation would be additive, but perhaps HOOMD is not recognizing this?

When I modified your scripts, I added five dihedral types, and set five corresponding coefficients. (I am not too familiar with snapshots, but I think that I have implemented this correctly.) Then I plotted the results of the HOOMD run, as you had done, and there seems to be a discrepancy between the energy of the dihedral calculated by HOOMD and the one we assumed (perhaps wrongly) would be implemented. The script and results are attached as well. I am not sure if my modified test is in error, or if we are unable to call the potential this way, but would appreciate your insight. In the meantime, I will be thoroughly vetting our analysis codes.

Thanks,

Nathan

[Eric Jankowski]

There are certainly inconsistent conventions between how people write down potentials and implement them. This discussion sounds like it might resolve with "Oh, we call 180 cis but they call it trans" like it did for some of our simulations a few months ago.

Eric

[Michael Howard]

Hi Nathan -- Thanks for running this additional test. I am confirming that this script also fails on my cpu system. But, I think I may have spotted the problem (I am not sure if this is by design or a bug, another developer will have to comment on how the dihedrals are implemented). Properly, you really only want to define one dihedral between those four atoms with multiple potentials acting on it, so I switched my original script so that the single harmonic dihedral is now replaced by multiple potentials acting on the dihedral of type "D":

d1 = dihedral.harmonic()

d1.set_coeff('D', k=2.363, d=-1, n=1)

d2 = dihedral.harmonic()

d2.set_coeff('D', k=1.578, d=1, n=2)

d3 = dihedral.harmonic()

d3.set_coeff('D', k=2.55, d=-1, n=3)

d4 = dihedral.harmonic()

d4.set_coeff('D', k=0.789, d=1, n=4)

d5 = dihedral.harmonic()

d5.set_coeff('D', k=0.4735, d=-1, n=5)

The resulting plot of the potential looks like what you were going for to me.

I don't know why the two aren't equivalent -- really they should be -- but maybe someone else can debug / answer that.

Regards... Mike

[Nathan Horst]

Thanks for the help! 

We will take this new syntax to our simulations, and see if that changes the results as we would expect, but at first glance it seems likely that it might.

Best,

Nathan

[Nathan Horst]

UPDATE: We ran this same test, where we set 5 separate harmonic dihedrals, and saw that indeed, as Mike showed, the system took the dihedral potential that we implemented. We then took this syntax to our simulations, in both HOOMD 1.3.3 and 2.0.0 and unfortunately, saw the same things that we did before--that there seems to be some problem with implementing dihedrals when the system is integrated.

We then ran some additional tests, and found some unexpected results, which I will try to clearly convey:

First we took a four bead system with a single dihedral and added bonds and angles to it, such that the beads started in the lowest energy state of all three potentials.

         At this point, if we run Mike's test, we see the same energy as he, since the bond and angle potentials are at their 0 point......no issues yet.

Next, we ran the nvt integrator at dt=0.0001 and kT=0.05 (what we understand to be ~20K in our units), dumped an xml file for the atom positions (which had not moved due to their existence in the minimum of the potential and a low kT value) and ran the same cycle through the dihedral angles to check the system energy. Strangely, despite the particles not moving, the energy began to change drastically with increasing simulation length. I have included the hoomd script (simlengthtest.py), final xml file (atoms10000000.xml), and plot of energies for the final state in simulations of increasing length (energy_by_simlength.png).

Any ideas as to what might be causing this? It certainly seems that HOOMD does not like the combination of dihedral potentials with integrators. Indeed, if we run this same system sans dihedrals, the system energy remains zero throughout when we cycle through the angles, as the particles are in the expected minima of the bond and angle potentials.

[Eric Irrgang]

I expect the energies are all wonky because you have a non-zero timestep while you are looping over angles but you are not restoring the positions of the first three particles after the run(1). You are flinging the molecule around by moving the fourth particle in a prescribed circle.

As to why the simulations of different lengths behave differently, I assume that has to do with the internal state of the thermostat that has been trying to add velocity for the whole simulation but can't because there is no velocity to scale.

[Michael Howard]

Eric is right, the file I supplied was to test the form of the potential, not to run a simulation. If you move the fourth particle but take a nonzero timestep, the other three particles have to move because the force of the dihedral is distributed among them (think of it like whipping a jump rope from one end). There is no reason to think this will give you back the potential. Things are even worse with the Nose-Hoover (NVT) thermostat, because it is trying to correct the velocity (smoothly), but the position of the atoms is changing discontinuously.

Since we already know the potential is right, the other thing to check is energy conservation. To test this, initialize the atoms from the eclipsed configuration (maximum torsional energy) with zero velocity and bonds and angles at their minima (say, unit bond length, 90 degree angles, which is the configuration I set up), and then run NVE. The total energy (potential + kinetic) should be conserved, and you should just see the ends of the molecule spin around. But, there will be numerical errors incurred that can lead to energy drift or deviations from this path. For example, with your choice of units, you have many very small or very large numbers (not recommended), and you likely need double precision to get good accuracy.

Regards... Mike

[Nathan Horst]

Thanks Eric and Mike for sorting that out, it seems I got caught up in trying things without fully realizing what I was doing.

As for the new NVE test from the eclipsed configuration, I don't think I am seeing what we would expect here either. Energy is conserved, but there is no transfer of energy from potential to kinetic at all and thus no spinning. I have included the runfile as well as the trajectory, xml dump and energy log for this problem, assuming that I initialized it correctly. I should note that this was done on HOOMD 2.0.0 rather than 1.3.3, purely because it was more convenient for me, but I can run the same test on v1.3.3 if needed.

Thanks,

Nathan

[Michael Howard]

No problem, and no need to run both 1.3.3 and 2.0, they are basically the same. I was a bit stupid, I should have said, "put it *almost* in the eclipsed position". The fully eclipsed position is a maximum of potential energy, and so there is also no force here (gradient is zero). Of course the maximum is unstable, and so any small perturbation will also cause them to start moving -- you could give the beads at either end of the molecule small, equal, and opposite tangential velocities so that the net momentum is zero -- no translation -- but you will start moving around the direction you expect the thing to twist.

Regards... Mike

[Nathan Horst]

Ok, I ran two separate tests, one where i gave the end particles equal but opposite tangential velocities and one where I initialized one particle just outside the eclipsed position. I have included the files, (I apologize for the number of attachments piling up on this thread) and though there is movement, it seems that the particles very much want to stay in the eclipsed position, and not the other way around. As always, it's quite possible that I've overlooked something here, but I would expect the particles to twist away from the eclipsed position?

Thanks,

Nathan

[Eric Jankowski]

If you pretend 0 degrees is 180, 1 degree is 181, and so on, do all
your results make sense?

[Eric Irrgang]

Eric is right. When you measured the potential, you swept phi from -pi to pi and established the minimum at phi=0 (where the 1st particle is at -0.5,1,0 and the 4th particle is at 0.5,1,0). In test_pos.py, you initialize the system very close to phi=0, so I would expect very little motion.

[Michael Howard]

Yup, the definition of torsional angle rears its ugly head again... you can probably do some trig to get your potential shifted by a factor of pi so that +/- pi is the minimum energy.

--Mike

[Nathan Horst]

We modified our potential accordingly, and took the time to run extensive tests with the shifted potential. We now believe that our dihedral angles seem to be implemented properly. Thanks to all the users who helped resolve this issue!

We do wonder if it may be prudent to include a line or two in the documentation mentioning that there seem to be multiple conventions as to what angle corresponds to the trans-state of a dihedral. The references we were following center the trans-state around 0 degrees, rather than the 180 degree convention found in HOOMD. We know now that this can be easily avoided if one is careful enough and knows some basic trig, but if we had known there were multiple conventions at the time of our first implementing the potential, we may have avoided the lost time spent diagnosing the problem. Of course, this is simply a suggestion tailored to other HOOMD-beginners, we will certainly know better in the future.

Thanks again,

Nathan

[Joshua Anderson]

I would be happy to include clarification in the documentation as to the specific angle convention that hoomd utilizes. Since you are familiar with the differences, could you draft the language and/or a figure?

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
Phone: 734-647-8244
http://www-personal.umich.edu/~joaander/

[Nathan Horst]

Hi Josh,

I have included a potential figure, as well as the following section of text shown in red, as a draft of possible additions to the documentation about dihedrals in HOOMD 2.0.0.

This discussion could follow the text already in the documentation stating "where phi is the angle between two sides of the dihedral" within the harmonic dihedral subsection.

Important:

There are multiple conventions pertaining to the dihedral angle (phi) in the literature. HOOMD utilizes the convention shown in the following figure, where vectors are defined from the central particles to the outer particles. These vectors correspond to a stretched state (phi=180 deg) when they are anti-parallel and a compact state (phi=0 deg) when they are parallel.

Of course, there are many ways that this clarification could be done, and any other ideas that are more succinct are certainly welcome. Thanks for your willingness to address this problem,

Nathan Horst