/* ----------------------------------------------------------------------------- * OpenMM(tm) HelloEthane example (May 2009) * ----------------------------------------------------------------------------- * This is a complete, self-contained "hello world" example demonstrating * GPU-accelerated simulation of a system with both bonded and nonbonded forces, * using ethane (H3-C-C-H3) in a vacuum as an example. A multi-frame PDB file is * written to stdout which can be read by VMD or other visualization tool to * produce an animation of the resulting trajectory. * * Pay particular attention to the handling of units in this example. Incorrect * handling of units is a very common error; this example shows how you can * continue to work with Amber-style units of Angstroms and kCals while correctly * communicating with OpenMM in nanometers and kJoules. * -------------------------------------------------------------------------- */ // Suppress irrelevant warnings from Microsoft's compiler. #ifdef _MSC_VER #pragma warning(disable:4996) // sprintf is unsafe #pragma warning(disable:4251) // no dll interface for some classes #endif #include "OpenMM.h" #include #include #include #include using namespace OpenMM; // ----------------------------------------------------------------------------- // MODELING AND SIMULATION PARAMETERS // ----------------------------------------------------------------------------- const bool UseConstraints = false; // Should we constrain C-H bonds? const double StepSizeInFs = 2; // integration step size (fs) const double ReportIntervalInFs = 10; // how often to generate PDB frame (fs) const double SimulationTimeInPs = 100; // total simulation time (ps) static void simulateEthane(); static void writePDB(const OpenMMContext&); // PDB file writer; see below. // ----------------------------------------------------------------------------- // MAIN PROGRAM // ----------------------------------------------------------------------------- int main() { // ALWAYS enclose all OpenMM calls with a try/catch block to make sure that // usage and runtime errors are caught and reported. try { // Load all available OpenMM plugins from their default location. Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory()); simulateEthane(); return 0; // Normal return from main. } // Catch and report usage and runtime errors detected by OpenMM and fail. catch(const std::exception& e) { printf("EXCEPTION: %s\n", e.what()); return 1; } } // ----------------------------------------------------------------------------- // FORCE FIELD DATA // ----------------------------------------------------------------------------- // These data structures are not part of OpenMM; they are a model of the kinds // of data structures an MD code uses to hold a set of force field parameters. // For this example we're using a tiny subset of the Amber99 force field. // We want to keep the data in the original unit system to avoid conversion // bugs; this requires conversion on the way in and out of OpenMM. // Amber reduces nonbonded forces between 1-4 bonded atoms. const double Coulomb14Scale = 0.5; const double LennardJones14Scale = 0.5; struct AtomType { double mass, charge, vdwRadiusInAngstroms, vdwEnergyInKcal; } atomType[] = {/*0 H*/ 1.008, 0.0605, 1.4870, 0.0157, /*1 C*/12.011, -.1815, 1.9080, 0.1094}; const int H = 0, C = 1; struct BondType { double nominalLengthInAngstroms, stiffnessInKcalPerAngstrom2; bool canConstrain; } bondType[] = {/*0 CC*/1.526, 310., false, /*1 CH*/1.09 , 340., true}; const int CC = 0, CH = 1; struct AngleType { double nominalAngleInDegrees, stiffnessInKcalPerRadian2; } angleType[] = {/*0 HCC*/109.5, 50., /*1 HCH*/109.5, 35.}; const int HCC = 0, HCH = 1; struct TorsionType { int periodicity; double phaseInDegrees, amplitudeInKcal; } torsionType[] = {/*0 HCCH*/3, 0., 0.150}; const int HCCH = 0; // ----------------------------------------------------------------------------- // MOLECULE DATA // ----------------------------------------------------------------------------- // Now describe an ethane molecule by listing its atoms, bonds, angles, and // torsions. We'll provide a default configuration which centers the molecule // at (0,0,0) with the C-C bond along the x axis. // Use this as an "end of list" marker so that we do not have to count; let the // computer do that! const int EndOfList=-1; struct AtomInfo { int type; const char* pdb; Vec3 initPosInAngstroms;} atoms[] = {/*0*/C, " C1 ", Vec3( -.7605, 0, 0 ), /*1*/C, " C2 ", Vec3( .7605, 0, 0 ), /*2*/H, "1H1 ", Vec3(-1.135, 1.03, 0 ), // bonded to C1 /*3*/H, "2H1 ", Vec3(-1.135, -.51, .89), /*4*/H, "3H1 ", Vec3(-1.135, -.51,-.89), /*5*/H, "1H2 ", Vec3( 1.135, 1.03, 0 ), // bonded to C2 /*6*/H, "2H2 ", Vec3( 1.135, -.51, .89), /*7*/H, "3H2 ", Vec3( 1.135, -.51,-.89), EndOfList}; static struct {int type; int atoms[2];} bonds[] = {CC,0,1, CH,0,2,CH,0,3,CH,0,4, // C1 methyl CH,1,5,CH,1,6,CH,1,7, // C2 methyl EndOfList}; static struct {int type; int atoms[3];} angles[] = {HCC,2,0,1,HCC,3,0,1,HCC,4,0,1, // C1 methyl HCH,2,0,3,HCH,2,0,4,HCH,3,0,4, HCC,5,1,0,HCC,6,1,0,HCC,7,1,0, // C2 methyl HCH,5,1,6,HCH,5,1,7,HCH,6,1,7, EndOfList}; static struct {int type; int atoms[4];} torsions[] = {HCCH,2,0,1,5,HCCH,2,0,1,6,HCCH,2,0,1,7, HCCH,3,0,1,5,HCCH,3,0,1,6,HCCH,3,0,1,7, HCCH,4,0,1,5,HCCH,4,0,1,6,HCCH,4,0,1,7, EndOfList}; // Add missing scalar product operators for Vec3. Vec3 operator*(const Vec3& v, double r) {return Vec3(v[0]*r, v[1]*r, v[2]*r);} Vec3 operator*(double r, const Vec3& v) {return Vec3(r*v[0], r*v[1], r*v[2]);} // This is the conversion factor that takes you from a van der Waals radius // (defined as 1/2 the minimum energy separation) to the related Lennard Jones // "sigma" parameter (defined as the zero crossing separation). static const double SigmaPerVdwRadius = 2*std::pow(2., -1./6.); // ----------------------------------------------------------------------------- // ETHANE SIMULATION // ----------------------------------------------------------------------------- static void simulateEthane() { // ------------------------------------------------------------------------- // Create a System and Force objects within the System. Retain a reference // to each force object so we can fill in the forces. Note: OpenMM owns // the objects and will take care of deleting them; don't do it yourself! // ------------------------------------------------------------------------- System system; NonbondedForce& nonbond = *new NonbondedForce(); HarmonicBondForce& bondStretch = *new HarmonicBondForce(); HarmonicAngleForce& bondBend = *new HarmonicAngleForce(); PeriodicTorsionForce& bondTorsion = *new PeriodicTorsionForce(); system.addForce(&nonbond); system.addForce(&bondStretch); system.addForce(&bondBend); system.addForce(&bondTorsion); // ------------------------------------------------------------------------- // Specify the atoms and their properties: // (1) System needs to know the masses. // (2) NonbondedForce needs charges,van der Waals properties (in MD units!). // (3) Collect default positions for initializing the simulation later. // ------------------------------------------------------------------------- std::vector initialPositions; for (int n=0; atoms[n].type != EndOfList; ++n) { const AtomType& atype = atomType[atoms[n].type]; system.addParticle(atype.mass); nonbond.addParticle(atype.charge, atype.vdwRadiusInAngstroms * NmPerAngstrom * SigmaPerVdwRadius, atype.vdwEnergyInKcal * KJPerKcal); initialPositions.push_back(atoms[n].initPosInAngstroms * NmPerAngstrom); } // ------------------------------------------------------------------------- // Process the bonds: // (1) HarmonicBondForce needs bond stretch parameters (in MD units!). // (2) If we're using constraints, tell System about constrainable bonds. // (3) Create a list of bonds for generating nonbond exclusions. // ------------------------------------------------------------------------- std::vector< std::pair > bondPairs; for (int i=0; bonds[i].type != EndOfList; ++i) { const int* atom = bonds[i].atoms; const BondType& bond = bondType[bonds[i].type]; // Note factor of 2 for stiffness below because Amber specifies the constant // as it is used in the harmonic energy term kx^2 with force 2kx; OpenMM wants // it as used in the force term kx, with energy kx^2/2. bondStretch.addBond(atom[0], atom[1], bond.nominalLengthInAngstroms * NmPerAngstrom, bond.stiffnessInKcalPerAngstrom2 * 2 * KJPerKcal * AngstromsPerNm * AngstromsPerNm); if (UseConstraints && bond.canConstrain) system.addConstraint(atom[0], atom[1], bond.nominalLengthInAngstroms * NmPerAngstrom); bondPairs.push_back(std::make_pair(atom[0], atom[1])); } // Exclude 1-2, 1-3 bonded atoms from nonbonded forces, and scale down 1-4 bonded atoms. nonbond.createExceptionsFromBonds(bondPairs, Coulomb14Scale, LennardJones14Scale); // ------------------------------------------------------------------------- // Create the 1-2-3 bond angle harmonic terms. // ------------------------------------------------------------------------- for (int i=0; angles[i].type != EndOfList; ++i) { const int* atom = angles[i].atoms; const AngleType& angle = angleType[angles[i].type]; // See note under bond stretch above regarding the factor of 2 here. bondBend.addAngle(atom[0],atom[1],atom[2], angle.nominalAngleInDegrees * RadiansPerDegree, angle.stiffnessInKcalPerRadian2 * 2 * KJPerKcal); } // ------------------------------------------------------------------------- // Create the 1-2-3-4 bond torsion (dihedral) terms. // ------------------------------------------------------------------------- for (int i=0; torsions[i].type != EndOfList; ++i) { const int* atom = torsions[i].atoms; const TorsionType& torsion = torsionType[torsions[i].type]; bondTorsion.addTorsion(atom[0],atom[1],atom[2],atom[3], torsion.periodicity, torsion.phaseInDegrees * RadiansPerDegree, torsion.amplitudeInKcal * KJPerKcal); } // ------------------------------------------------------------------------- // Choose an Integrator for advancing time, and a Context connecting the // System with the Integrator for simulation. Let the Context choose the // best available Platform. Initialize the configuration from the default // positions we collected above. Initial velocities will be zero. // ------------------------------------------------------------------------- VerletIntegrator integrator(StepSizeInFs * PsPerFs); OpenMMContext context(system, integrator); context.setPositions(initialPositions); // ------------------------------------------------------------------------- // Run the simulation: // (1) Write the first line of the PDB file and the initial configuration. // (2) Run silently entirely within OpenMM between reporting intervals. // (3) Write a PDB frame when the time comes. // ------------------------------------------------------------------------- printf("REMARK Using OpenMM platform %s\n", context.getPlatform().getName().c_str() ); writePDB(context); const int NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5); do { integrator.step(NumSilentSteps); writePDB(context); } while (context.getTime() < SimulationTimeInPs); } // ----------------------------------------------------------------------------- // PDB FILE WRITER // ----------------------------------------------------------------------------- static void writePDB(const OpenMMContext& context) { // Caution: at the moment asking for energy requires use of slow Reference // platform calculation. const State state = context.getState(State::Positions | State::Velocities | State::Energy); const double energy = state.getPotentialEnergy() + state.getKineticEnergy(); const std::vector& positions = state.getPositions(); static int modelFrameNumber = 0; // numbering for MODEL records in pdb output modelFrameNumber++; printf("MODEL %d\n", modelFrameNumber); printf("REMARK 250 time=%.3f picoseconds; Energy = %.3f kilojoules/mole\n", state.getTime(), energy); for (unsigned i=0; i < positions.size(); ++i) { const Vec3 pos = positions[i] * AngstromsPerNm; printf("ATOM %5d %4s ETH 1 %8.3f%8.3f%8.3f 1.00 0.00 \n", i+1, atoms[i].pdb, pos[0], pos[1], pos[2]); } printf("ENDMDL\n"); }