diff --git a/gamd/config.py b/gamd/config.py index ff1d2ab..3a89411 100644 --- a/gamd/config.py +++ b/gamd/config.py @@ -22,7 +22,7 @@ def assign_tag(root, tagname, value, attributes=None): xmlTag.text = str(value) for attribute in attributes: xmlTag.set(attribute, attributes[attribute]) - + return class SystemConfig: @@ -30,12 +30,16 @@ def __init__(self): self.nonbonded_method = "PME" self.nonbonded_cutoff = 0.9*unit.nanometer self.constraints = "HBonds" + self.switch_distance = 1.0*unit.nanometer + self.ewald_error_tolerance = 0.0005 return - + def serialize(self, root): assign_tag(root, "nonbonded-method", self.nonbonded_method) assign_tag(root, "nonbonded-cutoff", self.nonbonded_cutoff.value_in_unit(unit.nanometers)) assign_tag(root, "constraints", self.constraints) + assign_tag(root, "switch-distance", self.switch_distance.value_in_unit(unit.nanometers)) + assign_tag(root, "ewald-error-tolerance", self.ewald_error_tolerance) return class BarostatConfig: @@ -43,7 +47,7 @@ def __init__(self): self.pressure = 1.0 * unit.bar self.frequency = 25 return - + def serialize(self, root): assign_tag(root, "pressure", self.pressure.value_in_unit(unit.bar)) assign_tag(root, "frequency", self.frequency) @@ -54,7 +58,7 @@ def __init__(self): self.primary = 6.0 * unit.kilocalories_per_mole self.secondary = 6.0 * unit.kilocalories_per_mole return - + def serialize(self, root): assign_tag(root, "primary", self.primary.value_in_unit(unit.kilocalories_per_mole)) assign_tag(root, "secondary", self.secondary.value_in_unit(unit.kilocalories_per_mole)) @@ -70,7 +74,7 @@ def __init__(self): self.total_simulation_length = 0 self.averaging_window_interval = 0 return - + def serialize(self, root): assign_tag(root, "conventional-md-prep", self.conventional_md_prep) assign_tag(root, "conventional-md", self.conventional_md) @@ -80,7 +84,7 @@ def serialize(self, root): #assign_tag(root, "total-simulation-length", self.total_simulation_length) assign_tag(root, "averaging-window-interval", self.averaging_window_interval) return - + def compute_total_simulation_length(self): self.total_simulation_length = self.conventional_md \ + self.gamd_equilibration + self.gamd_production @@ -95,7 +99,7 @@ def __init__(self): self.friction_coefficient = 1.0 * unit.picoseconds ** -1 self.number_of_steps = IntegratorNumberOfStepsConfig() return - + def serialize(self, root): assign_tag(root, "algorithm", self.algorithm) assign_tag(root, "boost-type", self.boost_type) @@ -114,35 +118,59 @@ def __init__(self): self.coordinates = "" self.coordinates_filetype = "" return - + def serialize(self, root): assign_tag(root, "topology", self.topology) assign_tag(root, "coordinates", self.coordinates, {"type": self.coordinates_filetype}) return - class CharmmConfig: def __init__(self): self.topology = "" self.coordinates = "" + self.coordinates_filetype = "" self.parameters = [] + self.box_vectors = [] + self.alpha = 90 * unit.degree + self.beta = 90 * unit.degree + self.gamma = 90 * unit.degree + self.config_box_vector_defined = False return - + def serialize(self, root): assign_tag(root, "topology", self.topology) - assign_tag(root, "coordinates", self.coordinates) - xmlParams = ET.SubElement(root, 'parameters') - for params_filename in self.parameters: - assign_tag(xmlParams, "file", params_filename) + assign_tag(root, "coordinates", self.coordinates, {"type": self.coordinates_filetype}) + assign_tag(root, "parameters", self.parameters) + assign_tag(root, "box_vector", self.box_vectors) return + def get_box_vectors(self): + errorMessage = "The box vectors were only partially defined. " \ + "Box lengths 'a', 'b', and 'c' must be present, along with "\ + "box angles 'alpha', 'beta' and 'gamma'." + if not self.is_whole_box_defined(): + raise RuntimeError(errorMessage) + return [self.box_vectors["a"], self.box_vectors["b"], + self.box_vectors["c"], + self.alpha, self.beta, self.gamma] + + def is_whole_box_defined(self): + box_props = ["a", "b", "c", "alpha", "beta", "gamma"] + result = True + for property in box_props: + result = result and property in self.box_vectors + return result + + def is_config_box_vector_defined(self): + return self.config_box_vector_defined + class GromacsConfig: def __init__(self): self.topology = "" self.coordinates = "" self.include_dir = "" return - + def serialize(self, root): assign_tag(root, "topology", self.topology) assign_tag(root, "coordinates", self.coordinates) @@ -155,7 +183,7 @@ def __init__(self): self.forcefield_list_native = [] self.forcefield_list_external = [] return - + def serialize(self, root): assign_tag(root, "coordinates", self.coordinates) xmlForcefields = ET.SubElement(root, "forcefields") @@ -174,7 +202,7 @@ def __init__(self): self.gromacs = None self.forcefield = None return - + def serialize(self, root): if self.amber is not None: xml_amber_tags = ET.SubElement(root, "amber") @@ -198,13 +226,13 @@ def __init__(self): self.restart_checkpoint_interval = 50000 self.statistics_interval = 500 return - + def compute_chunk_size(self): gcd = np.gcd.reduce( [self.energy_interval, self.coordinates_interval, self.restart_checkpoint_interval, self.statistics_interval]) return gcd - + def serialize(self, root): xml_energy_tags = ET.SubElement(root, "energy") assign_tag(xml_energy_tags, "interval", self.energy_interval) @@ -220,7 +248,7 @@ def __init__(self): self.overwrite_output = True self.reporting = OutputsReportingConfig() return - + def serialize(self, root): assign_tag(root, "directory", self.directory) assign_tag(root, "overwrite-output", self.overwrite_output) @@ -231,7 +259,7 @@ def serialize(self, root): class Config: def __init__(self): # set all the default values - + self.temperature = 298.15 * unit.kelvin self.system = SystemConfig() self.barostat = None #BarostatConfig() @@ -239,7 +267,7 @@ def __init__(self): self.integrator = IntegratorConfig() self.input_files = InputFilesConfig() self.outputs = OutputsConfig() - + def serialize(self, filename): root = ET.Element('gamd') assign_tag(root, "temperature", self.temperature.value_in_unit(unit.kelvin)) @@ -255,13 +283,13 @@ def serialize(self, filename): self.input_files.serialize(xml_input_files) xml_outputs = ET.SubElement(root, "outputs") self.outputs.serialize(xml_outputs) - + xmlstr = minidom.parseString(ET.tostring(root)).toprettyxml( indent=" ") our_file=open(filename, 'w') our_file.write(xmlstr) our_file.close() - + if __name__ == "__main__": - pass \ No newline at end of file + pass diff --git a/gamd/gamdSimulation.py b/gamd/gamdSimulation.py index 28bdbf4..8467bdb 100644 --- a/gamd/gamdSimulation.py +++ b/gamd/gamdSimulation.py @@ -1,7 +1,7 @@ ''' Created on Oct 28, 2020 -Creates all OpenMM objects from Config() object that can be used in a +Creates all OpenMM objects from Config() object that can be used in a GaMD simulation. @author: lvotapka @@ -32,7 +32,7 @@ def load_pdb_positions_and_box_vectors(pdb_coords_filename, need_box): "found in {}. ".format(pdb_coords_filename) \ + "Box vectors for an anchor must be defined with a CRYST "\ "line within the PDB file." - + return positions, pdb_parmed.box_vectors @@ -53,53 +53,53 @@ def __init__(self): class GamdSimulationFactory: def __init__(self): return - + def createGamdSimulation(self, config, platform_name, device_index): need_box = True if config.system.nonbonded_method == "pme": nonbondedMethod = openmm_app.PME - + elif config.system.nonbonded_method == "nocutoff": nonbondedMethod = openmm_app.NoCutoff need_box = False - + elif config.system.nonbonded_method == "cutoffnonperiodic": nonbondedMethod = openmm_app.CutoffNonPeriodic - + elif config.system.nonbonded_method == "cutoffperiodic": nonbondedMethod = openmm_app.CutoffPeriodic - + elif config.system.nonbonded_method == "ewald": nonbondedMethod = openmm_app.Ewald - + else: - raise Exception("nonbonded method not found: %s", + raise Exception("nonbonded method not found: %s", config.system.nonbonded_method) - + if config.system.constraints == "none" \ or config.system.constraints is None: constraints = None - + elif config.system.constraints == "hbonds": constraints = openmm_app.HBonds - + elif config.system.constraints == "allbonds": constraints = openmm_app.AllBonds - + elif config.system.constraints == "hangles": constraints = openmm_app.HAngles - + else: - raise Exception("constraints not found: %s", + raise Exception("constraints not found: %s", config.system.constraints) - + box_vectors = None gamdSimulation = GamdSimulation() if config.input_files.amber is not None: prmtop = openmm_app.AmberPrmtopFile( config.input_files.amber.topology) topology = prmtop - if config.input_files.amber.coordinates_filetype in ["inpcrd", + if config.input_files.amber.coordinates_filetype in ["inpcrd", "rst7"]: positions = openmm_app.AmberInpcrdFile( config.input_files.amber.coordinates) @@ -112,25 +112,40 @@ def createGamdSimulation(self, config, platform_name, device_index): raise Exception("Invalid input type: %s. Allowed types are: "\ "'pdb' and 'rst7'/'inpcrd'.") gamdSimulation.system = prmtop.createSystem( - nonbondedMethod=nonbondedMethod, - nonbondedCutoff=config.system.nonbonded_cutoff, + nonbondedMethod=nonbondedMethod, + nonbondedCutoff=config.system.nonbonded_cutoff, constraints=constraints) - + + elif config.input_files.charmm is not None: - psf = openmm_app.CharmmPsfFile( - config.input_files.charmm.topology) - pdb_coords_filename = config.input_files.charmm.coordinates - positions, box_vectors = load_pdb_positions_and_box_vectors( - pdb_coords_filename, need_box) + psf = openmm_app.CharmmPsfFile(config.input_files.charmm.topology) + if config.input_files.charmm.coordinates_filetype == "crd": + positions = openmm_app.CharmmCrdFile( + config.input_files.charmm.coordinates) + elif config.input_files.charmm.coordinates_filetype == "pdb": + positions = openmm_app.PDBFile( + config.input_files.charmm.coordinates) + else: + raise Exception("Invalid input type: %s. Allowed types are: "\ + "'crd' and 'pdb'.") + + # Call a method to parse box vectors + if config.input_files.charmm.is_config_box_vector_defined: + psf.setBox(*config.input_files.charmm.get_box_vectors()) + + # Call a method to parse parameters files to be used params = openmm_app.CharmmParameterSet( *config.input_files.charmm.parameters) + topology = psf gamdSimulation.system = psf.createSystem( params=params, - nonbondedMethod=nonbondedMethod, - nonbondedCutoff=config.system.nonbonded_cutoff, + nonbondedMethod=nonbondedMethod, + nonbondedCutoff=config.system.nonbonded_cutoff, + switchDistance=config.system.switch_distance, + ewaldErrorTolerance = config.system.ewald_error_tolerance, constraints=constraints) - + elif config.input_files.gromacs is not None: gro = openmm_app.GromacsGroFile( config.input_files.gromacs.coordinates) @@ -142,10 +157,10 @@ def createGamdSimulation(self, config, platform_name, device_index): topology = top positions = gro gamdSimulation.system = top.createSystem( - nonbondedMethod=nonbondedMethod, - nonbondedCutoff=config.system.nonbonded_cutoff, + nonbondedMethod=nonbondedMethod, + nonbondedCutoff=config.system.nonbonded_cutoff, constraints=constraints) - + elif config.input_files.forcefield is not None: pdb_coords_filename = config.input_files.forcefield.coordinates positions, box_vectors = load_pdb_positions_and_box_vectors( @@ -157,47 +172,47 @@ def createGamdSimulation(self, config, platform_name, device_index): topology = positions gamdSimulation.system = forcefield.createSystem( topology.topology, - nonbondedMethod=nonbondedMethod, - nonbondedCutoff=config.system.nonbonded_cutoff, + nonbondedMethod=nonbondedMethod, + nonbondedCutoff=config.system.nonbonded_cutoff, constraints=constraints) - + else: raise Exception("No valid input files found. OpenMM simulation "\ "not made.") - + if config.integrator.algorithm == "langevin": boost_type_str = config.integrator.boost_type gamdIntegratorFactory = GamdIntegratorFactory() result = gamdIntegratorFactory.get_integrator( - boost_type_str, gamdSimulation.system, config.temperature, - config.integrator.dt, - config.integrator.number_of_steps.conventional_md_prep, - config.integrator.number_of_steps.conventional_md, - config.integrator.number_of_steps.gamd_equilibration_prep, - config.integrator.number_of_steps.gamd_equilibration, - config.integrator.number_of_steps.total_simulation_length, + boost_type_str, gamdSimulation.system, config.temperature, + config.integrator.dt, + config.integrator.number_of_steps.conventional_md_prep, + config.integrator.number_of_steps.conventional_md, + config.integrator.number_of_steps.gamd_equilibration_prep, + config.integrator.number_of_steps.gamd_equilibration, + config.integrator.number_of_steps.total_simulation_length, config.integrator.number_of_steps.averaging_window_interval, - sigma0p=config.integrator.sigma0.primary, + sigma0p=config.integrator.sigma0.primary, sigma0d=config.integrator.sigma0.secondary) - [gamdSimulation.first_boost_group, - gamdSimulation.second_boost_group, - integrator, gamdSimulation.first_boost_type, + [gamdSimulation.first_boost_group, + gamdSimulation.second_boost_group, + integrator, gamdSimulation.first_boost_type, gamdSimulation.second_boost_type] = result integrator.setRandomNumberSeed(config.integrator.random_seed) integrator.setFriction(config.integrator.friction_coefficient) gamdSimulation.integrator = integrator - + else: - raise Exception("Algorithm not implemented:", + raise Exception("Algorithm not implemented:", config.integrator.algorithm) - + if config.barostat is not None: barostat = openmm.MonteCarloBarostat( - config.barostat.pressure, + config.barostat.pressure, config.temperature, config.barostat.frequency) gamdSimulation.system.addForce(barostat) - + properties = {} user_platform_name = platform_name.lower() # @@ -210,7 +225,7 @@ def createGamdSimulation(self, config, platform_name, device_index): properties['CudaPrecision'] = 'mixed' properties['DeviceIndex'] = device_index gamdSimulation.simulation = openmm_app.Simulation( - topology.topology, gamdSimulation.system, + topology.topology, gamdSimulation.system, gamdSimulation.integrator, platform, properties) gamdSimulation.device_index = device_index gamdSimulation.platform = 'CUDA' @@ -218,35 +233,35 @@ def createGamdSimulation(self, config, platform_name, device_index): platform = openmm.Platform.getPlatformByName('OpenCL') properties['DeviceIndex'] = device_index gamdSimulation.simulation = openmm_app.Simulation( - topology.topology, gamdSimulation.system, + topology.topology, gamdSimulation.system, gamdSimulation.integrator, platform, properties) gamdSimulation.device_index = device_index gamdSimulation.platform = 'OpenCL' else: platform = openmm.Platform.getPlatformByName(platform_name) gamdSimulation.simulation = openmm_app.Simulation( - topology.topology, gamdSimulation.system, + topology.topology, gamdSimulation.system, gamdSimulation.integrator, platform) gamdSimulation.platform = platform_name - + gamdSimulation.simulation.context.setPositions(positions.positions) - if box_vectors is not None: + if box_vectors is not None and config.input_files.charmm is None: gamdSimulation.simulation.context.setPeriodicBoxVectors( *box_vectors) if config.run_minimization: gamdSimulation.simulation.minimizeEnergy() - + gamdSimulation.simulation.context.setVelocitiesToTemperature( config.temperature) - + if config.outputs.reporting.coordinates_file_type == "dcd": gamdSimulation.traj_reporter = openmm_app.DCDReporter - + elif config.outputs.reporting.coordinates_file_type == "pdb": gamdSimulation.traj_reporter = openmm_app.PDBReporter - + else: - raise Exception("Reporter type not found:", + raise Exception("Reporter type not found:", config.outputs.reporting.coordinates_file_type) return gamdSimulation diff --git a/gamd/parser.py b/gamd/parser.py index 7020d82..6f5c459 100644 --- a/gamd/parser.py +++ b/gamd/parser.py @@ -48,6 +48,21 @@ def assign_tag(tag, func, useunit=None): else: return None +def parse_and_assign_charmm_gui_toppar_file(charmm_config, xml_params_filename): + extlist = ['rtf', 'prm', 'str'] + parFiles = () + # this input file lists a series of parameter files to be parsed. This is + # the default approach used by CHARMM-GUI (in their toppar.str file) + for line in open(xml_params_filename.text, 'r'): + if '!' in line: line = line.split('!')[0] + parfile = line.strip() + if len(parfile) != 0: + ext = parfile.lower().split('.')[-1] + if not ext in extlist: continue + # Appending each file listed in inputfile to the existing + # list of parameters files to be read + charmm_config.parameters.append(parfile) + return charmm_config class Parser: def __init__(self): @@ -71,6 +86,13 @@ def parse_system_tag(tag): elif system_tag.tag == "constraints": system_config.constraints \ = assign_tag(system_tag, str).lower() + elif system_tag.tag == "switch-distance": + system_config.switch_distance \ + = assign_tag(system_tag, float, + useunit=unit.nanometer) + elif system_tag.tag == "ewald-error-tolerance": + system_config.ewald_error_tolerance \ + = assign_tag(system_tag, float) else: print("Warning: parameter in XML not found in system "\ "tag. Spelling error?", system_tag.tag) @@ -185,15 +207,40 @@ def parse_charmm_tag(input_files_tag): charmm_config.topology = assign_tag(charmm_tag, str) elif charmm_tag.tag == "coordinates": charmm_config.coordinates = assign_tag(charmm_tag, str) + if "type" in charmm_tag.attrib: + type_attrib = charmm_tag.attrib["type"] + charmm_config.coordinates_filetype = type_attrib + elif charmm_tag.tag == "box-vectors": + box_dict = {} + for box_info in charmm_tag: + charmm_config.config_box_vector_defined = True + if box_info.tag in ["a", "b", "c"]: + box_dict[box_info.tag] = (assign_tag(box_info, float, + useunit=unit.nanometer)) + elif box_info.tag in ["alpha", "beta", "gamma"]: + box_dict[box_info.tag] = (assign_tag(box_info, float, + useunit=unit.degree)) + else: + raise Exception("Unkown parameter '" +box_info.tag+ "'. "\ + "Accepted box-vector parameters are 'a', 'b', 'c', "\ + "'alpha', 'beta' and 'gamma'." ) + + charmm_config.box_vectors = box_dict + elif charmm_tag.tag == "parameters": charmm_config.parameters = [] for xml_params_filename in charmm_tag: - charmm_config.parameters.append( - assign_tag(xml_params_filename, str)) + if "type" in xml_params_filename.attrib: + if xml_params_filename.attrib["type"] == "charmm-gui-toppar": + # parsing list of parameter files like CHARMM-GUI does + charmm_config = parse_and_assign_charmm_gui_toppar_file( + charmm_config, xml_params_filename) + else: + charmm_config.parameters.append( + assign_tag(xml_params_filename, str)) else: print("Warning: parameter in XML not found in "\ - "charmm tag. Spelling error?", - charmm_tag.tag) + "charmm tag. Spelling error?", charmm_tag.tag) return charmm_config def parse_gromacs_tag(input_files_tag): @@ -388,4 +435,4 @@ def parse_file(self, input_file, input_type): if __name__ == "__main__": - pass \ No newline at end of file + pass