Fix remaining GeometryEngine issues #51

pgrinaway · 2015-12-02T01:57:14Z

As discussed in #49, there are some outstanding issues in the GeometryEngine. Currently:

Bond and angle logp have been converted to logq, so that the exponentiated version of it is the unnormalized probability density.
Deprecated the use of scipy.stats
Corrected the form of sigma in the proposals and logq
Switched to choosing from a discrete set of possible torsions instead of rejection sampling.

What remains:

Calculate reverse torsion probability (next, in this PR)
Cythonize the coordinate conversion (potentially another PR)

pgrinaway · 2015-12-02T02:10:25Z

Ok, I have now implemented the reverse torsion logp.

pgrinaway · 2015-12-02T02:14:05Z

Oops! stray imports from experimentation last night. Removed.

pgrinaway · 2015-12-02T02:16:33Z

I think there are also some unit errors... please stand by...

pgrinaway · 2015-12-02T02:33:08Z

Looks like I am calculating an excessive number of jacobians...

pgrinaway · 2015-12-02T02:35:54Z

Ok. It runs now without error. However, the final energy is pretty bad. I'm going to head home soon and take a look at a pdb of the output to see what is up with that.

jchodera · 2015-12-02T02:47:06Z

Do you mean this is slow in profiling, or you are just computing more than you need?

pgrinaway · 2015-12-02T02:48:10Z

I meant that I had accidentally deleted the calculate_jacobian=False keyword argument, so I was calculating Jacobians where it was unnecessary, and this caused the test to take a ridiculous amount of time.

pgrinaway · 2015-12-02T02:50:19Z

So, computing more than I needed. It's also slow, though.

pgrinaway · 2015-12-02T03:01:21Z

Hm, also it looks like there is an issue with integer division here. 1/17 is unfortunately 0, so I've changed 1 to 1.0 there.

jchodera · 2015-12-02T03:15:43Z

Reviewing PR now!

jchodera · 2015-12-02T03:20:37Z

perses/rjmc/geometry.py

-                logp_r = self._bond_logp(r_proposed, bond, top_proposal.beta)
+                r_proposed = self._propose_bond(bond, beta)
+                bond_k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2
+                sigma_r = units.sqrt(1/beta*bond_k)


This should be 1/(beta*bond_k), which is different from what you have which maps to (1/beta)*bond_k

pgrinaway · 2015-12-07T00:33:46Z

Removing the unit operations did seem to fix the coordinate_tools.py functions.

pgrinaway · 2015-12-07T20:52:46Z

Ok, update:

I have corrected the cartesian_to_internal code.
It looks like the internal_to_cartesian code is where the problem is.

jchodera · 2015-12-07T21:00:17Z

Progress!

jchodera · 2015-12-07T21:28:03Z

Could you check in your fixes when you get a chance?

pgrinaway · 2015-12-07T21:39:12Z

Yeah, sorry--forgot to do that. I still have to clean up the tests, but it looks like the coordinate conversion actually works now. I constructed a minimal OpenMM system of four particles, and added custom bond, angle, and torsion forces to them, where the potential was just the bond, angle, and torsion angle, respectively.

The _cartesian_to_internal code matches what OpenMM calculates, and the _internal_to_cartesian results in OpenMM seeing the same (r, theta, phi) as originally given.
Additionally, the code can reconstruct existing coordinates by going from cartesian to internal and back.

jchodera · 2015-12-07T22:01:10Z

Hooray! So things might work now for actual molecules?

pgrinaway · 2015-12-08T00:01:34Z

Unfortunately, that doesn't seem to be the case. I'm not sure what the deal is here--the coordinates are clearly converted correctly during the proposal, and the potentials seem to be implemented correctly. Is it possible that we can't capture the geometry of a methyl group by proposing only one atom at a time?

jchodera · 2015-12-08T00:15:13Z

Have you tried some totally planar molecules? Like toluene to phenol, where an -H is replaced with an -OH?

pgrinaway · 2015-12-08T01:22:14Z

I did try it, yeah. That is also problematic, though possibly due to ignoring improper torsions in the potential. We should probably, as you said earlier, add these in, but I think maybe a refactor should precede that.

jchodera · 2015-12-08T01:25:09Z

Ethene to ethene-1,2-diol would be free of the improper issue.

jchodera · 2015-12-08T16:49:11Z

We've still got a bug somewhere.

I am running perses/tests/test_geometry_engine.py and printing out the angles from _angle_logq() to other atoms as the torsion is cranked around phi (see below).

It looks like

The constant angle---the one we should have drawn, is 67 degrees, but should be within a few degrees of theta0 = 110 degrees.
Looking through the whole list, there should be some case where all the angles are ~108 degrees, but no such phi is present.

phi = 320.4 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 163.287622503 deg, theta0 = 108.350046723 deg, u_theta = 60.807370
     theta = 91.6814945347 deg, theta0 = 108.350046723 deg, u_theta = 5.597743
phi = 324.0 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 159.969868976 deg, theta0 = 108.350046723 deg, u_theta = 53.684664
     theta = 88.6167925071 deg, theta0 = 108.350046723 deg, u_theta = 7.845392
phi = 327.6 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 156.652525546 deg, theta0 = 108.350046723 deg, u_theta = 47.006300
     theta = 85.584753581 deg, theta0 = 108.350046723 deg, u_theta = 10.441520
phi = 331.2 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 153.337070837 deg, theta0 = 108.350046723 deg, u_theta = 40.774793
     theta = 82.5891419396 deg, theta0 = 108.350046723 deg, u_theta = 13.370248
phi = 334.8 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 150.024590927 deg, theta0 = 108.350046723 deg, u_theta = 34.991209
     theta = 79.6342293237 deg, theta0 = 108.350046723 deg, u_theta = 16.613441
phi = 338.4 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 146.715991429 deg, theta0 = 108.350046723 deg, u_theta = 29.655758
     theta = 76.7248732788 deg, theta0 = 108.350046723 deg, u_theta = 20.150372
phi = 342.0 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 143.412097073 deg, theta0 = 108.350046723 deg, u_theta = 24.768052
     theta = 73.8666064805 deg, theta0 = 108.350046723 deg, u_theta = 23.957330
phi = 345.6 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 140.113704504 deg, theta0 = 108.350046723 deg, u_theta = 20.327232
     theta = 71.0657376998 deg, theta0 = 108.350046723 deg, u_theta = 28.007183
phi = 349.2 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 136.821613765 deg, theta0 = 108.350046723 deg, u_theta = 16.332023
     theta = 68.3294642881 deg, theta0 = 108.350046723 deg, u_theta = 32.268891
phi = 352.8 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 133.536649488 deg, theta0 = 108.350046723 deg, u_theta = 12.780752
     theta = 65.6659948992 deg, theta0 = 108.350046723 deg, u_theta = 36.706968
phi = 356.4 deg
     theta = 67.1671869485 deg, theta0 = 110.050047451 deg, u_theta = 43.570751
     theta = 130.259677014 deg, theta0 = 108.350046723 deg, u_theta = 9.671355
     theta = 63.0846792812 deg, theta0 = 108.350046723 deg, u_theta = 41.280917

jchodera · 2015-12-08T16:55:19Z

Here's the first proposal, which shows that angle (1,3,2) is proposed as 114 degrees but then measured to be 65 degrees in the subsequent energy calculations. This is our primary problem.

PROPOSED ANGLE
    1     3     2 : theta = 114.442511463 deg, theta0 = 110.630047471 deg
sigma_theta = 3.93456639833 deg, theta_proposed = 114.442511463 deg
     theta = 114.442511463 deg, theta0 = 110.630047471 deg, u_theta = 0.469448
phi = 0.0 deg
    1     3     2 :       theta = 65.5574885365 deg, theta0 = 110.630047471 deg, u_theta = 65.614629
    1     3    12 :       theta = 122.094925262 deg, theta0 = 110.050047451 deg, u_theta = 3.437418
    1     3    13 :       theta = 128.08734777 deg, theta0 = 110.050047451 deg, u_theta = 7.708513
phi = 3.6 deg
    1     3     2 :       theta = 65.5574885365 deg, theta0 = 110.630047471 deg, u_theta = 65.614629
    1     3    12 :       theta = 125.33750595 deg, theta0 = 110.050047451 deg, u_theta = 5.537301
    1     3    13 :       theta = 124.846376294 deg, theta0 = 110.050047451 deg, u_theta = 5.187230
phi = 7.2 deg
    1     3     2 :       theta = 65.5574885365 deg, theta0 = 110.630047471 deg, u_theta = 65.614629
    1     3    12 :       theta = 128.589171328 deg, theta0 = 110.050047451 deg, u_theta = 8.143404
    1     3    13 :       theta = 121.615624212 deg, theta0 = 110.050047451 deg, u_theta = 3.169291

jchodera · 2015-12-08T16:57:51Z

perses/rjmc/geometry.py

@@ -650,10 +646,15 @@ def _calculate_angle(self, atom_position, bond_atom_position, angle_atom_positio
        b = angle_atom_position - bond_atom_position
        a_u = a / np.linalg.norm(a)
        b_u = b / np.linalg.norm(b)
-        theta = np.arccos(np.dot(a_u, b_u))
+        cos_theta = np.dot(a_u, b_u)


This is the bug. Because of the way you define a and b vectors, this should be

cos_theta = np.dot(-a_u, b_u)

pgrinaway · 2015-12-08T19:22:17Z

Looks like some unit thing is a problem in one of the tests. Working on that.

jchodera · 2015-12-08T19:26:30Z

perses/tests/test_geometry_engine.py

+        bond_position = units.Quantity(np.array([ 0.0765,  0.1  ,  -0.4005]), unit=units.nanometers)
+        angle_position = units.Quantity(np.array([ 0.0829 , 0.0952 ,-0.2479]) ,unit=units.nanometers)
+        torsion_position = units.Quantity(np.array([-0.057 ,  0.0951 ,-0.1863] ) ,unit=units.nanometers)
+        rtp, detJ = geometry_engine._cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position)


Your documentation for _cartesian_to_internal says:

Returns without units!

but then you feed the unitless r into _internal_to_cartesian which immediately removes units from r, assuming it is unit-bearing.

We will need some way to be more consistent with units in a refactoring PR.

We will need some way to be more consistent with units in a refactoring PR.

Agreed. The error in the current build is due to a different issue with units, though this is also a problem

pgrinaway · 2015-12-08T19:42:28Z

Ok, right now it looks like the only problem in the tests is the test_rotation_matrix. I'll rewrite that test then commit, and merge this PR. Next PR will clean things up and add the rest of the unit tests.

I should clarify that the problem is an exception caused by an incorrectly written test.

pgrinaway · 2015-12-08T19:44:26Z

Also need to restore the Jacobian calculation--leaving this here as a note to myself.

pgrinaway · 2015-12-08T20:52:22Z

Passes! Going to merge.

Fix remaining GeometryEngine issues

pgrinaway added 10 commits December 1, 2015 14:09

replace rejection sampling with selection from precomputed values

a6b5636

fixing probability calculations

8dbb9c1

moved old geometry engine to the attic

b11e958

add 1/2 prefactor

a1b3e69

calculate logZ separately

24ed133

fixed bond and angle proposals

1bb11e4

make sure beta is in the right units

48ffa76

fixed angle thing

38ac293

removed unused method

1344f99

calculate reverse torsion_logp

d2adbb6

removed unnecessary imports

bd08b30

pgrinaway added 3 commits December 1, 2015 21:18

fixed some unit stuff

25aedf7

remove loop

9346fff

remove another loop

52aa9d2

make jacobian optional again

cf92673

prevent integer division

ec3e218

jchodera reviewed Dec 2, 2015
View reviewed changes

fixes

0ba012b

jchodera reviewed Dec 8, 2015
View reviewed changes

pgrinaway added 3 commits December 8, 2015 13:44

fix angles

806adc0

test angle computation

d956d88

remove plotting

5ca5c9f

jchodera reviewed Dec 8, 2015
View reviewed changes

temprorary unit fix

cba5b45

pgrinaway added 3 commits December 8, 2015 15:39

change name of test to make run

390efd1

calculate determinant of jacobian

d465db1

fix unit

6d5788e

pgrinaway changed the title ~~[WIP] Fix remaining GeometryEngine issues~~ Fix remaining GeometryEngine issues Dec 8, 2015

pgrinaway added a commit that referenced this pull request Dec 8, 2015

Merge pull request #51 from choderalab/fix-geometry

d540b53

Fix remaining GeometryEngine issues

pgrinaway merged commit d540b53 into master Dec 8, 2015

pgrinaway mentioned this pull request Feb 12, 2016

Bond constraints #50

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix remaining GeometryEngine issues #51

Fix remaining GeometryEngine issues #51

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

jchodera commented Dec 2, 2015 via email

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

jchodera commented Dec 2, 2015

jchodera Dec 2, 2015

pgrinaway commented Dec 7, 2015

pgrinaway commented Dec 7, 2015

jchodera commented Dec 7, 2015

jchodera commented Dec 7, 2015

pgrinaway commented Dec 7, 2015

jchodera commented Dec 7, 2015

pgrinaway commented Dec 8, 2015

jchodera commented Dec 8, 2015 via email

pgrinaway commented Dec 8, 2015

jchodera commented Dec 8, 2015 via email

jchodera commented Dec 8, 2015

jchodera commented Dec 8, 2015

jchodera Dec 8, 2015

pgrinaway commented Dec 8, 2015

jchodera Dec 8, 2015

pgrinaway Dec 8, 2015

pgrinaway commented Dec 8, 2015

pgrinaway commented Dec 8, 2015

pgrinaway commented Dec 8, 2015

Fix remaining GeometryEngine issues #51

Fix remaining GeometryEngine issues #51

Conversation

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

jchodera commented Dec 2, 2015 via email

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

pgrinaway commented Dec 2, 2015

jchodera commented Dec 2, 2015

jchodera Dec 2, 2015

Choose a reason for hiding this comment

pgrinaway commented Dec 7, 2015

pgrinaway commented Dec 7, 2015

jchodera commented Dec 7, 2015

jchodera commented Dec 7, 2015

pgrinaway commented Dec 7, 2015

jchodera commented Dec 7, 2015

pgrinaway commented Dec 8, 2015

jchodera commented Dec 8, 2015 via email

pgrinaway commented Dec 8, 2015

jchodera commented Dec 8, 2015 via email

jchodera commented Dec 8, 2015

jchodera commented Dec 8, 2015

jchodera Dec 8, 2015

Choose a reason for hiding this comment

pgrinaway commented Dec 8, 2015

jchodera Dec 8, 2015

Choose a reason for hiding this comment

pgrinaway Dec 8, 2015

Choose a reason for hiding this comment

pgrinaway commented Dec 8, 2015

pgrinaway commented Dec 8, 2015

pgrinaway commented Dec 8, 2015