Snapshot converter, hoomd simulation module #622

ahy3nz · 2019-10-16T22:53:29Z

PR Summary:

Attempt at #619 and #570 . Given a parameterized structure, create a hoomd simulation, with all topological and force field information successfully preserved.

import mbuild as mb
from foyer import Forcefield
from mbuild.examples import Ethane
from mbuild.formats.hoomd_simulation import create_hoomd_simulation

cmpd = Ethane()
ff = Forcefield(name='oplsaa')
struc = ff.apply(cmpd)
# like struc.createSystem(), but now for hoomd and not openmm
create_hoomd_simulation(struc) # this is where all the hoomd information is set up

import hoomd
import hoomd.md
import hoomd.group
# now we just have to specify the simulation run parameters
all = hoomd.group.all()
hoomd.md.integrate.mode_standard(dt=0.001)
integrator = hoomd.md.integrate.nvt(group=all, kT=1, tau=1)
hoomd.dump.dcd('traj.dcd', period=10, group=all, phase=0, overwrite=True,
    unwrap_full=False)

hoomd.run(1e3)

This code works all the way until the you hit hoomd.run when the simulation crashes. But it looks like we're successfully parametrizing bond and nonbonded interactions:

LJ and charge interactions based on atomtypes, nbfixes, and a mixing rule check
14 interactions based on struc.adjusts
harmonic bonds, harmonic angles
periodic torsions, rb torsions

An intermediate step was generating hoomd.data.Snapshot from parmed.Structure

PR Checklist

Includes appropriate unit test(s)
Appropriate docstring(s) are added/updated
Code is (approximately) PEP8 compliant
Issue(s) raised/addressed?

lgtm-com · 2019-10-16T23:11:10Z

This pull request introduces 8 alerts when merging ad0353a into ed96251 - view on LGTM.com

new alerts:

8 for Unused local variable

codecov · 2019-10-16T23:18:47Z

Codecov Report

Merging #622 into master will decrease coverage by 8.25%.
The diff coverage is 0.51%.

@@            Coverage Diff             @@
##           master     #622      +/-   ##
==========================================
- Coverage   90.81%   82.56%   -8.26%     
==========================================
  Files          50       52       +2     
  Lines        3898     4290     +392     
==========================================
+ Hits         3540     3542       +2     
- Misses        358      748     +390

Impacted Files	Coverage Δ
mbuild/formats/gsdwriter.py	`95.86% <ø> (ø)`	⬆️
mbuild/formats/hoomd_snapshot.py	`0% <0%> (ø)`
mbuild/formats/hoomd_simulation.py	`0% <0%> (ø)`
mbuild/utils/io.py	`77.61% <50%> (-1.76%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 68ecd75...643a2a4. Read the comment docs.

justinGilmer

Just a typo fix and changing the hoomd imports to use _import.

I can try this on a few UA systems later when you think its ready to test.

justinGilmer · 2019-10-17T14:59:39Z

mbuild/formats/hoomd_simulation.py

+              ref_energy=1.0, mixing_rule='lorentz', r_cut=1.2, 
+              snapshot_kwargs={}, 
+              pppm_kwargs={'Nx':1, 'Ny':1, 'Nz':1, 'order':4}):
+    import hoomd


I guess we should probably do an _import here, right? Since hoomd is not a requirement.

justinGilmer · 2019-10-17T15:05:45Z

mbuild/formats/hoomd_simulation.py

+    return harmonic_angle
+
+def _init_hoomd_dihedrals(structure, ref_energy=1.0):
+    # Identify the uniqeu dihedral types before setting


Suggested change

# Identify the uniqeu dihedral types before setting

# Identify the unique dihedral types before setting

justinGilmer · 2019-10-17T15:06:54Z

mbuild/formats/hoomd_snapshot.py

+              ref_energy=1.0, rigid_bodies=None, shift_coords=True,
+              write_special_pairs=True):
+    """Convert parmed.Structure to hoomd.data.Snapshot"""
+    import hoomd


same here for the _import

…rmed structure

… and structure

…ization, docuemtnatoin

…ng FF translations, passing bad inputs

ahy3nz · 2019-10-17T17:26:17Z

Barring some other issues with dependency syntax breaking foyer (netwworkx), I think this PR is ready to be looked at

Converting compounds and structures to hoomd.Snapshot
Converting structures to hoomd.SimulationContext
Unit tests for input-error-checking, proper data conversions to snapshots, proper data conversions to simulations

import mbuild as mb
from foyer import Forcefield
from mbuild.examples import Ethane

compound = Ethane()
mixture_box = mb.fill_box(compound, n_compounds=10, box=[10,10,10])

ff = Forcefield(name='oplsaa')
struc = ff.apply(mixture_box)
struc.save('stuff.mol2', overwrite=True)

from mbuild.formats.hoomd_simulation import create_hoomd_simulation
create_hoomd_simulation(struc, ref_distance=10, r_cut=1.2, snapshot_kwargs={'shift_coords':True})

import hoomd
import hoomd.md
import hoomd.group
all = hoomd.group.all()
hoomd.md.integrate.mode_standard(dt=0.0001)
integrator = hoomd.md.integrate.nvt(group=all, kT=2.5, tau=1)
hoomd.dump.dcd('query.dcd', period=10, group=all, phase=0, overwrite=True,
    unwrap_full=False)

hoomd.run(1e3)

lgtm-com · 2019-10-17T17:34:36Z

This pull request introduces 8 alerts when merging c924855 into ed96251 - view on LGTM.com

new alerts:

8 for Unused local variable

ahy3nz · 2019-10-17T21:50:14Z

import mbuild as mb
import parmed as pmd
import foyer
from mbuild.examples import Ethane
import hoomd

eth = Ethane()
cmpd = mb.fill_box(eth, n_compounds=10, box=[10,10,10])
ff = foyer.Forcefield(name='oplsaa')
structure = ff.apply(cmpd)
for atom in structure.atoms:
    atom.charge = 0

from mbuild.formats.hoomd_simulation import create_hoomd_simulation
create_hoomd_simulation(structure, ref_distance=10)

all = hoomd.group.all()
hoomd.md.integrate.mode_standard(dt=0.0000001)
hoomd.md.integrate.nve(all)
hoomd.run(1)

total_energy = 0
for force in hoomd.context.current.forces:
    force_energy = force.get_energy(all)
    print((force, force_energy))
    total_energy += force_energy

print("Hoomd energy:{}".format(total_energy))

import simtk.openmm as openmm
import simtk.unit as unit
omm_system = structure.createSystem()

integrator = openmm.VerletIntegrator(1.0)
omm_context = openmm.Context(omm_system, integrator)
omm_context.setPositions(structure.positions)
omm_state = omm_context.getState(getEnergy=True)
omm_energy = omm_state.getPotentialEnergy().in_units_of(unit.kilocalories_per_mole)
print("Openmm energy:{}".format(omm_energy))

Hoomd energy:352.7569344857175
Openmm energy:352.75510711378166 kcal/mol

If charges are included, I get a lot of malloc errors

lgtm-com · 2019-10-17T22:01:51Z

This pull request introduces 8 alerts when merging 7ebe624 into ed96251 - view on LGTM.com

new alerts:

8 for Unused local variable

ahy3nz · 2019-10-18T23:46:43Z

For a box of two chargeless 30-anes (alkane with length 30):

Hoomd energy:66036.05035805175
Openmm energy:66052.78125 kJ/mol
E_Group         Hoomd     Openmm     % diff
Bond       1104.877   1104.853   0.002
Angle      35956.915  35952.957  0.011
Dihedral   338.851    402.129    -15.736
Nonbond    28635.408  28592.844  0.149

From this, it looks like bonds, angles, and LJs are converted correctly.
Someone can double check my scratchwork, but the dihedral discrepancy is unknown. Even though openmm uses RB coeffs and hoomd uses OPLS coeffs, converting the openmm RB coeffs to OPLS yields the same coeffs. Maybe it's very slight numerical differences that propagate out to lots of dihedrals? Maybe it's how each engine is computing these energies? As it seems to me, even though the energies evaluated differently, I think we did our job in correctly processing dihedral coefficients.

Engine	F1	F2	F3
openmm	0	0	1.2552
openmm	5.4392	-0.2092	0.8368
hoomd	0	0	1.2552
hoomd	5.4392	-0.2092	0.8368

lgtm-com · 2019-10-18T23:46:45Z

This pull request introduces 8 alerts when merging ab79fdc into d6ad4db - view on LGTM.com

new alerts:

8 for Unused local variable

uppittu11 · 2019-10-21T14:49:16Z

@ahy3nz Can you put a script to replicate the 30-anes results in a gist or something so we don't have to rewrite it? The 10 ethane simulation in your previous comment seemed to have agreeing total energies, but was the dihedral energy off in that one too?

ahy3nz · 2019-10-21T15:01:02Z

energy_comparison.zip

ahy3nz · 2019-10-31T22:46:41Z

For time sensitivity and calling it a day with this PR, I think this is ready to go. We're doing well translating bonded-parameters (with hoomd 2.8.0, and we also check for hoomd version). There might be some issues with nonbonded interactions, but diagnosing those issues might take longer. So for now I think this is ready for review before we do rigorous energy comparisons that could take a while

If anyone's interested, try adapting this code:

import mbuild as mb
from foyer import Forcefield
from mbuild.examples import Ethane

compound = Ethane()
mixture_box = mb.fill_box(compound, n_compounds=10, box=[10,10,10])

ff = Forcefield(name='oplsaa')
struc = ff.apply(mixture_box)
struc.save('stuff.mol2', overwrite=True)

from mbuild.formats.hoomd_simulation import create_hoomd_simulation
create_hoomd_simulation(struc, ref_distance=10, r_cut=1.2, snapshot_kwargs={'shift_coords':True})

import hoomd
import hoomd.md
import hoomd.group
all = hoomd.group.all()
hoomd.md.integrate.mode_standard(dt=0.0001)
integrator = hoomd.md.integrate.nvt(group=all, kT=2.5, tau=1)
hoomd.dump.dcd('query.dcd', period=10, group=all, phase=0, overwrite=True,
    unwrap_full=False)

hoomd.run(1e3)

ahy3nz · 2019-10-31T22:47:15Z

@mosdef-hub/mosdef-contributors

mikemhenry · 2019-11-01T15:09:45Z

I will test this PR with one of our systems and if it works I will then review it! Thanks for all the hard work @ahy3nz I'm really excited for your energy validation project.

…into hoomd_integration

mikemhenry · 2019-11-05T00:16:11Z

Our lab has used this for atomistic and coarse-grain systems without issue, ~~I think the only thing missing is a way to get the values post autoscale~~ Never mind! It is returned.

mattwthompson · 2019-11-05T00:30:40Z

Could one of the UMich folks have a look at this? I am less familiar with the internals of HOOMD so I am hesitant to review it myself.

joaander

Thanks! This looks great.

I don't see constrain.distance. Is this not supported, or are you not targeting any force fields that use hard distance constraints?

joaander · 2019-11-05T12:53:19Z

mbuild/formats/hoomd_simulation.py

+
+def _check_hoomd_version():
+    version = hoomd.__version__
+    version_numbers = version.split('.')


Consider adding a check and error for major version >= 3.

…md snapshot

ahy3nz · 2019-11-05T20:26:31Z

With the current foyer Forcefield schema and parmed.Structure, it appears information about distance constraints are not codified or present in the foyer Forcefield or resultant parmed.Structure. Something like hoomd.md.constrain.distance, at this current state of Foyer and ParmEd, would require a user to specify these constraints.

As for constrain.rigid, I'm working through to see how easy it would be to specify rigid bodies from an existing mbuild compound/parmed structure

mikemhenry

Looks good! I wouldn't worry about the rigid body code for this PR.

ahy3nz · 2019-11-06T00:34:08Z

Okay! Will deal with rigid bodies in another PR - I think this is good to go if someone wants to try putting together some mbuild + foyer + hoomd systems

mikemhenry · 2019-11-06T01:10:44Z

It works pretty well: https://github.com/cmelab/CG-Tutorial

Add HOOMD snapshot and simulation converters

Snapshot converter, hoomd simulation module

ad0353a

justinGilmer requested changes Oct 17, 2019

View reviewed changes

ahy3nz added 12 commits October 17, 2019 10:56

Use safer imports, fix typo in docstring

312fc17

Safer imports, update neighborlist if 1-4 adjustments are found in pa…

7f5ea37

…rmed structure

Add has_hoomd for imports

349e2ae

Move hoomdxml tests into test_hoomd, test hoomd snapshot for compound…

8cfcb8b

… and structure

Remove testhoomdxml since it's now in test_hoomd

0f238b3

Fix imports, input type checking

df8f79c

Modify gsd docstring

f324137

More documentation

559a62a

Eliminate rcut conversions, leave in dimensionless units. Input sanit…

55ea0f3

…ization, docuemtnatoin

Extra quotes

7638f57

Typo in instance checking

85c2a72

More exhaustive tests for hoomdsnapshots and hoomdsimulations, checki…

c924855

…ng FF translations, passing bad inputs

ahy3nz marked this pull request as ready for review October 17, 2019 17:26

Don't initialize electrostatics if no charged group

7ebe624

ahy3nz added 2 commits October 18, 2019 16:01

Merge remote-tracking branch 'upstream/master' into hoomd_integration

12afdbd

More docstring

ab79fdc

Adjust typos in code

3f2584c

Check for non-parametrized bonds as constraints

0bdb7fb

ahy3nz added 6 commits October 30, 2019 12:26

Chekc hoomd version if using periodic/charmm dihedrals

02fc434

make list of hoomd objects and return them

6a316ce

Update mixing rules

a04b0b1

Don't initialize hoomd context if already initialized

d73612b

Autoscale for detecting reference units

4ad8aa7

Typo in arg

0210024

mikemhenry self-requested a review November 1, 2019 15:08

added import that was missing

09c6ad1

ahy3nz mentioned this pull request Nov 4, 2019

UA custom element example mosdef-hub/mosdef-workflows#10

Merged

ahy3nz added 2 commits November 4, 2019 09:52

Merge remote-tracking branch 'upstream/master' into hoomd_integration

86504c8

Merge branch 'hoomd_integration' of https://github.com/ahy3nz/mbuild …

4cebda0

…into hoomd_integration

ahy3nz added 2 commits November 4, 2019 18:24

Return ref values from hoomd

31b54ab

update docs

9e9a219

mattwthompson requested review from joaander and vyasr November 5, 2019 00:30

joaander approved these changes Nov 5, 2019

View reviewed changes

ahy3nz added 2 commits November 5, 2019 10:02

More version cehcking

d9402f0

Update neighborlist exclusions if rigid bodies are present in the hoo…

daaea95

…md snapshot

mikemhenry approved these changes Nov 6, 2019

View reviewed changes

Don't update neighborlist for rigid bodies

643a2a4

mattwthompson merged commit 7e15c3c into mosdef-hub:master Nov 6, 2019

umesh-timalsina referenced this pull request in GOMC-WSU/MoSDeF-GOMC Mar 22, 2022

Merge pull request #622 from ahy3nz/hoomd_integration

7f7189b

Add HOOMD snapshot and simulation converters

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Snapshot converter, hoomd simulation module #622

Snapshot converter, hoomd simulation module #622

ahy3nz commented Oct 16, 2019 •

edited

Loading

lgtm-com bot commented Oct 16, 2019

codecov bot commented Oct 16, 2019 •

edited

Loading

justinGilmer left a comment

justinGilmer Oct 17, 2019

justinGilmer Oct 17, 2019

justinGilmer Oct 17, 2019

ahy3nz commented Oct 17, 2019

lgtm-com bot commented Oct 17, 2019

ahy3nz commented Oct 17, 2019

lgtm-com bot commented Oct 17, 2019

ahy3nz commented Oct 18, 2019

lgtm-com bot commented Oct 18, 2019

uppittu11 commented Oct 21, 2019

ahy3nz commented Oct 21, 2019 •

edited

Loading

ahy3nz commented Oct 31, 2019

ahy3nz commented Oct 31, 2019

mikemhenry commented Nov 1, 2019

mikemhenry commented Nov 5, 2019 •

edited

Loading

mattwthompson commented Nov 5, 2019

joaander left a comment

joaander Nov 5, 2019

ahy3nz commented Nov 5, 2019

mikemhenry left a comment

ahy3nz commented Nov 6, 2019

mikemhenry commented Nov 6, 2019

	# Identify the uniqeu dihedral types before setting
	# Identify the unique dihedral types before setting

Snapshot converter, hoomd simulation module #622

Snapshot converter, hoomd simulation module #622

Conversation

ahy3nz commented Oct 16, 2019 • edited Loading

PR Summary:

PR Checklist

lgtm-com bot commented Oct 16, 2019

codecov bot commented Oct 16, 2019 • edited Loading

Codecov Report

justinGilmer left a comment

Choose a reason for hiding this comment

justinGilmer Oct 17, 2019

Choose a reason for hiding this comment

justinGilmer Oct 17, 2019

Choose a reason for hiding this comment

justinGilmer Oct 17, 2019

Choose a reason for hiding this comment

ahy3nz commented Oct 17, 2019

lgtm-com bot commented Oct 17, 2019

ahy3nz commented Oct 17, 2019

lgtm-com bot commented Oct 17, 2019

ahy3nz commented Oct 18, 2019

lgtm-com bot commented Oct 18, 2019

uppittu11 commented Oct 21, 2019

ahy3nz commented Oct 21, 2019 • edited Loading

ahy3nz commented Oct 31, 2019

ahy3nz commented Oct 31, 2019

mikemhenry commented Nov 1, 2019

mikemhenry commented Nov 5, 2019 • edited Loading

mattwthompson commented Nov 5, 2019

joaander left a comment

Choose a reason for hiding this comment

joaander Nov 5, 2019

Choose a reason for hiding this comment

ahy3nz commented Nov 5, 2019

mikemhenry left a comment

Choose a reason for hiding this comment

ahy3nz commented Nov 6, 2019

mikemhenry commented Nov 6, 2019

ahy3nz commented Oct 16, 2019 •

edited

Loading

codecov bot commented Oct 16, 2019 •

edited

Loading

ahy3nz commented Oct 21, 2019 •

edited

Loading

mikemhenry commented Nov 5, 2019 •

edited

Loading