def _identifyRestrainedAtoms(self, coordinates, minimum_mass=10.0*units.amu, natoms_to_pick=5, minimum_distance=5.0*units.angstroms, maximum_distance=15.0*units.angstroms, verbose=False): """ Identify atoms to restraint using a simple automated algorithm. ARGUMENTS restraint_filename restraint filename system system object to which restraints are to be added positions reference configuration from which to extract equilibrium distances receptor_atoms atom(s) in receptor (zero-indexed) ligand_atoms atom(s) in ligand (zero-indexed) OPTIONAL ARGUMENTS minimum_mass (simtk.unit.Quantity compatible with simtk.unit.amu) - minimum atomic mass for heavy atoms natoms_to_pick (int) - number of receptor atoms to pick minimum_distance (simtk.unit.Quantity compatible with simtk.unit.angstrom) - minimum distance from ligand atom to receptor atoms to use as restraints maximum_distance (simtk.unit.Quantity compatible with simtk.unit.angstrom) - maximum distance from ligand atom to receptor atoms to use as restraints verbose (boolean) - if True, more debug information is printed during standard execution (default: False) RETURNS restrained_receptor_atoms (set integers) restrained_ligand_atoms (set of integers) NOTES The strategy used here is to first choose a ligand heavy atom near the centroid of the initial ligand geometry. All protein heavy atoms within a given distance range of this designated ligand atom (e.g. between 5 and 15 A) are then identified, and a subset (e.g. 5 atoms) are selected that are far apart from each other. TESTS >>> # Create a test system. >>> import simtk.pyopenmm.extras.testsystems as testsystems >>> [system, coordinates] = testsystems.LysozymeImplicit() >>> # Initialize receptor-ligand restraints. >>> receptor_atoms = range(0,2603) >>> ligand_atoms = range(2603,2621) >>> sigma = 1.0 * units.nanometer >>> temperature = 298.0 * units.kelvin >>> restraints = ReceptorLigandRestraints(system, ligand_atoms, receptor_atoms, temperature) >>> # Identify atoms to be restrained. >>> [restrained_receptor_atoms, restrained_ligand_atoms] = restraints._identifyRestrainedAtoms(coordinates) >>> print restrained_receptor_atoms >>> print restrained_ligand_atoms """ # PARAMETERS minimum_mass = 10.0 * units.amu # minimum atomic mass for heavy atoms natoms_to_pick = 5 # number of receptor atoms to pick minimum_distance = 5.0 * units.angstroms # minimum distance from ligand atom to receptor atoms to use as restraints maximum_distance = 15.0 * units.angstroms # maximum distance from ligand atom to receptor atoms to use as restraints # Compute centroid of ligand. centroid = positions[ligand_atoms,:].mean(axis=0) if verbose: print "Ligand centroid is at (%8.3f, %8.3f, %8.3f) nm\n" % (centroid[0] / units.nanometers, centroid[1] / units.nanometers, centroid[2] / units.nanometers) restrained_ligand_atoms = set() restrained_receptor_atoms = set() # Pick a ligand heavy atom closest to the ligand center of geometry. closest_distance = None for ligand_atom in ligand_atoms: delta = positions[ligand_atom,:] - centroid distance = units.sqrt(delta**2) if ((closest_distance is None) or (distance < closest_distance)) and (system.getParticleMass(ligand_atom) >= minimum_mass): closest_distance = distance closest_atom = ligand_atom restrained_ligand_atoms.add(closest_atom) if verbose: print "Identified ligand atom %d as atom closest to centroid (distance %.3f nm)\n" % (closest_atom, closest_distance / units.nanometers) # First pass: Identify all receptor heavy atoms within specified distance range of chosen ligand atom. atom_pool = list() # pool of all receptor atoms that fit criteria ligand_atom_position = positions[closest_atom,:] # ligand atom position for receptor_atom in receptor_atoms: delta = ligand_atom_position - positions[receptor_atom,:]; distance = units.sqrt(delta**2) if ((distance >= minimum_distance) and (distance < maximum_distance) and (system.getParticleMass(receptor_atom) >= minimum_mass)): atom_pool.append(receptor_atom); # Sanity check: Did we identify the expected number of receptor atoms? if (len(atom_pool) < natoms_to_pick): raise Exception("identifyRestrainedAtoms: Not enough receptor atoms match criterion. Expected at least %d, found only %lu.\n" % (natoms_to_pick, len(atom_pool))) # Identify a given number of receptor heavy atoms that are close to the designated ligand atom but well-separated from each other. These will be restrained. # TODO: Make this less ugly, using simplified routines to find closest and farthest atoms in sets. restrained_receptor_atoms = [atom_index.pop(0)] # list of atoms in receptor where restraints will terminate while (len(restrained_receptor_atoms) < natoms_to_pick): # Find atom in 'atom_pool' that is farthest away from all atoms in 'restrained_receptor_atoms'. farthest_distance = None farthest_atom = None for receptor_atom in atom_pool: closest_distance = None closest_atom = None for restrained_atom in restrained_receptor_atoms: delta = positions[receptor_atom,:] - positions[restrained_atom,:] distance = units.sqrt(delta**2) if ((closest_distance is None) or (distance < closest_distance)): closest_distance = distance closest_atom = restrained_atom if (fathest_distance is None) or (closest_distance > farthest_distance): fathest_distance = closest_distance farthest_atom = receptor_atom # Add the farthest receptor atom to the set of restrained atoms. restrained_receptor_atoms.append(farthest_atom) atom_pool.remove(farthest_atom) if verbose: print "Identified atoms to restrain:" print " receptor: %s\n" % str(restrained_receptor_atoms) print " ligand: %s\n" % str(restrined_ligand_atoms) return [restrained_receptor_atoms, restrained_ligand_atoms] def createRestraintsAutomatically(reference_positions): """ Create restraints automatically based on the given reference positions. ARGUMENTS reference_positions (simtk.unit.Quantity array compatible with simtk.unit.nanometers) - positions to use in determining restraints """ [restrained_receptor_atoms, restrained_ligand_atoms] = self.identifyRestrainedAtoms(reference_positions) return self.imposeRestraints(restrained_receptor_atoms, restrained_ligand_atoms) #============================================================================================= # Harmonic protein-ligand restraint. #============================================================================================= class HarmonicReceptorLigandRestraints(object): """ """ def __init__(self, system, ligand_atoms, receptor_atoms, temperature): ReceptorLigandRestraints.__init__(self, system, ligand_atoms, receptor_atoms, temperature) def getStandardStateCorrection(): """ Compute the standard state correction for a single harmonic distance restraint with nonzero equilibrium length. ARGUMENTS beta (simtk.unit.Quantity) - inverse temperature (1/energy) K (simtk.unit.Quantity) - spring constant (length/distance) r0 (simtk.unit.Quantity) - equilibrium spring distance (length) OPTIONAL ARGUMENTS verbose (bool) - if True, will print extra debug information (default: False) RETURN VALUES DeltaG (float) - computed standard-state correction in dimensionless units (kT); Equivalent to the free energy of releasing restraints and confining into a box of standard state size. EXAMPLE Compute standard state correction for spring with K = 10.0 kcal/mol/nm, r0 = 1.0 nm, at 298 K. >>> temperature = 298.0 * units.kelvin # temperature >>> kB = units.BOLTZMANN_CONSTANT_kB # Boltzmann constant >>> kT = kB * temperature # thermal energy >>> beta = 1.0 / kT # inverse temperature >>> K = 10.0 * units.kilocalories_per_mole / units.nanometers**2 # spring constant >>> r0 = 1.0 * units.nanometers # equilibrium distance >>> DeltaG = _standard_state_correction(beta, K, r0) # compute standard state correction in dimensionless units (kT) >>> print 'DeltaG = %8.3f kT' % DeltaG DeltaG = 4.620 kT """ # Compute standard deviation of equilibrium distribution. # \sigma^2 = 1 / (\beta K) sigma = 1.0 / units.sqrt(beta * K) # standard deviation of equilibrium distribution in 1D if verbose: print "sigma = %f.3 nm" % (sigma / units.nanometers) # Use numerical quadrature to integrate partition function (equal to effective volume of sampled shell). r_min = 0.0 * units.nanometers # initial value for integration r_max = r0 + 10.0 * sigma # maximum radius to integrate to integrand = lambda r : 4.0 * math.pi * r**2 * math.exp(-beta * (K/2.0) * (r*units.nanometers - r0)**2) (shell_volume, shell_volume_error) = scipy.integrate.quad(lambda r : integrand(r), r_min / units.nanometers, r_max / units.nanometers) * units.nanometers**3 # integrate shell volume if verbose: print "shell_volume = %f nm^3" % (shell_volume / units.nanometers**3) # Compute standard-state volume for a single molecule in a box of size (1 L) / (avogadros number) liter = 1000.0 * units.centimeters**3 # one liter box_volume = liter / units.AVOGADRO_NUMBER_NA # standard state volume if verbose: print "box_volume = %f nm^3" % (box_volume / units.nanometers**3) # Compute standard state correction for releasing shell restraints into standard-state box (in units of kT). DeltaG = - math.log(box_volume / shell_volume) if verbose: print "Standard state correction: %.3f kT" % DeltaG # Return standard state correction (in kT). return DeltaG #============================================================================================= # Compute standard state correction. #============================================================================================= @accepts_compatible_units(1.0 / units.kilocalories_per_mole, units.kilocalories_per_mole / units.nanometers**2, units.nanometers) def _standard_state_correction(beta, K, r0, verbose=False): """ Compute the standard state correction for a single harmonic distance restraint with nonzero equilibrium length. ARGUMENTS beta (simtk.unit.Quantity) - inverse temperature (1/energy) K (simtk.unit.Quantity) - spring constant (length/distance) r0 (simtk.unit.Quantity) - equilibrium spring distance (length) OPTIONAL ARGUMENTS verbose (bool) - if True, will print extra debug information (default: False) RETURN VALUES DeltaG (float) - computed standard-state correction in dimensionless units (kT); Equivalent to the free energy of releasing restraints and confining into a box of standard state size. EXAMPLE Compute standard state correction for spring with K = 10.0 kcal/mol/nm, r0 = 1.0 nm, at 298 K. >>> temperature = 298.0 * units.kelvin # temperature >>> kB = units.BOLTZMANN_CONSTANT_kB # Boltzmann constant >>> kT = kB * temperature # thermal energy >>> beta = 1.0 / kT # inverse temperature >>> K = 10.0 * units.kilocalories_per_mole / units.nanometers**2 # spring constant >>> r0 = 1.0 * units.nanometers # equilibrium distance >>> DeltaG = _standard_state_correction(beta, K, r0) # compute standard state correction in dimensionless units (kT) >>> print 'DeltaG = %8.3f kT' % DeltaG DeltaG = 4.620 kT """ # Compute standard deviation of equilibrium distribution. # \sigma^2 = 1 / (\beta K) sigma = 1.0 / units.sqrt(beta * K) # standard deviation of equilibrium distribution in 1D if verbose: print "sigma = %f.3 nm" % (sigma / units.nanometers) # Use numerical quadrature to integrate partition function (equal to effective volume of sampled shell). r_min = 0.0 * units.nanometers # initial value for integration r_max = r0 + 10.0 * sigma # maximum radius to integrate to integrand = lambda r : 4.0 * math.pi * r**2 * math.exp(-beta * (K/2.0) * (r*units.nanometers - r0)**2) (shell_volume, shell_volume_error) = scipy.integrate.quad(lambda r : integrand(r), r_min / units.nanometers, r_max / units.nanometers) * units.nanometers**3 # integrate shell volume if verbose: print "shell_volume = %f nm^3" % (shell_volume / units.nanometers**3) # Compute standard-state volume for a single molecule in a box of size (1 L) / (avogadros number) liter = 1000.0 * units.centimeters**3 # one liter box_volume = liter / units.AVOGADRO_NUMBER_NA # standard state volume if verbose: print "box_volume = %f nm^3" % (box_volume / units.nanometers**3) # Compute standard state correction for releasing shell restraints into standard-state box (in units of kT). DeltaG = - math.log(box_volume / shell_volume) if verbose: print "Standard state correction: %.3f kT" % DeltaG # Return standard state correction (in kT). return DeltaG #============================================================================================= # Method for creating receptor-ligand restraints. #============================================================================================= @accepts_compatible_units(None, units.nanometers, None, None, units.nanometers, units.kelvin, None, None) def create_binding_restraints(reference_system, reference_positions, receptor_atoms, ligand_atoms, sigma, temperature, verbose=False, mm=simtk.chem.openmm): """ Automatically choose restraints between receptor and ligand. ARGUMENTS reference_system (pyopenmm.System) - OpenMM System object containing receptor and ligand reference_positions (numpy array [natoms x 3]) - reference coordinates to use for selecting atoms involved in restraints receptor_atoms (Python list of int) - list of atoms in receptor ligand_atoms (Python list of int) - list of atoms in ligand sigma (float) - desired standard deviation to allow (in nm) tempreature (float) - desired temperature (in K) OPTIONAL ARGUMENTS verbose (boolean) - if True, detailed information about the restraint selection process will be emitted (default: False) mm - OpenMM implementation (default:simtk.chem.openmm) RETURNS system (pyopenmm.System) - System object with restraints added between receptor and ligand NOTES A simple approach to adding protein-ligand restraints to ensure the ligand stays near the binding pocket at intermediate alchemical states, while being able to explore alternative orientations and binding geometries when fully decoupled. The strategy used here is to first choose a ligand heavy atom near the centroid of the initial ligand geometry. All protein heavy atoms within a given distance range of this designated ligand atom (e.g. between 5 and 15 A) are then identified, and a subset (e.g. 5 atoms) are selected that are far apart from each other. Harmonic restraints are then applied between all such designated receptor atoms and the designated ligand atom with an equilibrium length equal to the initial heavy atom-ligand atom separation, and a spring constant that is automatically computed so as to provide an approximate standard deviation (when the ligand is noninteracting with the protein) specified by the user. Because the designated protein atoms are more numerous than three and far apart from each other, this should define a single 'soft pocket' for the ligand to reside in, while still allowing reorientation. Note that these restraints bind together parts of the protein, so they must be turned off after the ligand is made noninteracting, and an additional spring added to ensure that the ligand occupies a known volume and an analytical standard-state correction can be applied. TODO * Is this approach optimal, or even a good idea? Let's think about this some more. * Use Python unit. FUTURE * Let's make this more general and allow us to 'plug in' different schemes for identifying and imposing restraints. @author John D. Chodera EXAMPLES Create restraints between atoms in a simple test system. >>> import simtk.pyopenmm.extras.testsystems as testsystems >>> [system, coordinates] = testsystems.LysozymeImplicit() >>> receptor_atoms = range(0,2603) >>> ligand_atoms = range(2603,2621) >>> sigma = 1.0 * units.nanometer >>> temperature = 298.0 * units.kelvin >>> create_binding_restraints(system, coordinates, receptor_atoms, ligand_atoms, sigma, temperature) """ # Compute spring constant from desired standard deviation. # We want an effective spring constant of K = 1 / (\beta \sigma^2), and there will be n springs, so we use K = 1 / (n \beta \sigma^2). kB = units.BOLTZMANN_CONSTANT_kB kT = (kB * temperature) # thermal energy, kJ/mol beta = 1.0 / kT # inverse temperature, 1/(kJ/mol) K = 1.0 / (len(restrained_receptor_atoms) * len(restrained_ligand_atoms) * beta * sigma**2); # Add restrants. # TODO: Are we allowed only one HarmonicBondForce per System? bonds = mm.HarmonicBondForce() if verbose: print "%16s %16s %24s %24s\n" % ("rec atm", "lig atm", "r0 (nm)", "K (kcal/mol/nm^2)") for receptor_atom in restrained_receptor_atoms: for ligand_atom in restrained_ligand_atoms: delta = coordinates[receptor_atom,:] - coordinates[ligand_atom,:] r0 = units.sqrt(delta**2) bonds.addBond(receptor_atom, ligand_atom, r0, K) if verbose: print "%16d %16d %24.3f %24.3f\n" % (receptor_atom, ligand_atom, r0, K); system.addForce(bonds) # Compute and return standard-state correction. DeltaG = _standard_state_correction(beta, K, r0, verbose=verbose) return DeltaG @accepts_compatible_units(None, None, units.nanometers, None, None, None) def identify_restrained_atoms(restraint_filename, system, positions, receptor_atoms, ligand_atoms, verbose=False): """ Pick atoms in ligand and receptor to be restrained. The strategy used here is to first choose a ligand heavy atom near the centroid of the initial ligand geometry. All protein heavy atoms within a given distance range of this designated ligand atom (e.g. between 5 and 15 A) are then identified, and a subset (e.g. 5 atoms) are selected that are far apart from each other. TODO: Return the restrained_receptor_atoms and restraint_ligand_atoms as a pair<> type instead of making them arguments? ARGUMENTS restraint_filename restraint filename system system object to which restraints are to be added positions reference configuration from which to extract equilibrium distances receptor_atoms atom(s) in receptor (zero-indexed) ligand_atoms atom(s) in ligand (zero-indexed) OPTIONAL ARGUMENTS verbose (boolean) - if True, more debug information is printed during standard execution (default: False) RETURNS restrained_receptor_atoms (set integers) restrained_ligand_atoms (set of integers) EXAMPLES Pick atoms to be restrained in a simple test system. >>> # Create a test system. >>> import simtk.pyopenmm.extras.testsystems as testsystems >>> [system, coordinates] = testsystems.LysozymeImplicit() >>> # Create temporary file for storing output. >>> import tempfile >>> file = tempfile.NamedTemporaryFile() # temporary file for testing >>> restraint_filename = file.name >>> # Identify atoms to be restrained. >>> receptor_atoms = range(0,2603) >>> ligand_atoms = range(2603,2621) >>> sigma = 1.0 * units.nanometer >>> temperature = 298.0 * units.kelvin >>> [restrained_receptor_atoms, restrained_ligand_atoms] = create_binding_restraints(system, coordinates, receptor_atoms, ligand_atoms, sigma, temperature) """ # PARAMETERS minimum_mass = 10.0 * units.amu # minimum atomic mass for heavy atoms natoms_to_pick = 5 # number of receptor atoms to pick minimum_distance = 5.0 * units.angstroms # minimum distance from ligand atom to receptor atoms to use as restraints maximum_distance = 15.0 * units.angstroms # maximum distance from ligand atom to receptor atoms to use as restraints # Compute centroid of ligand. centroid = positions[ligand_atoms,:].mean(axis=0) if verbose: print "Ligand centroid is at (%8.3f, %8.3f, %8.3f) nm\n" % (centroid[0] / units.nanometers, centroid[1] / units.nanometers, centroid[2] / units.nanometers) restrained_ligand_atoms = set() restrained_receptor_atoms = set() # Pick a ligand heavy atom closest to the ligand center of geometry. closest_distance = None for ligand_atom in ligand_atoms: delta = positions[ligand_atom,:] - centroid distance = units.sqrt(delta**2) if ((closest_distance is None) or (distance < closest_distance)) and (system.getParticleMass(ligand_atom) >= minimum_mass): closest_distance = distance closest_atom = ligand_atom restrained_ligand_atoms.add(closest_atom) if verbose: print "Identified ligand atom %d as atom closest to centroid (distance %.3f nm)\n" % (closest_atom, closest_distance / units.nanometers) # First pass: Identify all receptor heavy atoms within specified distance range of chosen ligand atom. atom_pool = list() # pool of all receptor atoms that fit criteria ligand_atom_position = positions[closest_atom,:] # ligand atom position for receptor_atom in receptor_atoms: delta = ligand_atom_position - positions[receptor_atom,:]; distance = units.sqrt(delta**2) if ((distance >= minimum_distance) and (distance < maximum_distance) and (system.getParticleMass(receptor_atom) >= minimum_mass)): atom_pool.append(receptor_atom); # Sanity check: Did we identify the expected number of receptor atoms? if (len(atom_pool) < natoms_to_pick): raise Exception("identifyRestrainedAtoms: Not enough receptor atoms match criterion. Expected at least %d, found only %lu.\n" % (natoms_to_pick, len(atom_pool))) # Identify a given number of receptor heavy atoms that are close to the designated ligand atom but well-separated from each other. These will be restrained. # TODO: Make this less ugly, using simplified routines ot find closest and farthest atoms in sets. restrained_receptor_atoms = [atom_index.pop(0)] # list of atoms in receptor where restraints will terminate while (len(restrained_receptor_atoms) < natoms_to_pick): # Find atom in 'atom_pool' that is farthest away from all atoms in 'restrained_receptor_atoms'. farthest_distance = None farthest_atom = None for receptor_atom in atom_pool: closest_distance = None closest_atom = None for restrained_atom in restrained_receptor_atoms: delta = positions[receptor_atom,:] - positions[restrained_atom,:] distance = units.sqrt(delta**2) if ((closest_distance is None) or (distance < closest_distance)): closest_distance = distance closest_atom = restrained_atom if (fathest_distance is None) or (closest_distance > farthest_distance): fathest_distance = closest_distance farthest_atom = receptor_atom # Add the farthest receptor atom to the set of restrained atoms. restrained_receptor_atoms.append(farthest_atom) atom_pool.remove(farthest_atom) if verbose: print "Identified receptor atoms to restrain:\n" for atom in restrained_receptor_atoms: print "%8d" % atom, print "" return @accepts_compatible_units(None, None, units.nanometers, None, None, units.nanometers, units.kelvin, None) def create_restraints(restraint_filename, system, positions, receptor_atoms, ligand_atoms, sigma, temperature, verbose=False): """ Create receptor-ligand distance restraints to keep ligand in binding pocket. If a file exists specifying which atoms to restrain, use these. Otherwise, identify receptor atoms and a ligand atom and compute standard-state correction. ARGUMENTS restraint_filename restraint filename system system object to which restraints are to be added positions reference configuration from which to extract equilibrium distances receptor_atoms atom(s) in receptor (zero-indexed) ligand_atoms atom(s) in ligand (zero-indexed) sigma the desired standard deviation to allow, for computation of spring constants (nm) temperature temperature (K) EXAMPLES Pick atoms to be restrained in a simple test system. >>> # Create a test system. >>> import simtk.pyopenmm.extras.testsystems as testsystems >>> [system, coordinates] = testsystems.LysozymeImplicit() >>> # Create temporary file for storing output. >>> import tempfile >>> file = tempfile.NamedTemporaryFile() # temporary file for testing >>> restraint_filename = file.name >>> # Identify atoms to be restrained. >>> receptor_atoms = range(0,2603) >>> ligand_atoms = range(2603,2621) >>> sigma = 1.0 * units.nanometer >>> temperature = 298.0 * units.kelvin >>> create_restraints(restraint_filename, system, coordinates, receptor_atoms, ligand_atoms, sigma, temperature) """ # Storage for restrained atoms. restrained_receptor_atoms = set() restrained_ligand_atoms = set() standard_state_correction = None # Attempt to open the restraint file. restraint_file = open(restraint_filename, 'r') if (restraint_file is None): # Restraint file does not exist. Pick restraints and write these to restraint file. if verbose: print "Restraint file does not exist. Will create it." # Pick atoms to restrain, returning standard state correction. identify_restrained_atoms(restraint_filename, system, positions, receptor_atoms, ligand_atoms, restrained_receptor_atoms, restrained_ligand_atoms) # Impose restraints, and store computed standard-state correction for retaining first distance restraint. standard_state_correction = impose_harmonic_restraints(system, positions, restrained_receptor_atoms, restrained_ligand_atoms, sigma, temperature) # Write restraints to file. # Open file. restraint_file = open(restraint_filename, 'w') # write number of restrained ligand atoms on a line restraint_file.write("%8d\n" % (len(restrained_ligand_atoms))) # write restrained ligand atoms on a line for restrained_atom in restrained_ligand_atoms: restraint_file.write('%8d' % restrained_atom) restraint_file.write("\n") # write number of restrained receptor atoms on a line fprintf(restraint_file, "%8lu\n", (long unsigned)(restrained_receptor_atoms.size())); // write restrained receptor atoms on a line for(std::set::iterator restrained_atom = restrained_receptor_atoms.begin(); restrained_atom != restrained_receptor_atoms.end(); restrained_atom++) fprintf(restraint_file, "%8d", *restrained_atom); fprintf(restraint_file, "\n"); // write standard-state correction fprintf(restraint_file, "%16.8f\n", standard_state_correction); // close the file fclose(restraint_file); else: # Restraint file exists. Read restraints from file. if verbose: print "Restraint file exists. Reading it." # Read the restraint file. int nrestrained, restrained_atom; // read number of restrained ligand atoms fscanf(restraint_file, "%8d\n", &nrestrained); // read restrained ligand atoms for(int restrained_count = 0; restrained_count < nrestrained; restrained_count++) { fscanf(restraint_file, "%8d", &restrained_atom); restrained_ligand_atoms.insert(restrained_atom); } // read number of restrained receptor atoms fscanf(restraint_file, "%8d\n", &nrestrained); // read restrained receptor for(int restrained_count = 0; restrained_count < nrestrained; restrained_count++) { fscanf(restraint_file, "%8d", &restrained_atom); restrained_receptor_atoms.insert(restrained_atom); } // Read standard-state correction. fscanf(restraint_file, "%16.8lf\n", &standard_state_correction); // close the file fclose(restraint_file); // Impose restraints, and store computed standard-state correction for retaining first distance restraint. standard_state_correction = impose_harmonic_restraints(system, positions, restrained_receptor_atoms, restrained_ligand_atoms, sigma, temperature); } if (verbose) { printf("Standard state correction: %.3f kT\n", standard_state_correction); printf("Restrained ligand atoms:\n"); for(std::set::iterator atom = restrained_ligand_atoms.begin(); atom != restrained_ligand_atoms.end(); atom++) printf("%8d", *atom); printf("\n"); printf("Restrained receptor atoms:\n"); for(std::set::iterator atom = restrained_receptor_atoms.begin(); atom != restrained_receptor_atoms.end(); atom++) printf("%8d", *atom); printf("\n"); } } return def CreateFlatBottomedRestraints(): """ """ ned_receptor_atoms, std::set & restrained_ligand_atoms) { bool verbose = true; OpenMM::CustomBondForce* flatbonds = new OpenMM::CustomBondForce("step(r-r0)*scale*K*(r-r0)^2"); flatbonds->addPerBondParameter("r0"); flatbonds->addPerBondParameter("K"); flatbonds->addGlobalParameter("scale",0.5); // Calculate r0, which is half of the longest distance in protein plus 5 A. double r00 = 0.00; OpenMM::Vec3 distance(0.0, 0.0, 0.0); for(std::set< int >::iterator receptor_atom1 = receptor_atoms.begin(); receptor_atom1 != receptor_atoms.end(); receptor_atom1++){ for(std::set< int >::iterator receptor_atom2 = receptor_atoms.begin(); receptor_atom2 != receptor_atoms.end(); receptor_atom2++){ \ distance=positions[*receptor_atom2]-positions[*receptor_atom1]; double Distance = sqrt(distance.dot(distance)); if (Distance > r00) r00=Distance; } } double r0 = r00/2.0 + 0.5; //double r0 = 2.5; double K = 2.0 * 0.0083144621 * 5.0 * 298.0 *100 ; double standard_state_correction(double beta, double K, double r0) { // Compute standard deviation of equilibrium distribution. // \sigma^2 = 1 / (\beta K) double sigma = 1.0 / sqrt(beta * K); // Use trapezoidal quadrature to integrate. double dr = sigma / 10000.0; double r_max = r0 + 10.0 * sigma; double shell_volume = 0.0; double r = 0.0; shell_volume += (dr/2.0) * 4.0 * M_PI * (r*r) * exp(-beta * (K/2.0) * (r-r0)*(r-r0)); while(r < r_max) { shell_volume += dr * 4.0 * M_PI * (r*r) * exp(-beta * (K/2.0) * (r-r0)*(r-r0)); r += dr; } shell_volume += (dr/2.0) * 4.0 * M_PI * (r*r) * exp(-beta * (K/2.0) * (r-r0)*(r-r0)); // Compute standard-state volume for a single molecule in a box of size (1 L) / (avogadros number) double NA = 6.0221415e23; // Avogadro's number double liter = 1.0e24; // one liter, nm^3 double box_volume = liter / NA; // standard state, nm^3 printf("shell_volume = %f\n", shell_volume); // Compute standard state correction for releasing shell restraints into standard-state box in units of kT. double DeltaG = - log(box_volume / shell_volume); printf("Standard state correction: %.3f kT\n", DeltaG); // Return standard state correction. return DeltaG; }