Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Interactive simulation #129

Open
wants to merge 50 commits into
base: master
Choose a base branch
from
Open

Interactive simulation #129

wants to merge 50 commits into from

Conversation

dane1122
Copy link

@dane1122 dane1122 commented Feb 1, 2017

Pull request!

Copy link
Contributor

@avirshup avirshup left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like a great start! Just a few more things need to get hooked up before we're ready to start running it.

@exports
class LAMMPSNvt(ConstantTemperatureBase):
def __init__(self, *args, **kwargs):
super(LAMMPSNvt, self).__init__(*args, **kwargs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can omit this __init__ method - the superclass __init__ will get called automatically if you do.

"""
Users won't call this directly - instead, use mol.run
Propagate position, momentum by a single timestep using velocity verlet
:param run_for: number of timesteps OR amount of time to run for
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like this docstring came from a piece of code I haven't touched in a while :)

We've switched over to Google's docstring format - see http://google.github.io/styleguide/pyguide.html?showone=Comments#Comments

So you'll want to make sure to write your docstrings in that format (or update mine, as the case may be)

# if istep + 1 >= next_trajectory_frame:
# self.traj.new_frame()
# next_trajectory_frame += self.params.frame_interval
# return self.traj
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Once we're ready to get started with the run method, it will need to:

  1. Send the molecule's current positions and momenta to LAMMPS
  2. Create a Trajectory object to store the results.
  3. Run dynamics in LAMMPS, storing new a frame to the Trajectory object at fixed simulation time intervals (defined in self.params.frame_interval )

lammps_system.command(nvt_command)

# TODO:
if self.params.constrain_hbonds && self.model.group_hbond == True :
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No && or || in python - it's just and and or

def calculate(self, requests):

lammps_system = self.lammps_system
lammps_system.run(500)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This part should be in the integrator - we won't do any running in the EnergyModel

# Store the PRMTOP file if ff are assigned and use it to create lammps molecule data
self.mol.ff.amber_params.prmtop.put('/tmp/tmp.prmtop')

import parmed as med
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should go at the top of the function declaration

parmedmol = med.load_file('/tmp/tmp.prmtop')

lmp = lammps()
L = PyLammps(ptr=lmp)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable name should be lower case - something like pylmp would work here

self.group_water = False

self.lammps_system = L
print 'Created LAMMPS object'
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We generally want to keep console logging like this pretty light, except in the following situations:

  1. We made potentially unexpected changes to the Molecule object.
  2. The function has a somewhat ambiguous effect that needs to be logged
  3. We're in the middle of a long simulation and need to provide status updates

For OpenMMPotential, we log the system creation because it falls into case 2). OpenMM systems bind to a specific hardware interface - the CPU, an nVidia GPU, or a generic GPU - so we need to tell the user which resource the simulation will be using. But here, there's nothing ambiguous, so we don't need to log it.

print 'Created LAMMPS object'


def _group_hbonds(parmedmol):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Everyone is bad at this, me especially, but it's still good to write at least a one-line docstring, even for private functions, unless it's really really really obvious

# See http://parmed.github.io/ParmEd/html/index.html for ParmEd units

def _create_lammps_data(parmedmol):

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You'll of course need actually open a lammps_data file-like object here. It would be great if we could pass this file directly through the LAMMPS API so that we don't need to write anything to disk, but that's not always possible unfortunately.

@dane1122
Copy link
Author

dane1122 commented Apr 7, 2017

Hey Aaron, could we merge this into master so we can use it for interactive simulation MST app workflow? It would be great if you can review the code for this!

…ds to ensure that molecule can move around
Copy link
Contributor

@avirshup avirshup left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Day - the LAMMPSPotential class in general looks good! A few Python style comments here (will tackle the integrator and the interactive model next)

@@ -0,0 +1,945 @@
#! /usr/bin/evn python2
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's move this file to moldesign/external and import it from there, just for consistency

import tempfile
import os
import sys
from .amber2lammps import *
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

import * is discouraged in most contexts because it makes it hard to track down variables. This should be refactored to from ..external import amber2lammps, and all variables from that module should be referred to as amber2lammps.Lammps etc.

self.lammps_system = None # Basic LAMMPS system for this molecule
self.unit_system = None

def __exit__(self, exc_type, exc_value, traceback):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you checked to see if this method gets executed? I think you'll want to rename it to __del__ if you want it to run when this object is garbage-collected.

Copy link
Contributor

@avirshup avirshup left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few things to look at for the LAMMPS integrator

lmps.command("timestep " + str(self.params.timestep.value_in(u.fs)))

# Set NVT settings
nvt_command = "fix {0} all nvt temp {1} {2} {3}".format(self.NAME_NVT_SIM,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's make sure to apply @ktbolt's fix for the Langevin integrator here (i.e., use nve but apply a Langevin thermostat)

lmps.command("unfix " + fix.get('name'))

# Create trajectory object
self.traj = Trajectory(self.mol, unit_system=self._model.unit_system)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't specify unit_system here - it's actually for output, not input, so MDT will set it automatically.

self._prepped = True

def _configure_system(self):
# Get lammps object from model
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We'll need to make sure to tell LAMMPS to handle geometry constraints. Here are the 3 things that the interface will need to handle:

  1. if self.params.get('constrain_hbonds', False) == True, then you'll need to set up SHAKE to constrain all bonds involving a hydrogen atom
  2. if self.params.get('constrain_water', False) == True, then you'll need to set up SHAKE to constrain all bond lengths in any water molecules (residue name HOH or WAT)
  3. You'll need to examine the list of constraints in self.mol.constraints. This is a list of objects that inherit from GeometryConstraint (see moldesign/geom/constraints.py), each of which describes a constraint that needs to be satisfied during MD.

Alternatively, you can just raise NotImplementedError in any/all of these cases :)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed with ae97b6b

@@ -0,0 +1,1913 @@
{
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please go ahead and move this notebook to moldesign/_notebooks, and then you can delete the rest of the moldesign-examples directory

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

moved in ae97b6b

@avirshup
Copy link
Contributor

avirshup commented Apr 11, 2017

@dane1122 - a couple issues I'm running into:

  • The LAMMPS integrator throws an error when using the base LAMMPS potential model, e.g.:
import moldesign as mdt
from moldesign import units as u
mol = mdt.from_pdb('1yu8')
protein = mdt.assign_forcefield(mol)
protein.set_energy_model(mdt.models.lammps_model.LAMMPSPotential)
protein.set_integrator(mdt.integrators.LAMMPSNvt)
protein.run(5.0*u.ps)

Raises AttributeError: 'LAMMPSPotential' object has no attribute 'affected_atoms'

It looks like it's looking for an attribute that's in the Interactive potential class.

  • We'll need to disable LAMMPS' built-in logging when debugging == False so that the notebook doesn't get swamped with messages like the following:
Default LAMMPS logging (click to embiggen)

LAMMPS output is captured by PyLammps wrapper

Welcome to amber2lammps, a program to convert Amber files to Lammps format!

Reading /tmp/tmp.crd... default_name
12
0
done.
Reading /tmp/tmp.top... default_name
12
Reading Charges...
Reading Atomic Number...
Reading Atomic Masses...
Reading Atom Types...
Reading Excluded Atoms...
Reading Non-bonded Parameter Index...
Reading Residue Labels...
Reading Residues Starting Pointers...
Reading Bond Force Constants...
Reading Equilibrium Bond Values...
Reading Angle Force Constants...
Reading Equilibrium Angle Values...
Reading Dihedral Force Constants...
Reading Dihedral Periodicity...
Reading Dihedral Phase...
Reading 1-4 Electrostatic Scaling Factor...
Reading 1-4 Van der Waals Scaling Factor...
Reading Solty...
Reading LJ A Coefficient...
Reading LJ B Coefficient...
Reading Bonds which include hydrogen...
Reading Bonds which dont include hydrogen...
Reading Angles which include hydrogen...
Reading Angles which dont include hydrogen...
Reading Dihedrals which include hydrogen...
Reading Dihedrals which dont include hydrogen...
Reading Excluded Atom List...
Reading H-Bond A Coefficient, corresponding to r12 term for all possible types...
Reading H-Bond B Coefficient, corresponding to r
10 term for all possible types...
Reading H-Bond Cut...
Reading Amber Atom Types for each atom...
Reading Tree Chain Classification...
Reading Join Array: Tree joining information
Reading IRotate...
(warning: File too large!) done.
Converting... (warning: Guessing at periodic box!) done.
Writing /var/folders/8d/k0qh4ksd5nv5b2vp6plfx6240000gp/T/tmpDdrFfg/data.lammps_mol... done.


@exports
class LAMMPSInteractive(EnergyModelBase):

Copy link
Contributor

@avirshup avirshup Apr 11, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the purposes of getting the PR merged into MDT, we don't need to worry about interactive potentials (we'll deal with that in the web app). So let's leave this class out of the current PR, and merge it in later.

Copy link

@ktbolt ktbolt Apr 11, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should the "pair_style lj/charmm/coul/long 8.0 10.0 10.0" be changed to "pair_style lj/charmm/coul/charmm/implicit 8.0 10.0" or is this for something else?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. How do I leave it out in the current PR, do you know?

Dayeong Lee added 2 commits April 10, 2017 17:26
…zhi coords to ensure that molecule can move around"

This reverts commit bc562c3.

Revert "Lammps model and integrator"

This reverts commit 4454121.

Format using parmed instead of using amber2lammps.py since it is being maintained
@dane1122
Copy link
Author

I think I fixed the bugs that you listed. Could you please verify for me? All the fixes should be in either one of the commits:
ae97b6b
418f022

@avirshup
Copy link
Contributor

Code is looking good! I'll play with it a bit when I get into the office. In the meantime, could you check on these two test failures?
test_lammps_shake_constrain_hbonds
test_minimization_reduces_energy

@dane1122
Copy link
Author

Update test: 0e4f1f4

@avirshup
Copy link
Contributor

avirshup commented Apr 18, 2017

@dane1122 looks like we're getting a NameError at lammps_model.py:147

Oh right, this because LAMMPS isn't installed on Travis, and we're not building the docker container yet (future PR). So never mind on this one

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants