diff --git a/DockerMakefiles/DockerMake.yml b/DockerMakefiles/DockerMake.yml index 7e3979a..62be7d7 100644 --- a/DockerMakefiles/DockerMake.yml +++ b/DockerMakefiles/DockerMake.yml @@ -10,7 +10,7 @@ _SOURCES_: - PySCF.yml - PythonTools.yml - SymMol.yml - + - LAMMPS.yml _ALL_: - ambertools @@ -24,7 +24,7 @@ _ALL_: - moldesign_complete - moldesign_notebook - moldesign_minimal - + - lammps ################################################## # Base image definitions diff --git a/DockerMakefiles/LAMMPS.yml b/DockerMakefiles/LAMMPS.yml new file mode 100644 index 0000000..7f4f6de --- /dev/null +++ b/DockerMakefiles/LAMMPS.yml @@ -0,0 +1,52 @@ +lammps_requirements: + description: Base libraries for BOTH running and building LAMMPS + build: | + RUN apt-get -y update + RUN apt-get -y install \ + bc \ + ssh \ + mpich2 \ + python-dev \ + libz-dev \ + fftw3-dev + + ENV PATH=${PATH}:/usr/lib + ENV PYTHONPATH=${PYTHONPATH}:/opt/lammps/python + ENV LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/opt/lammps/src + +lammps_build: + requires: + - buildbase + - lammps_requirements + description: Build image for LAMMPS. Based on https://github.com/malramsay64/lammps-docker/blob/master/Dockerfile + build: | + + RUN git clone https://github.com/lammps/lammps.git /opt/lammps && \ + cd /opt/lammps && \ + git fetch && \ + git checkout stable + + RUN cd /opt/lammps/src && \ + make yes-ASPHERE yes-BODY yes-CLASS2 yes-COLLOID yes-COMPRESS yes-CORESHELL yes-DIPOLE yes-GRANULAR yes-KSPACE yes-MANYBODY && \ + make yes-MC yes-MOLECULE yes-OPT yes-PERI yes-PERI yes-PYTHON yes-REPLICA yes-RIGID yes-SHOCK yes-SNAP yes-SRD + + # RUN cd /opt/lammps/src && \ + # make serial mode=shlib -j4 + + RUN cd /opt/lammps/src && \ + make mpi mode=shlib -j4 MPI_INC="-DMPICH_SKIP_MPICXX -I/usr/include/mpich" MPI_LIB="-Wl,-rpath,/usr/lib/mpich/lib -L/usr/lib/mpich/lib -lmpl -lmpich" + + RUN cd /opt/lammps/python && \ + python install.py + +lammps: + description: Deployable LAMMPS image + build_directory: buildfiles/lammps/ + requires: + - python_deploy_base + - lammps_requirements + copy_from: + lammps_build: + /usr/local/lib/python2.7/dist-packages: /usr/local/lib/python2.7/ + /usr/lib: /usr + /opt/lammps: /opt diff --git a/moldesign/_notebooks/Example 6. LAMMPS energy model and integrator.ipynb b/moldesign/_notebooks/Example 6. LAMMPS energy model and integrator.ipynb new file mode 100644 index 0000000..9901b39 --- /dev/null +++ b/moldesign/_notebooks/Example 6. LAMMPS energy model and integrator.ipynb @@ -0,0 +1,150 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import moldesign as mdt\n", + "from moldesign import units as u\n", + "import ipywidgets as ipy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from moldesign.models.lammps_model import LAMMPSPotential\n", + "import parmed as med\n", + "# mdt.configure()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mol = mdt.from_pdb('1yu8')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mol = mdt.Molecule(mol.residues[0:10])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mdmol = mdt.assign_forcefield(mol)\n", + "atom = mdmol.atoms[13]\n", + "print \"Done parameterizing\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print mdmol.write('pdb')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "######### Create LAMMPSInteractive model and run NVT simulation #########\n", + "mdmol.set_energy_model(mdt.models.LAMMPSInteractive)\n", + "\n", + "data = mdmol.energy_model.get_lammps_data()\n", + "print" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mdmol.calculate()\n", + "mdmol.set_integrator(integrator=mdt.integrators.LAMMPSLangevin, timestep=1.0*u.fs,\n", + " temperature=300.0*u.kelvin, frame_interval=50.0*u.fs, constrain_hbonds=False, constrain_water=False)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "mdmol.energy_model.apply_force(mdmol.atoms[11:20], [10.0, 0.0, 3.0] * u.kcalpermol/u.angstrom)\n", + "mytraj = mdmol.run(200 * u.fs)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/moldesign/_tests/test_mm.py b/moldesign/_tests/test_mm.py index 38bbd6d..bcfbca7 100644 --- a/moldesign/_tests/test_mm.py +++ b/moldesign/_tests/test_mm.py @@ -122,3 +122,54 @@ def test_minimization_reduces_energy(objkey, request): mol = request.getfuncargvalue(objkey) traj = mol.minimize() assert mol.calculate_potential_energy() < e1 + +@typedfixture('hasmodel') +def protein_stripped_amber_forcefield_lammps(): + mol = mdt.from_pdb('1yu8') + mol = mdt.Molecule([res for res in mol.residues if res.type == "protein"]) + newmol = mdt.assign_forcefield(mol) + + # Set Potential energy model + newmol.set_energy_model(mdt.models.LAMMPSPotential) + newmol.calculate() + + # Set Interactive energy model + newmol.set_energy_model(mdt.models.LAMMPSInteractive) + newmol.calculate() + newmol.energy_model.apply_force(atoms=newmol.atoms[11:20], vector=[10.0, 10.0, 0.0] * u.kcalpermol / u.angstrom) + newmol.energy_model.reset_force() + return newmol + +@typedfixture('hasmodel') +def protein_amber_forcefield_lammps(): + mol = mdt.from_pdb('1yu8') + newmol = mdt.assign_forcefield(mol) + + # Set Potential energy model + newmol.set_energy_model(mdt.models.LAMMPSPotential) + newmol.calculate() + + # Set Interactive energy model + newmol.set_energy_model(mdt.models.LAMMPSInteractive) + newmol.calculate() + newmol.energy_model.apply_force(atoms=newmol.atoms[11:20], vector=[10.0, 10.0, 0.0] * u.kcalpermol / u.angstrom) + newmol.energy_model.reset_force() + return newmol + +def test_lammps_shake_constrain_hbonds(protein_stripped_amber_forcefield_lammps): + mdmol = protein_stripped_amber_forcefield_lammps + mol = mdt.from_pdb('1yu8') + newmol = mdt.assign_forcefield(mol) + newmol.set_energy_model(mdt.models.LAMMPSPotential) + + # Test without constraints + mdmol.set_integrator(mdt.integrators.LAMMPSLangevin, timestep=2.0*u.fs, temperature=300*u.kelvin, frame_interval=10.0*u.fs, constrain_hbonds=False, constrain_water=False) + traj = mdmol.run(1000*u.fs) + + assert len(traj) > 0 + + # Test with constraints. Should raise NotImplementedError + mdmol.set_integrator(mdt.integrators.LAMMPSLangevin, timestep=2.0*u.fs, temperature=300*u.kelvin, frame_interval=10.0*u.fs, constrain_hbonds=True, constrain_water=True) + with pytest.raises(NotImplementedError): + traj = mdmol.run(1000*u.fs); + diff --git a/moldesign/integrators/__init__.py b/moldesign/integrators/__init__.py index 4637c0a..7678266 100644 --- a/moldesign/integrators/__init__.py +++ b/moldesign/integrators/__init__.py @@ -1,3 +1,4 @@ from . import base from .verlet import * from .openmm import * +from .lammps_integrator import * diff --git a/moldesign/integrators/lammps_integrator.py b/moldesign/integrators/lammps_integrator.py new file mode 100644 index 0000000..c8384a8 --- /dev/null +++ b/moldesign/integrators/lammps_integrator.py @@ -0,0 +1,201 @@ +# Copyright 2016 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from __future__ import absolute_import + +import moldesign as mdt +from moldesign import units as u +from moldesign.molecules import Trajectory +import numpy as np + +from .base import IntegratorBase, ConstantTemperatureBase + +try: + from lammps import lammps, PyLammps + # WARNING: this is the real library, not our interface - this works because of absolute + # imports. We should probably rename this interface +except ImportError: + print 'LAMMPS could not be imported; using remote docker container' + force_remote = True +else: # this should be configurable + force_remote = False # debugging + +import tempfile +import os +from copy import deepcopy +from itertools import islice + +def exports(o): + __all__.append(o.__name__) + return o +__all__ = [] + + +@exports +class LAMMPSLangevin(ConstantTemperatureBase): + + NAME_RESULT = "result" # Name used for dump + NAME_AFFECTED_ATOMS = "affected_atoms" # Name used for grouped atoms + NAME_LANGEVIN = "langevin_sim" # Name used for Langevin + NAME_ADDFORCE = "add_user_force" # Name used for addforce + + # TODO: raise exception if any constraints are requested ... + def __init__(self, *args, **kwargs): + super(LAMMPSLangevin, self).__init__(*args, **kwargs) + self._prepped = False # is the model prepped? + self._model = None + + self.lammps_system = None + self.total_iter = None # How many iterations of timesteps before end of simulation + self.output_iter = None # How often we dump the positions of the atoms during the simulation + self.traj = None + + def run(self, run_for): + """ + Users won't call this directly - instead, use mol.run + Propagate position, momentum by a single timestep using LAMMPS Langevin + + Arguments: + run_for (int): number of timesteps OR amount of time to run for + self.params.timestep (float): timestep length + + """ + self.total_iter = self.time_to_steps(run_for, self.params.timestep) + self.output_iter = self.time_to_steps(self.params.frame_interval, self.params.timestep) + + if self.total_iter < self.output_iter: + raise ValueError("Run duration {0} " + "can\'t be smaller than frame interval {1}".format(self.total_iter, self.output_iter)) + + # Prepare lammps system + self.prep() + self._configure_system() + lmps = self.lammps_system + + predictable_filename = "dump.result" + dump_filepath = os.path.join(self.mol.energy_model.tmpdir, predictable_filename) + self.mol.energy_model.files.append(dump_filepath) + + dump_command = "dump {0} all custom {1} {2} id x y z vx vy vz".format(self.NAME_RESULT, + self.output_iter, dump_filepath) + + lmps.command(dump_command) + lmps.command("dump_modify {0} sort id".format(self.NAME_RESULT)) + + # run simulation + lmps.run(self.total_iter, "post no") + + # Reset lammps system by undumping and unfixing + for dmp in lmps.dumps: + lmps.command("undump " + dmp.get('name')) + + for fix in lmps.fixes: + lmps.command("unfix " + fix.get('name')) + + # Create trajectory object + self.traj = Trajectory(self.mol) + self.mol.calculate() + + # Dynamics loop over the dump file + dump = open(dump_filepath, "rw+") + content = dump.readlines() + for istep in xrange(self.total_iter/self.output_iter): + self.step(istep, content) + self.traj.new_frame() + dump.close() + + # Set time + self.mol.time += self.time + + return self.traj + + def prep(self): + """ + Prepare the LAMMPS system by configuring it for Langevin simulation + + Arguments: + run_for (int): number of timesteps OR amount of time to run for + self.params.timestep (float): timestep length + + """ + # prepare lammps model prior to preparing the integrator + self._model = self.mol.energy_model + self._model.prep() + + self.time = 0.0 * self.params.timestep + + self._prepped = True + + def _configure_system(self): + # Get lammps object from model + lmps = self._model.lammps_system + + # Check if we need to constrain hbonds + if self.params.get('constrain_hbonds', False): + raise NotImplementedError('SHAKE not implemented') + + # Check if we need to constrain water + if self.params.get('constrain_water', False): + raise NotImplementedError('SHAKE not implemented') + + # Set timestep + lmps.command("timestep " + str(self.params.timestep.value_in(u.fs))) + + # Set Langevin settings + lmps.command("fix lange_nve all nve") + langevin_command = "fix {0} all langevin {1} {2} {3} 48279".format(self.NAME_LANGEVIN, + self.params.temperature.value_in(u.kelvin), + self.params.temperature.value_in(u.kelvin), 100.0) + lmps.command(langevin_command) + + # Apply force to group of atoms + if hasattr(self._model, 'affected_atoms') and hasattr(self._model, 'force_vector'): + if self._model.affected_atoms is not None and self._model.force_vector is not None: + # Group selected atoms + atoms = self._model.affected_atoms + + group_atom_cmd = "group {0} id ".format(self.NAME_AFFECTED_ATOMS) + for atom in atoms: + group_atom_cmd += "{0} ".format(atom.index+1) + lmps.command(group_atom_cmd.rstrip()) + + # Apply force to selected atoms + force_in_real_units = self._model.force_vector + lmps.command("fix {0} all addforce {1} {2} {3}".format(self.NAME_ADDFORCE, force_in_real_units[0], + force_in_real_units[1], force_in_real_units[2])) + + # Set thermo style + lmps.command("thermo_style custom step") + self.lammps_system = lmps # Save lammps configuration + + def step(self, idx, results): + + """ + Update molecule's position and velocity at each step of the simulation + + """ + lmps = self.lammps_system + # start at 9, stop at 9+L.atoms.natoms + # start at 14+L.atoms.natoms + + tmp = list(islice(results, 9*(idx+1) + lmps.atoms.natoms*idx, (lmps.atoms.natoms+9)*(idx+1))) + pos = list(item.split()[1:4] for item in tmp) + vel = list(item.split()[4:] for item in tmp) + + pos_array = np.array(pos).astype(np.float) + vel_array = np.array(vel).astype(np.float) + + self.mol.positions = pos_array * u.angstrom + self.mol.velocities = vel_array * u.angstrom / u.fs + + self.time += self.params.frame_interval \ No newline at end of file diff --git a/moldesign/models/__init__.py b/moldesign/models/__init__.py index 3288a85..4a0cf52 100644 --- a/moldesign/models/__init__.py +++ b/moldesign/models/__init__.py @@ -5,3 +5,4 @@ from .amber import * from .openbabel import * from .nwchem import * +from .lammps_model import * diff --git a/moldesign/models/lammps_model.py b/moldesign/models/lammps_model.py new file mode 100644 index 0000000..45d9b69 --- /dev/null +++ b/moldesign/models/lammps_model.py @@ -0,0 +1,462 @@ +# Copyright 2016 Autodesk Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from __future__ import absolute_import + +import moldesign as mdt +import moldesign.units as u + +from moldesign import compute +from moldesign import forcefields as ff + +from .base import EnergyModelBase + +try: + from lammps import lammps, PyLammps + # WARNING: this is the real library, not our interface - this works because of absolute + # imports. We should probably rename this interface +except ImportError: + print 'LAMMPS could not be imported; using remote docker container' + force_remote = True +else: # this should be configurable + force_remote = False # debugging + +debugging = True # debugging + +import parmed as med +import tempfile +import os +import sys + +def exports(o): + __all__.append(o.__name__) + return o +__all__ = [] + + +@exports +class LAMMPSPotential(EnergyModelBase): + """Creates a LAMMPS system to drive energy calculations. + """ + # NEWFEATURE: need to set/get platform (and properties, e.g. number of threads) + DEFAULT_PROPERTIES = ['potential_energy', 'forces'] + + def __init__(self, **kwargs): + super(LAMMPSPotential, self).__init__(**kwargs) + self._prepped = False # is the model prepped? + self._lammps_data_path = None + + # For temporary directory + # create temporary file system + # Create temporary filesystem + self.tmpdir = tempfile.mkdtemp() + self.saved_umask = os.umask(0077) + self.files = [] + + self.lammps_system = None # Basic LAMMPS system for this molecule + self.unit_system = None + + # TODO: + # def __del__(self, exc_type, exc_value, traceback): + # print "Deleting LAMMPSPotential" + + # # Unlink all files + # for f in self.files: + # os.unlink(f) + # os.remove(f) + + # # securely remove temporary file system + # os.umask(self.saved_umask) + # os.rmdir(self.tmpdir) + + # # Delete lammps object + # del self.lammps_system + + def prep(self): + """ + Drive the construction of the LAMMPS simulation + This will rebuild this OpenMM simulation if: A) it's not built yet, or B) + there's a new integrator + """ + + if not self._prepped or self.lammps_system is None: + # Create lammps + self._create_lammps() + self.unit_system = u.UnitSystem(length=u.angstrom, force=u.kcalpermol / u.fs, energy=u.kcalpermol, + time=u.fs, + mass=u.gpermol, charge=u.electron_charge) + + self._prepped = True + + # calculate potential energy and force + def calculate(self, requests=None): + + # Recreate system if molecule changed + results = self._evaluate_molecule() + return results + + ################################################# + # "Private" methods for managing LAMMPS are below + + def _evaluate_molecule(self): + # Recreate system if molecule changed + self.prep() + + # Run for 0 fs to evaluate system + my_lmps = self.lammps_system + my_lmps.command("thermo_style one") + my_lmps.run(0) + + # Evaluate potential energy and forces + pe = my_lmps.eval("pe") + force_array = [] + for i in range(0, my_lmps.atoms.natoms): + force_array.append(my_lmps.atoms[i].force) + return {'potential_energy': pe * u.kcalpermol, + 'forces': force_array * u.kcalpermol / u.angstrom} + + # NOTE: Used by workflows to create lammps formatted data file + def get_lammps_data_path(self): + self.prep() + + # initialize lammps object + return self._lammps_data_path + + def _create_lammps(self): + + """ + Create a LAMMPS system. Use MDT molecule object to construct a LAMMPS system + + """ + + # Ensure force fields are assigned + if self.mol.ff is None or self.mol.ff.parmed_obj is None: + raise NotImplementedError('Assign forcefield parameters to the system') + + # initialize lammps object + lmp = lammps() + pylmp = PyLammps(ptr=lmp) + + # Create temperary data file + predictable_filename = 'data.lammps_mol' + data_path = os.path.join(self.tmpdir, predictable_filename) + + # Add file paths to list of files + self.files.append(data_path) + self._lammps_data_path = data_path + + with open(data_path, "w") as lammps_data: + formatted_data = self._get_lammps_data() + lammps_data.write(formatted_data) + + # Set basic parameters. + pylmp.command("units real") + pylmp.command("dimension 3") + pylmp.command("atom_style full") + pylmp.command("pair_style lj/charmm/coul/charmm/implicit 8.0 10.0") + + pylmp.command("bond_style harmonic") + pylmp.command("angle_style harmonic") + pylmp.command("dihedral_style harmonic") + pylmp.command("improper_style harmonic") + + pylmp.command("dielectric 4.0") + pylmp.command("read_data " + data_path) + + # Assign lammps system + self.lammps_system = pylmp + + def _get_lammps_data(self): + """ + Parse molecule and force fields info to generate a formatted data + + """ + + parmedmol = self.mol.ff.parmed_obj + + # Get non-bond atom types + parmedmol.fill_LJ() + max_idx = 0 + for nonbond_idx in parmedmol.LJ_types.itervalues(): + if nonbond_idx > max_idx: + max_idx = nonbond_idx + + datalines = "" + + datalines += "LAMMPS Description\r\n\r\n" + datalines += "{0} atoms\r\n".format(len(parmedmol.atoms)) + datalines += "{0} bonds\r\n".format(len(parmedmol.bonds)) + datalines += "{0} angles\r\n".format(len(parmedmol.angles)) + datalines += "{0} dihedrals\r\n".format(len(parmedmol.dihedrals)) + datalines += "{0} impropers\r\n".format(len(parmedmol.impropers)) + datalines += "\r\n" + + datalines += "{0} atom types\r\n".format(max_idx) + if len(parmedmol.bond_types) > 0: + datalines += "{0} bond types\r\n".format(len(parmedmol.bond_types)) + if len(parmedmol.angle_types) > 0: + datalines += "{0} angle types\r\n".format(len(parmedmol.angle_types)) + if len(parmedmol.dihedral_types) > 0: + datalines += "{0} dihedral types\r\n".format(len(parmedmol.dihedral_types)) + if len(parmedmol.improper_types) > 0: + datalines += "{0} improper types\r\n".format(len(parmedmol.improper_types)) + + datalines += "\r\n" + + # TODO: calculate lo and hi coordinates for Box size + xlo = xhi = ylo = yhi = zlo = zhi = None + for atom in self.mol.atoms: + if xlo == None: + xlo = atom.x.value_in(u.angstrom) + xhi = atom.x.value_in(u.angstrom) + ylo = atom.y.value_in(u.angstrom) + yhi = atom.y.value_in(u.angstrom) + zlo = atom.z.value_in(u.angstrom) + zhi = atom.z.value_in(u.angstrom) + else: + xlo = min(xlo, atom.x.value_in(u.angstrom)) + xhi = max(xhi, atom.x.value_in(u.angstrom)) + ylo = min(ylo, atom.y.value_in(u.angstrom)) + yhi = max(yhi, atom.y.value_in(u.angstrom)) + zlo = min(zlo, atom.z.value_in(u.angstrom)) + zhi = max(zhi, atom.z.value_in(u.angstrom)) + + datalines += "{0} {1} xlo xhi\r\n".format(xlo-50, xhi+50) + datalines += "{0} {1} ylo yhi\r\n".format(ylo-50, yhi+50) + datalines += "{0} {1} zlo zhi\r\n".format(zlo-50, zhi+50) + + datalines += "\r\n" + + # Masses - atom types + datalines += "Masses\r\n\r\n" + for i in range(1, max_idx+1): + for atom in parmedmol.atoms: + if atom.nb_idx == i: + datalines += "{0} {1}\r\n".format(i, atom.mass) + break + datalines += "\r\n" + + # Pair coefficients + datalines += "Pair Coeffs\r\n\r\n" + for i in range(1, max_idx+1): + for atom in parmedmol.atoms: + if atom.nb_idx == i: + datalines += "{0} {1} {2} {3} {4}\r\n".format(i, atom.epsilon, atom.sigma, atom.epsilon_14, atom.sigma_14) + break + datalines += "\r\n" + + # Bond Coeffs - Harmonic + + if parmedmol.bond_types: + datalines += "Bond Coeffs\r\n\r\n" + for bt in parmedmol.bond_types: + datalines += "{0} {1} {2}\r\n".format(bt.idx+1, bt.k, bt.req) + + datalines += "\r\n" + + # Angle Coeffs - Harmonic + if parmedmol.angle_types: + datalines += "Angle Coeffs\r\n\r\n" + for at in parmedmol.angle_types: + datalines +="{0} {1} {2}\r\n".format(at.idx+1, at.k, at.theteq) + + datalines +="\r\n" + + + # Dihedral Coeffs - Charmm + if parmedmol.dihedral_types: + datalines += "Dihedral Coeffs\r\n\r\n" + for dt in parmedmol.dihedral_types: + datalines += "{0} {1} {2} {3}\r\n".format(dt.idx+1, dt.phi_k, (-1 if dt.phase >= 180 else 1), dt.per) + + datalines += "\r\n" + + + # Improper Coeffs - Harmonic + if parmedmol.improper_types: + datalines += "Improper Coeffs\r\n\r\n" + for it in parmedmol.improper_types: + datalines += "{0} {1} {2}\r\n".format(it.idx+1, it.psi_k, it.psi_eq) + + datalines += "\r\n" + + + # print atoms + datalines += "Atoms\r\n\r\n" + for atom in parmedmol.atoms: + datalines += "{0} {1} {2} {3} ".format(atom.idx+1, atom.residue.idx+1, atom.nb_idx, atom.charge) + if not self.mol.positions.any(): + datalines += "{0} {1} {2} ".format(0.00, 0.00, 0.00) + else: + pos = self.mol.positions[atom.idx] + datalines += "{0} {1} {2} ".format(pos[0].value_in(u.angstrom), pos[1].value_in(u.angstrom), + pos[2].value_in(u.angstrom)) + datalines += "{0} {1} {2}\r\n".format(0, 0, 0) # nx, ny, nz not defined + + datalines += "\r\n" + + + # print velocities + datalines += "Velocities\r\n\r\n" + for atom in self.mol.atoms: + xvel = atom.vx.value_in(u.angstrom/u.fs) + yvel = atom.vy.value_in(u.angstrom/u.fs) + zvel = atom.vz.value_in(u.angstrom/u.fs) + datalines += "{0} {1} {2} {3}\r\n".format(atom.index+1, xvel, yvel, zvel) + + datalines += "\r\n" + + + # print bonds + datalines += "Bonds\r\n\r\n" + for idx, bond in enumerate(parmedmol.bonds): + datalines += "{0} {1} {2} {3}\r\n".format(idx+1, bond.type.idx+1, + bond.atom1.idx+1, bond.atom2.idx+1) + + datalines += "\r\n" + + + # print angles + datalines += "Angles\r\n\r\n" + for idx, angle in enumerate(parmedmol.angles): + datalines += "{0} {1} {2} {3} {4}\r\n".format(idx+1, angle.type.idx+1, + angle.atom1.idx+1, + angle.atom2.idx+1, angle.atom3.idx+1) + + datalines += "\r\n" + + + # print dihedrals + if parmedmol.dihedrals: + datalines += "Dihedrals\r\n\r\n" + for idx, di in enumerate(parmedmol.dihedrals): + datalines += "{0} {1} {2} {3} {4} {5}\r\n".format(idx+1, di.type.idx+1, + di.atom1.idx+1, di.atom2.idx+1, + di.atom3.idx+1, di.atom4.idx+1) + + datalines += "\r\n" + + + # print impropers + if parmedmol.impropers: + datalines += "Impropers\r\n\r\n" + for im in parmedmol.impropers: + datalines += "{0} {1} {2} {3} {4} {5}\r\n".format(im.idx+1, im.type.idx+1, + im.atom1.idx+1, im.atom2.idx+1, + im.atom3.idx+1, im.atom4.idx+1) + + datalines += "\r\n" + + # Need both atom types and LAMMPS formatted data for web + return datalines + +@exports +class LAMMPSInteractive(EnergyModelBase): + + def __init__(self, **kwargs): + super(LAMMPSInteractive, self).__init__(**kwargs) + self._prepped = False # is the model prepped? + + # LAMMPSPotential model object + self._potential_model = None + self.lammps_system = None # Get from potential model + self.unit_system = None # Get from potential model + self.tmpdir = None + self.saved_umask = None + self.files = None + + self._calculated_properties = None + self.affected_atoms = None # atom that the user is directly interacting with + self.force_vector = None # force that acts on the molecule + + # TODO: + # def __del__(self, exc_type, exc_value, traceback): + # print "Deleting LAMMPSInteractive" + + # # Delete LAMMPS potential model + # del self._potential_model + + def prep(self): + """ + Drive the construction of the LAMMPS simulation + This will rebuild this OpenMM simulation if: A) it's not built yet, or B) + there's a new integrator + """ + + if not self._prepped or self._potential_model is None or self.lammps_system is None: + self._potential_model = mdt.models.LAMMPSPotential() + self._potential_model.mol = self.mol + + self._potential_model.prep() + + self._calculated_properties = None + + # Get settings from LAMMPSPotential model + self.lammps_system = self._potential_model.lammps_system + self.unit_system = self._potential_model.unit_system + self.tmpdir = self._potential_model.tmpdir + self.saved_umask = self._potential_model.saved_umask + self.files = self._potential_model.files + + self._prepped = True + + def calculate(self, requests): + ''' + After assigning this energy model to the molecule, call this function as soon as possible, so + the molecule is ready for the interactive simulation + ''' + + self.prep() + + # Recreate system if molecule changed + if self._calculated_properties is None: + self._calculated_properties = self._potential_model.calculate() + + return self._calculated_properties + + def apply_force(self, atoms, vector): + ''' + :param atoms: atoms to apply force + :param vector: amount of force to apply to atoms in kcal/mol/angstrom + ''' + if atoms is None or vector is None: + raise ValueError('One or more argument(s) missing!') + + vector_real_units = vector.value_in(u.kcalpermol / u.angstrom) + + # Set internal variables + self.affected_atoms = atoms + self.force_vector = vector_real_units + + def reset_force(self): + ''' + Reset user-added force + ''' + self.affected_atoms = None + self.force_vector = None + + def get_lammps_data(self): + self.prep() + lammps_data_path = self._potential_model.get_lammps_data_path() + + export_data_string = '' + with open(lammps_data_path) as data: + lines = data.readlines() + for line in lines: + export_data_string += line + data.close() + + return export_data_string + + + diff --git a/moldesign/molecules/molecule.py b/moldesign/molecules/molecule.py index 95956d7..bb3fb14 100644 --- a/moldesign/molecules/molecule.py +++ b/moldesign/molecules/molecule.py @@ -699,6 +699,23 @@ def run(self, run_for): print 'Done - integrated "%s" from %s to %s' % (self, init_time, self.time) return traj + def run_simulation(self, force_vector): + """ Starts the integrator's default integration + + Args: + run_for (int or [time]): number of steps or amount of time to run for + + Returns: + moldesign.trajectory.Trajectory + """ + if type(self.integrator) is not mdt.integrators.LAMMPSNvt: + raise ValueError('Cannot simulate; no integrator set for %s' % self) + + init_time = self.time + traj = self.integrator.run_simulation(force_vector) + print 'Done - integrated "%s" from %s to %s' % (self, init_time, self.time) + return traj + def calculate(self, requests=None, wait=True, use_cache=True): """ Runs a potential energy calculation on the current geometry, returning the requested quantities.