|
| 1 | +#!/usr/bin/env python3 |
| 2 | + |
| 3 | +import sys |
| 4 | +import argparse |
| 5 | +import importlib |
| 6 | +import numpy as np |
| 7 | + |
| 8 | +from pysages.colvars import DihedralAngle |
| 9 | +from pysages.methods import UmbrellaIntegration, SerialExecutor |
| 10 | +from pysages.utils import try_import |
| 11 | + |
| 12 | +import pysages |
| 13 | + |
| 14 | +openmm = try_import("openmm", "simtk.openmm") |
| 15 | +unit = try_import("openmm.unit", "simtk.unit") |
| 16 | +app = try_import("openmm.app", "simtk.openmm.app") |
| 17 | + |
| 18 | + |
| 19 | +def generate_simulation(**kwargs): |
| 20 | + """ |
| 21 | + Generates a simulation context, we pass this function to `pysages.run`. |
| 22 | + """ |
| 23 | + pdb_filename = "adp-explicit.pdb" |
| 24 | + T = 298.15 * unit.kelvin |
| 25 | + dt = 2.0 * unit.femtoseconds |
| 26 | + pdb = app.PDBFile(pdb_filename) |
| 27 | + |
| 28 | + ff = app.ForceField("amber99sb.xml", "tip3p.xml") |
| 29 | + cutoff_distance = 1.0 * unit.nanometer |
| 30 | + topology = pdb.topology |
| 31 | + system = ff.createSystem( |
| 32 | + topology, |
| 33 | + constraints=app.HBonds, |
| 34 | + nonbondedMethod=app.NoCutoff, |
| 35 | + nonbondedCutoff=cutoff_distance, |
| 36 | + ) |
| 37 | + |
| 38 | + positions = pdb.getPositions(asNumpy=True) |
| 39 | + |
| 40 | + integrator = openmm.LangevinIntegrator(T, 1 / unit.picosecond, dt) |
| 41 | + |
| 42 | + simulation = app.Simulation(topology, system, integrator) |
| 43 | + simulation.context.setPositions(positions) |
| 44 | + simulation.minimizeEnergy() |
| 45 | + |
| 46 | + return simulation |
| 47 | + |
| 48 | + |
| 49 | +def get_args(argv): |
| 50 | + available_args = [ |
| 51 | + ("k-spring", "k", float, 50, "Spring constant for each replica"), |
| 52 | + ("replicas", "N", int, 25, "Number of replicas along the path"), |
| 53 | + ("time-steps", "t", int, 1e4, "Number of simulation steps for each replica"), |
| 54 | + ("log-period", "l", int, 100, "Frequency of logging the CVs into each histogram"), |
| 55 | + ("log-delay", "d", int, 0, "Number of timesteps to discard before logging"), |
| 56 | + ] |
| 57 | + parser = argparse.ArgumentParser(description="Example script to run umbrella integration") |
| 58 | + for (name, short, T, val, doc) in available_args: |
| 59 | + parser.add_argument("--" + name, "-" + short, type=T, default=T(val), help=doc) |
| 60 | + parser.add_argument("--mpi", action="store_true", help="Use MPI executor") |
| 61 | + args = parser.parse_args(argv) |
| 62 | + return args |
| 63 | + |
| 64 | + |
| 65 | +def get_executor(args): |
| 66 | + if args.mpi: |
| 67 | + futures = importlib.import_module("mpi4py.futures") |
| 68 | + return futures.MPIPoolExecutor() |
| 69 | + return SerialExecutor() |
| 70 | + |
| 71 | + |
| 72 | +def main(argv): |
| 73 | + args = get_args(argv) |
| 74 | + |
| 75 | + cvs = (DihedralAngle((4, 6, 8, 14)), DihedralAngle((6, 8, 14, 16))) |
| 76 | + centers = [] |
| 77 | + center_pos = np.linspace(+0.45 * np.pi, -0.45 * np.pi, args.replicas) |
| 78 | + for pos in center_pos: |
| 79 | + centers.append((pos, pos)) |
| 80 | + method = UmbrellaIntegration(cvs, args.k_spring, centers, args.log_period, args.log_delay) |
| 81 | + |
| 82 | + raw_result = pysages.run( |
| 83 | + method, |
| 84 | + generate_simulation, |
| 85 | + args.time_steps, |
| 86 | + # post_run_action=post_run_action, |
| 87 | + executor=get_executor(args), |
| 88 | + ) |
| 89 | + result = pysages.analyze(raw_result) |
| 90 | + print(result) |
| 91 | + |
| 92 | + return result |
| 93 | + |
| 94 | + |
| 95 | +if __name__ == "__main__": |
| 96 | + main(sys.argv[1:]) |
0 commit comments