Facade

foldFacade.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+from simtk.openmm.app import *
+from simtk.openmm import *
+from simtk.unit import *
+from sys import stdout
+import os
+import sys
+import argparse
+
+from parser import Parser
+from unfold import Unfold
+from simulate import Simulate
+from modeller import ModellerF
+
+args = Parser(argparse).parse()
+
+try:
+    platform = Platform.getPlatformByName("CUDA")
+except Exception:
+    platform = Platform.getPlatformByName("OpenCL")
+
+pdb = Unfold(args).unfold
+forcefield = ForceField('amber03.xml', 'amber03_obc.xml')
+
+ModellerF(forcefield).model()
+
+Simulate(forcefield, modeller, args).simulate()

modeller.py

@@ -0,0 +1,12 @@
+import sys
+
+
+class ModellerF:
+    def __init__(self, forcefield):
+        self.forcefield = forcefield
+
+    def model(self):
+       #forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
+        modeller = Modeller(pdb.topology, pdb.positions)
+        modeller.addHydrogens(self.forcefield)
+        print(modeller.topology)

simulate.py

@@ -0,0 +1,29 @@
+import sys
+
+
+class Simulate:
+    def __init__(self, forcefield, modeller, args):
+        self.forcefield = forcefield
+        self.modeller = modeller
+        self.args = args
+
+    def simulate(self):
+        system = self.forcefield.createSystem(self.modeller.topology,
+                                              # matches https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2980750/#bib39
+                                              implicitSolvent=OBC2,
+                                              nonbondedMethod=NoCutoff, nonbondedCutoff=1*nanometer,
+                                              constraints=HBonds)
+        integrator = LangevinIntegrator(
+            self.args.temp*kelvin, 1/picosecond, 2*femtoseconds)
+        simulation = Simulation(
+            self.modeller.topology, system, integrator, platform)
+        simulation.context.setPositions(self.modeller.positions)
+        simulation.minimizeEnergy()
+
+        steps = self.args.steps
+        steps_write = max(1, steps//self.args.writes)
+        print("writing every %d steps" % steps_write)
+        simulation.reporters.append(PDBReporter(args.out, steps_write))
+        simulation.reporters.append(StateDataReporter(
+            stdout, steps_write, step=True, potentialEnergy=True, temperature=True))
+        simulation.step(steps)

unfold.py

@@ -0,0 +1,23 @@
+import sys
+
+
+class Unfold:
+    def __init__(self, args):
+        self.args = args
+
+    def unfold(self):
+        if self.args.scratch:
+            # unfolded protein
+            if self.args.fasta is not None:
+                fasta = self.args.fasta
+            else:
+                protein_fasta = "proteins/villin/1vii.fasta"
+                fasta = open(protein_fasta).read().split("\n")[1]
+            print("folding %s" % fasta)
+            from lib import write_unfolded
+            write_unfolded(fasta, "/tmp/unfolded.pdb")
+            pdb = PDBFile("/tmp/unfolded.pdb")
+        else:
+            # already folded protein
+            pdb = PDBFile(self.args.pdb)
+        return pdb