Skip to content

Conversation

@jchodera
Copy link
Member

@jchodera jchodera commented Jul 10, 2022

Bugfixes:

  • Fix critical bug introduced in Remove debug prints. #198 that resulted in torsions being omitted from template-generated systems

Minor improvements:

  • Change how we measure relative deviation in vectors to avoid underflow.

@jchodera
Copy link
Member Author

Hm, I'm not sure why this is happening. There is a significant energy discrepancy for '[H]C(=O)[H]' in openmmforcefields/tests/test_template_generators.py::TestSMIRNOFFTemplateGenerator::test_energies, but none of the individual components has a significant difference:

Force components:
component                      relative deviation
NonbondedForce                       0.0000000000
HarmonicBondForce                    0.0000000000
HarmonicAngleForce                   0.0000000000
TOTAL                                0.0000858405

@jchodera
Copy link
Member Author

The System created by the SMIRNOFFTemplateGenerator appears to be missing impropers.
Investigating further.

@jchodera
Copy link
Member Author

jchodera commented Jul 10, 2022

OK, what is currently happening is that the OpenFF toolkit generates an OpenMM System containing a SMIRNOFF improper:

<?xml version="1.0" ?>
<System openmmVersion="7.7" type="System" version="1">
	<PeriodicBoxVectors>
		<A x="2" y="0" z="0"/>
		<B x="0" y="2" z="0"/>
		<C x="0" y="0" z="2"/>
	</PeriodicBoxVectors>
	<Particles>
		<Particle mass="12.01078"/>
		<Particle mass="15.99943"/>
		<Particle mass="1.007947"/>
		<Particle mass="1.007947"/>
	</Particles>
	<Constraints/>
	<Forces>
		<Force forceGroup="0" name="PeriodicTorsionForce" type="PeriodicTorsionForce" usesPeriodic="0" version="2">
			<Torsions>
				<Torsion k="1.5341333333333336" p1="0" p2="1" p3="2" p4="3" periodicity="2" phase="3.141592653589793"/>
				<Torsion k="1.5341333333333336" p1="0" p2="2" p3="3" p4="1" periodicity="2" phase="3.141592653589793"/>
				<Torsion k="1.5341333333333336" p1="0" p2="3" p3="1" p4="2" periodicity="2" phase="3.141592653589793"/>
			</Torsions>
		</Force>
		<Force alpha="0" cutoff=".9" dispersionCorrection="1" ewaldTolerance=".0005" exceptionsUsePeriodic="0" forceGroup="0" includeDirectSpace="1" ljAlpha="0" ljnx="0" ljny="0" ljnz="0" method="0" name="NonbondedForce" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="78.3" switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" version="4">
			<GlobalParameters/>
			<ParticleOffsets/>
			<ExceptionOffsets/>
			<Particles>
				<Particle eps=".359824" q=".5632600113749504" sig=".3399669508423535"/>
				<Particle eps=".87864" q="-.5147399976849556" sig=".2959921901149463"/>
				<Particle eps=".06276" q="-.024260006844997406" sig=".2510552587719476"/>
				<Particle eps=".06276" q="-.024260006844997406" sig=".2510552587719476"/>
			</Particles>
			<Exceptions>
				<Exception eps="0" p1="0" p2="1" q="0" sig="1"/>
				<Exception eps="0" p1="0" p2="2" q="0" sig="1"/>
				<Exception eps="0" p1="1" p2="2" q="0" sig="1"/>
				<Exception eps="0" p1="0" p2="3" q="0" sig="1"/>
				<Exception eps="0" p1="1" p2="3" q="0" sig="1"/>
				<Exception eps="0" p1="2" p2="3" q="0" sig="1"/>
			</Exceptions>
		</Force>
		<Force forceGroup="0" name="HarmonicBondForce" type="HarmonicBondForce" usesPeriodic="0" version="2">
			<Bonds>
				<Bond d=".12290000000000001" k="476976.00000000006" p1="0" p2="1"/>
				<Bond d=".10800000000000001" k="307105.60000000003" p1="0" p2="2"/>
				<Bond d=".10800000000000001" k="307105.60000000003" p1="0" p2="3"/>
			</Bonds>
		</Force>
		<Force forceGroup="0" name="HarmonicAngleForce" type="HarmonicAngleForce" usesPeriodic="0" version="2">
			<Angles>
				<Angle a="2.0943951023931953" k="418.40000000000003" p1="1" p2="0" p3="2"/>
				<Angle a="2.0943951023931953" k="418.40000000000003" p1="1" p2="0" p3="3"/>
				<Angle a="2.0943951023931953" k="292.88" p1="2" p2="0" p3="3"/>
			</Angles>
		</Force>
	</Forces>
</System>

But the torsion processing part of convert_system_to_ffxml() is not generating torsion terms for the improper:

<ForceField>
  <AtomTypes>
    <Type name="[H]C(=O)[H]$C1#0" element="C" mass="12.01078" class="[H]C(=O)[H]$C1#0"/>
    <Type name="[H]C(=O)[H]$O1#1" element="O" mass="15.99943" class="[H]C(=O)[H]$O1#1"/>
    <Type name="[H]C(=O)[H]$H1#2" element="H" mass="1.007947" class="[H]C(=O)[H]$H1#2"/>
    <Type name="[H]C(=O)[H]$H2#3" element="H" mass="1.007947" class="[H]C(=O)[H]$H2#3"/>
  </AtomTypes>
  <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
    <UseAttributeFromResidue name="charge"/>
    <Atom sigma="0.3399669508423535" epsilon="0.359824" class="[H]C(=O)[H]$C1#0"/>
    <Atom sigma="0.2959921901149463" epsilon="0.87864" class="[H]C(=O)[H]$O1#1"/>
    <Atom sigma="0.2510552587719476" epsilon="0.06276" class="[H]C(=O)[H]$H1#2"/>
    <Atom sigma="0.2510552587719476" epsilon="0.06276" class="[H]C(=O)[H]$H2#3"/>
  </NonbondedForce>
  <HarmonicBondForce>
    <Bond class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$O1#1" length="0.12290000000000001" k="476976.00000000006"/>
    <Bond class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$H1#2" length="0.10800000000000001" k="307105.60000000003"/>
    <Bond class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$H2#3" length="0.10800000000000001" k="307105.60000000003"/>
  </HarmonicBondForce>
  <HarmonicAngleForce>
    <Angle class1="[H]C(=O)[H]$O1#1" class2="[H]C(=O)[H]$C1#0" class3="[H]C(=O)[H]$H1#2" angle="2.0943951023931953" k="418.40000000000003"/>
    <Angle class1="[H]C(=O)[H]$O1#1" class2="[H]C(=O)[H]$C1#0" class3="[H]C(=O)[H]$H2#3" angle="2.0943951023931953" k="418.40000000000003"/>
    <Angle class1="[H]C(=O)[H]$H1#2" class2="[H]C(=O)[H]$C1#0" class3="[H]C(=O)[H]$H2#3" angle="2.0943951023931953" k="292.88"/>
  </HarmonicAngleForce>
  <Residues>
    <Residue name="[H]C(=O)[H]">
      <Atom name="C1" type="[H]C(=O)[H]$C1#0" charge="0.5632600113749504"/>
      <Atom name="O1" type="[H]C(=O)[H]$O1#1" charge="-0.5147399976849556"/>
      <Atom name="H1" type="[H]C(=O)[H]$H1#2" charge="-0.024260006844997406"/>
      <Atom name="H2" type="[H]C(=O)[H]$H2#3" charge="-0.024260006844997406"/>
      <Bond atomName1="C1" atomName2="O1"/>
      <Bond atomName1="C1" atomName2="H1"/>
      <Bond atomName1="C1" atomName2="H2"/>
    </Residue>
  </Residues>
</ForceField>

@jchodera
Copy link
Member Author

OK, it looks like some critical code was removed in #198, a PR that was intended to eliminate extra debug printing but instead may have eliminated critical code needed to correctly convert all torsions for SMIRNOFF and Espaloma. Specifically, this change was one of the important deletions.

@jchodera jchodera changed the title Fix test failures Critical bugfix: Torsions (propers and impropers) were not being added to template-generated systems Jul 10, 2022
@jchodera jchodera added this to the 0.11.1 milestone Jul 10, 2022
@jchodera jchodera self-assigned this Jul 10, 2022
@jchodera
Copy link
Member Author

Fortunately, it looks like this bug didn't make it into a release!

@jchodera
Copy link
Member Author

@yuanqing-wang: I'm running into a type error with the tests parameterizing acetone with espaloma 0.2.2:

2022-07-10T22:14:21.2517660Z openmmforcefields/tests/test_template_generators.py:136: 
2022-07-10T22:14:21.2518497Z _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
2022-07-10T22:14:21.2519397Z openmmforcefields/tests/test_template_generators.py:131: in test_add_molecules
2022-07-10T22:14:21.2520339Z     system = forcefield.createSystem(openmm_topology, nonbondedMethod=NoCutoff)
2022-07-10T22:14:21.2521471Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmm/app/forcefield.py:1212: in createSystem
2022-07-10T22:14:21.2522587Z     templateForResidue = self._matchAllResiduesToTemplates(data, topology, residueTemplates, ignoreExternalBonds)
2022-07-10T22:14:21.2523839Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmm/app/forcefield.py:1417: in _matchAllResiduesToTemplates
2022-07-10T22:14:21.2524908Z     if generator(self, res):
2022-07-10T22:14:21.2526141Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmmforcefields-0+untagged.1.g7f40898-py3.10.egg/openmmforcefields/generators/template_generators.py:304: in generator
2022-07-10T22:14:21.2534583Z     ffxml_contents = self.generate_residue_template(molecule)
2022-07-10T22:14:21.2536082Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/openmmforcefields-0+untagged.1.g7f40898-py3.10.egg/openmmforcefields/generators/template_generators.py:1629: in generate_residue_template
2022-07-10T22:14:21.2537140Z     molecule_graph = esp.Graph(molecule)
2022-07-10T22:14:21.2538121Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/espaloma/graphs/graph.py:62: in __init__
2022-07-10T22:14:21.2539044Z     heterograph = self.get_heterograph_from_graph_and_mol(
2022-07-10T22:14:21.2540284Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/espaloma/graphs/graph.py:126: in get_heterograph_from_graph_and_mol
2022-07-10T22:14:21.2541312Z     heterograph = esp.graphs.utils.read_heterogeneous_graph.from_homogeneous_and_mol(
2022-07-10T22:14:21.2542775Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/espaloma/graphs/utils/read_heterogeneous_graph.py:272: in from_homogeneous_and_mol
2022-07-10T22:14:21.2543789Z     hg = dgl.heterograph({key: list(value) for key, value in hg.items()})
2022-07-10T22:14:21.2544827Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/dgl/convert.py:350: in heterograph
2022-07-10T22:14:21.2545792Z     (sparse_fmt, arrays), urange, vrange = utils.graphdata2tensors(
2022-07-10T22:14:21.2546856Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/dgl/utils/data.py:179: in graphdata2tensors
2022-07-10T22:14:21.2547732Z     src, dst = elist2tensor(data, idtype)
2022-07-10T22:14:21.2548688Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/dgl/utils/data.py:32: in elist2tensor
2022-07-10T22:14:21.2549584Z     return F.tensor(u, idtype), F.tensor(v, idtype)
2022-07-10T22:14:21.2550381Z _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
2022-07-10T22:14:21.2550830Z 
2022-07-10T22:14:21.2551192Z data = [0.0, 1.0, 2.0, 3.0], dtype = torch.int64
2022-07-10T22:14:21.2551608Z 
2022-07-10T22:14:21.2551951Z     def tensor(data, dtype=None):
2022-07-10T22:14:21.2552696Z         if isinstance(data, numbers.Number):
2022-07-10T22:14:21.2553407Z             data = [data]
2022-07-10T22:14:21.2554221Z         if isinstance(data, list) and len(data) > 0 and isinstance(data[0], th.Tensor):
2022-07-10T22:14:21.2555140Z             # prevent GPU->CPU->GPU copies
2022-07-10T22:14:21.2555870Z             if data[0].ndim == 0:
2022-07-10T22:14:21.2556613Z                 # zero dimenion scalar tensors
2022-07-10T22:14:21.2557372Z                 return th.stack(data)
2022-07-10T22:14:21.2558123Z         if isinstance(data, th.Tensor):
2022-07-10T22:14:21.2558928Z             return th.as_tensor(data, dtype=dtype, device=data.device)
2022-07-10T22:14:21.2559660Z         else:
2022-07-10T22:14:21.2560380Z >           return th.as_tensor(data, dtype=dtype)
2022-07-10T22:14:21.2561307Z E           TypeError: 'numpy.float64' object cannot be interpreted as an integer
2022-07-10T22:14:21.2561790Z 
2022-07-10T22:14:21.2562375Z /usr/share/miniconda/envs/test/lib/python3.10/site-packages/dgl/backend/pytorch/tensor.py:46: TypeError

Could this be due to a version issue? Here's how we are installing espaloma and its dependencies in the unit tests:
https://github.com/openmm/openmmforcefields/blob/master/devtools/conda-envs/test_env.yaml#L41-L50

@yuanqing-wang
Copy link

yuanqing-wang commented Jul 11, 2022 via email

@yuanqing-wang
Copy link

@jchodera could you list the pytorch and dgl versions used to produce that error?

@jchodera
Copy link
Member Author

Here are the versions that were installed in the conda environment:

    + pytorch                              1.11.0  cpu_py310h75c9ab6_2      conda-forge/linux-64       47 MB
    + dgl                                   0.8.2  py310_0                  dglteam/linux-64            5 MB

@jchodera jchodera requested a review from mikemhenry July 11, 2022 03:57
@jchodera jchodera changed the title Critical bugfix: Torsions (propers and impropers) were not being added to template-generated systems Fix accidental elimination of torsions after 0.11.0 release Jul 11, 2022
@jchodera
Copy link
Member Author

@yuanqing-wang : Do you think we need to pin the DGL version for now?

@yuanqing-wang
Copy link

yuanqing-wang commented Jul 11, 2022 via email

@yuanqing-wang
Copy link

@jchodera fixed. could you try with the current espaloma master branch

@yuanqing-wang
Copy link

I accidentally had a float as graph index which was tolerated by older dgl but not newer versions.

@jchodera
Copy link
Member Author

@yuanqing-wang : Thanks! It looks like this works, though there is a segfault on osx.
Can you cut a new release with this bugfix?

We can debug the osx failure separately.

@jchodera jchodera marked this pull request as ready for review July 12, 2022 05:24
@jchodera jchodera merged commit 815dfd1 into master Jul 12, 2022
@jchodera jchodera deleted the fix-tests branch July 12, 2022 05:26
@yuanqing-wang
Copy link

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants