-
Notifications
You must be signed in to change notification settings - Fork 37
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
base: master
Are you sure you want to change the base?
Conversation
There was a problem hiding this 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.
moldesign/integrators/lammps.py
Outdated
@exports | ||
class LAMMPSNvt(ConstantTemperatureBase): | ||
def __init__(self, *args, **kwargs): | ||
super(LAMMPSNvt, self).__init__(*args, **kwargs) |
There was a problem hiding this comment.
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.
moldesign/integrators/lammps.py
Outdated
""" | ||
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 |
There was a problem hiding this comment.
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)
moldesign/integrators/lammps.py
Outdated
# if istep + 1 >= next_trajectory_frame: | ||
# self.traj.new_frame() | ||
# next_trajectory_frame += self.params.frame_interval | ||
# return self.traj |
There was a problem hiding this comment.
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:
- Send the molecule's current positions and momenta to LAMMPS
- Create a
Trajectory
object to store the results. - Run dynamics in LAMMPS, storing new a frame to the
Trajectory
object at fixed simulation time intervals (defined inself.params.frame_interval
)
moldesign/integrators/lammps.py
Outdated
lammps_system.command(nvt_command) | ||
|
||
# TODO: | ||
if self.params.constrain_hbonds && self.model.group_hbond == True : |
There was a problem hiding this comment.
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
moldesign/models/lammps.py
Outdated
def calculate(self, requests): | ||
|
||
lammps_system = self.lammps_system | ||
lammps_system.run(500) |
There was a problem hiding this comment.
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
moldesign/models/lammps.py
Outdated
# 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 |
There was a problem hiding this comment.
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
moldesign/models/lammps.py
Outdated
parmedmol = med.load_file('/tmp/tmp.prmtop') | ||
|
||
lmp = lammps() | ||
L = PyLammps(ptr=lmp) |
There was a problem hiding this comment.
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
moldesign/models/lammps.py
Outdated
self.group_water = False | ||
|
||
self.lammps_system = L | ||
print 'Created LAMMPS object' |
There was a problem hiding this comment.
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:
- We made potentially unexpected changes to the
Molecule
object. - The function has a somewhat ambiguous effect that needs to be logged
- 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.
moldesign/models/lammps.py
Outdated
print 'Created LAMMPS object' | ||
|
||
|
||
def _group_hbonds(parmedmol): |
There was a problem hiding this comment.
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
moldesign/models/lammps.py
Outdated
# See http://parmed.github.io/ParmEd/html/index.html for ParmEd units | ||
|
||
def _create_lammps_data(parmedmol): | ||
|
There was a problem hiding this comment.
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.
Fix visualization problems
…ems like it's broken
In an attemp to fix Trajectory, merging master
…nd due to small box size
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
There was a problem hiding this 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)
moldesign/models/amber2lammps.py
Outdated
@@ -0,0 +1,945 @@ | |||
#! /usr/bin/evn python2 |
There was a problem hiding this comment.
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
moldesign/models/lammps_model.py
Outdated
import tempfile | ||
import os | ||
import sys | ||
from .amber2lammps import * |
There was a problem hiding this comment.
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.
moldesign/models/lammps_model.py
Outdated
self.lammps_system = None # Basic LAMMPS system for this molecule | ||
self.unit_system = None | ||
|
||
def __exit__(self, exc_type, exc_value, traceback): |
There was a problem hiding this comment.
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.
There was a problem hiding this 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, |
There was a problem hiding this comment.
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) |
There was a problem hiding this comment.
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 |
There was a problem hiding this comment.
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:
- if
self.params.get('constrain_hbonds', False) == True
, then you'll need to set up SHAKE to constrain all bonds involving a hydrogen atom - 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 nameHOH
orWAT
) - You'll need to examine the list of constraints in
self.mol.constraints
. This is a list of objects that inherit fromGeometryConstraint
(seemoldesign/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 :)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed with ae97b6b
moldesign-examples/Untitled1.ipynb
Outdated
@@ -0,0 +1,1913 @@ | |||
{ |
There was a problem hiding this comment.
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
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
moved in ae97b6b
@dane1122 - a couple issues I'm running into:
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 It looks like it's looking for an attribute that's in the Interactive potential class.
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 |
moldesign/models/lammps_model.py
Outdated
|
||
@exports | ||
class LAMMPSInteractive(EnergyModelBase): | ||
|
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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?
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? |
…esign/_notebooks directory
Update test: 0e4f1f4 |
Now directly access parmed object.
@dane1122 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 |
Pull request!