diff --git a/proteinbenchmark/gmx_simulation.py b/proteinbenchmark/gmx_simulation.py index 63b73cd..2c78b3a 100644 --- a/proteinbenchmark/gmx_simulation.py +++ b/proteinbenchmark/gmx_simulation.py @@ -5,8 +5,6 @@ import numpy from openff.units import unit -from proteinbenchmark.utilities import exists_and_not_empty, read_xml - class GMXSimulation: """A class representing a simulation in GROMACS.""" @@ -16,17 +14,16 @@ def __init__( gmx_executable: str, parametrized_system: str, initial_coords_file: str, - save_state_prefix: str, - setup_prefix: str, + output_prefix: str, + nonbonded_cutoff: unit.Quantity, + vdw_switch_width: unit.Quantity, temperature: unit.Quantity, pressure: unit.Quantity, + langevin_friction: unit.Quantity, barostat_constant: unit.Quantity, - thermostat_constant: unit.Quantity, timestep: unit.Quantity, traj_length: unit.Quantity, frame_length: unit.Quantity, - restraints_present: bool, - load_state_prefix: str = None, ): """ Initializes the simulation parameters and checks units. @@ -40,16 +37,22 @@ def __init__( The path to the parametrized system as a top file. initial_coords_file Path to GRO file used to set initial coordinates. - save_state_prefix - Path prefix for files to write serialized simulation states. + output_prefix + Path prefix for GMX output files. + nonbonded_cutoff + The cutoff for the Lennard-Jones potential and PME direct space + summation. + vdw_switch_width + The distance from the nonbonded cutoff at which to apply the + switching function. temperature The target temperature of the Langevin thermostat. pressure The target pressure of the Monte Carlo barostat. + langevin_friction + The collision frequency of the Langevin integrator. barostat_constant The time constant for pressure coupling - thermostat_constant - The time constant for temperature coupling timestep The timestep for the Langevin integrator. traj_length @@ -64,17 +67,19 @@ def __init__( self.gmx_executable = gmx_executable self.parametrized_system = parametrized_system self.initial_coords_file = initial_coords_file - self.save_state_prefix = save_state_prefix - self.setup_prefix = setup_prefix - self.load_prefix = load_state_prefix + self.output_prefix = output_prefix # Check units of arguments + if not nonbonded_cutoff.is_compatible_with(unit.nanometer): + raise ValueError("nonbonded_cutoff does not have units of Length") + if not vdw_switch_width.is_compatible_with(unit.nanometer): + raise ValueError("vdw_switch_width does not have units of Length") if not temperature.is_compatible_with(unit.kelvin): raise ValueError("temperature does not have units of Temperature") - if not pressure.is_compatible_with(unit.atmosphere): raise ValueError("pressure does not have units of Mass Length^-1 Time^-2") - + if not langevin_friction.is_compatible_with(unit.picosecond**-1): + raise ValueError("langevin_friction does not have units of Time^-1") if not timestep.is_compatible_with(unit.picosecond): raise ValueError("timestep does not have units of Time") if not traj_length.is_compatible_with(unit.picosecond): @@ -83,59 +88,82 @@ def __init__( raise ValueError("frame_length does not have units of Time") # Ensure proper units for output + self.nonbonded_cutoff = nonbonded_cutoff.m_as(unit.nanometer) + self.switch_distance = (nonbonded_cutoff - vdw_switch_width).m_as( + unit.nanometer + ) self.temperature = temperature.m_as(unit.kelvin) # Kelvin - self.pressure = pressure.m_as(unit.atmosphere) # atm - self.thermostat_constant = thermostat_constant.m_as(unit.picosecond) # ps + self.pressure = pressure.m_as(unit.bar) # atm + self.thermostat_constant = (1.0 / langevin_friction).m_as(unit.picosecond) # ps self.barostat_constant = barostat_constant.m_as(unit.picosecond) # ps self.timestep = timestep.m_as(unit.picosecond) # fs self.n_steps = int(numpy.round(traj_length / timestep)) # steps self.output_frequency = int(numpy.round(frame_length / timestep)) # steps - self.restraints_present = restraints_present # str - def setup_simulation(self): + def setup_simulation(self, continuation: bool = False): """ - Set up a GROMACS simulation. + Set up a GROMACS simulation by writing the MDP file. + + Parameters + ---------- + continuation + Continue from an existing simulation. """ # create mdp file - mdp_file = self.save_state_prefix + ".mdp" + mdp_file = self.output_prefix + ".mdp" # Write settings to MDP File with open(mdp_file, "w") as mdp_file_w: - # Apply position restraints for equilibration - if self.restraints_present == "NPT": - mdp_file_w.write("define=-DPOSRES\n" + "continuation=no\n") + if continuation: + mdp_file_w.write("continuation=yes\ngen_vel = no\n") + else: - mdp_file_w.write("continuation=yes\n") + mdp_file_w.write( + "continuation=no\n" + "gen_vel = yes\n" + f"gen_temp = {self.temperature}\n" + "gen_seed = -1\n" + ) mdp_file_w.write( - f"integrator = md\nnsteps={self.n_steps}\n" - f"dt = {self.timestep}\nnstenergy = 5000\nnstlog = 5000\n" + f"integrator = sd\n" + f"nsteps={self.n_steps}\n" + f"dt = {self.timestep}\n" + f"nstenergy = {self.output_frequency}\n" + f"nstlog = {self.output_frequency}\n" f"nstxout-compressed = {self.output_frequency}\n" - "constraint_algorithm = lincs\nlincs-warnangle = 45\n" - "constraints = h-bonds\nlincs_iter = 2\nlincs_order = 4\n" - "cutoff-scheme = Verlet\nns_type = grid\nnstlist = 40\n" - "rlist = 0.9\nvdwtype = cutoff\nvdw-modifier = force-switch\n" - "rvdw-switch = 0.88\nrvdw = 0.9\ncoulombtype = PME\n" - "rcoulomb = 0.9\npme_order = 4\nfourierspacing = 0.16\n" - "tcoupl = V-rescale\ntc-grps = System\n" + "constraint_algorithm = lincs\n" + "lincs-warnangle = 45\n" + "constraints = h-bonds\n" + "lincs_iter = 2\n" + "lincs_order = 4\n" + "cutoff-scheme = Verlet\n" + "nstlist = 40\n" + "vdwtype = cut-off\n" + "vdw-modifier = potential-switch\n" + f"rvdw = {self.nonbonded_cutoff}\n" + f"rvdw-switch = {self.switch_distance}\n" + "coulombtype = PME\n" + f"rcoulomb = {self.nonbonded_cutoff}\n" + "pme_order = 4\n" + "fourierspacing = 0.16\n" + "tc-grps = System\n" f"tau_t = {self.thermostat_constant}\n" - f"ref_t = {self.temperature}\npbc = xyz\nDispCorr = EnerPres\n" - f"pcoupltype = isotropic\ntau_p = {self.barostat_constant}\n" - f"ref_p = {self.pressure}\ncompressibility = 4.5e-5\n" + f"ref_t = {self.temperature}\n" + "pbc = xyz\n" + "DispCorr = EnerPres\n" + "pcoupl = C-rescale\n" + "pcoupltype = isotropic\n" + f"tau_p = {self.barostat_constant}\n" + f"ref_p = {self.pressure}\n" + "compressibility = 4.5e-5\n" "refcoord_scaling = com\n" ) - if self.restraints_present == "NPT": - mdp_file_w.write( - f"gen_vel = yes\ngen_temp = {self.temperature}\n" - "gen_seed = -1\npcoupl = C-rescale\n" - ) - else: - mdp_file_w.write("gen_vel = no\npcoupl = Parrinello-Rahman\n") return mdp_file - def run(self): + def start_from_pdb(self): """ Start a new simulation initializing positions from a PDB and velocities to random samples from a Boltzmann distribution at the simulation @@ -143,102 +171,104 @@ def run(self): """ # Generate MDP - mdp_file = self.setup_simulation() + mdp_file = self.setup_simulation(continuation=False) # Generate TPR - if self.restraints_present == "NPT": - grompp = subprocess.run( - [ - self.gmx_executable, - "grompp", - "-f", - mdp_file, - "-p", - self.parametrized_system, - "-c", - self.initial_coords_file, - "-r", - self.initial_coords_file, - "-maxwarn", - "2", - "-o", - str(self.save_state_prefix), - ] - ) - - else: - print(self.gmx_executable) - grompp = subprocess.run( - [ - self.gmx_executable, - "grompp", - "-f", - mdp_file, - "-p", - self.parametrized_system, - "-c", - self.initial_coords_file, - "-t", - f"{self.load_prefix}.cpt", - "-maxwarn", - "2", - "-o", - str(self.save_state_prefix), - ] - ) + grompp = subprocess.run( + [ + self.gmx_executable, + "grompp", + "-f", + mdp_file, + "-p", + self.parametrized_system, + "-c", + self.initial_coords_file, + "-maxwarn", + "2", + "-o", + str(self.output_prefix), + ] + ) - # Run Simulation - if self.restraints_present == "NPT": - simulation = subprocess.run( - [ - self.gmx_executable, - "mdrun", - "-deffnm", - str(self.save_state_prefix), - "-ntmpi", - "1", - ] - ) - # simulation = gmx.mdrun( - # input=grompp.output.file['-o'], runtime_args={'-ntmpi': '1'} - # ) - else: - simulation = subprocess.run( - [ - self.gmx_executable, - "mdrun", - "-deffnm", - str(self.save_state_prefix), - "-ntomp", - "1", - ] - ) - # simulation = gmx.mdrun( - # input=grompp.output.file['-o'], runtime_args={'-ntomp': '1'} - # ) + # Run simulation + simulation = subprocess.run( + [ + self.gmx_executable, + "mdrun", + "-deffnm", + str(self.output_prefix), + "-ntmpi", + "1", + ] + ) def start_from_save_state( self, - save_state_file: str, + checkpoint_file: str, ): """ Start a new simulation initializing positions and velocities from a - serialized OpenMM state. + GMX checkpoint. Parameters ---------- - save_state_file - Path to the serialized simulation state. + checkpoint_file + Path to the checkpoint file. """ + + # Generate MDP + mdp_file = self.setup_simulation(continuation=True) + + # Generate TPR + grompp = subprocess.run( + [ + self.gmx_executable, + "grompp", + "-f", + mdp_file, + "-p", + self.parametrized_system, + "-c", + self.initial_coords_file, + "-t", + checkpoint_file, + "-maxwarn", + "2", + "-o", + str(self.output_prefix), + ] + ) + + # Run Simulation + simulation = subprocess.run( + [ + self.gmx_executable, + "mdrun", + "-deffnm", + str(self.output_prefix), + "-ntomp", + "1", + ] + ) + # simulation = gmx.mdrun( + # input=grompp.output.file['-o'], runtime_args={'-ntomp': '1'} + # ) + + def resume_from_checkpoint(self): + """ + Resume an existing simulation from a GMX checkpoint. + """ + # Run Simulation simulation = subprocess.run( [ self.gmx_executable, "mdrun", "-deffnm", - str(self.save_state_prefix), + str(self.output_prefix), "-cpi", - save_state_file, + f"{self.output_prefix}.cpt", "-ntomp", "1", ] diff --git a/proteinbenchmark/openmm_simulation.py b/proteinbenchmark/openmm_simulation.py index 053ed6f..4b0c22a 100644 --- a/proteinbenchmark/openmm_simulation.py +++ b/proteinbenchmark/openmm_simulation.py @@ -78,10 +78,8 @@ def __init__( # Check units of arguments if not temperature.unit.is_compatible(unit.kelvin): raise ValueError("temperature does not have units of Temperature") - if not pressure.unit.is_compatible(unit.atmosphere): raise ValueError("pressure does not have units of Mass Length^-1 Time^-2") - if not langevin_friction.unit.is_compatible(unit.picosecond**-1): raise ValueError("langevin_friction does not have units of Time^-1") if not timestep.unit.is_compatible(unit.picosecond): diff --git a/proteinbenchmark/protein_system.py b/proteinbenchmark/protein_system.py index bb6bf09..1c191f6 100644 --- a/proteinbenchmark/protein_system.py +++ b/proteinbenchmark/protein_system.py @@ -1,7 +1,5 @@ from pathlib import Path -import numpy - from proteinbenchmark.analysis import ( align_trajectory, assign_dihedral_clusters, @@ -87,9 +85,20 @@ def __init__( f'"{quantity}"' ) - self.system_name = f"{target_name}-{force_field_name}" + # Set up nonbonded parameters. Do this here so these values are + # persistent. + if "nonbonded_cutoff" in self.target_parameters: + self.nonbonded_cutoff = self.target_parameters["nonbonded_cutoff"] + else: + self.nonbonded_cutoff = NONBONDED_CUTOFF + + if "vdw_switch_width" in self.target_parameters: + self.vdw_switch_width = self.target_parameters["vdw_switch_width"] + else: + self.vdw_switch_width = VDW_SWITCH_WIDTH # Create a directory to store results for this benchmark system + self.system_name = f"{target_name}-{force_field_name}" self.base_path = Path(result_directory, self.system_name) # File paths for setup @@ -181,16 +190,6 @@ def setup(self): else: solvent_padding = SOLVENT_PADDING - if "nonbonded_cutoff" in self.target_parameters: - nonbonded_cutoff = self.target_parameters["nonbonded_cutoff"] - else: - nonbonded_cutoff = NONBONDED_CUTOFF - - if "vdw_switch_width" in self.target_parameters: - vdw_switch_width = self.target_parameters["vdw_switch_width"] - else: - vdw_switch_width = VDW_SWITCH_WIDTH - solvate( ionic_strength=self.target_parameters["ionic_strength"], protonated_pdb_file=self.protonated_pdb, @@ -203,20 +202,10 @@ def setup(self): if not exists_and_not_empty(self.parametrized_system): print(f"Parametrizing system {self.system_name}") - if "nonbonded_cutoff" in self.target_parameters: - nonbonded_cutoff = self.target_parameters["nonbonded_cutoff"] - else: - nonbonded_cutoff = NONBONDED_CUTOFF - - if "vdw_switch_width" in self.target_parameters: - vdw_switch_width = self.target_parameters["vdw_switch_width"] - else: - vdw_switch_width = VDW_SWITCH_WIDTH - assign_parameters( simulation_platform=self.simulation_platform, - nonbonded_cutoff=nonbonded_cutoff, - vdw_switch_width=vdw_switch_width, + nonbonded_cutoff=self.nonbonded_cutoff, + vdw_switch_width=self.vdw_switch_width, protonated_pdb_file=self.protonated_pdb, solvated_pdb_file=solvated_pdb, parametrized_system=self.parametrized_system, @@ -231,36 +220,38 @@ def setup(self): if not exists_and_not_empty(self.minimized_coords): print(f"Minimizing energy for system {self.system_name}") - if self.simulation_platform == "openmm": - if "restraint_energy_constant" in self.target_parameters: - restraint_energy_constant = self.target_parameters[ - "restraint_energy_constant" - ] + if "restraint_energy_constant" in self.target_parameters: + restraint_energy_constant = self.target_parameters[ + "restraint_energy_constant" + ] + else: + restraint_energy_constant = RESTRAINT_ENERGY_CONSTANT - else: - restraint_energy_constant = RESTRAINT_ENERGY_CONSTANT + if "force_tolerance" in self.target_parameters: + force_tolerance = self.target_parameters["force_tolerance"] + else: + force_tolerance = FORCE_TOLERANCE + if self.simulation_platform == "openmm": minimize_openmm( parametrized_system=self.parametrized_system, solvated_pdb_file=solvated_pdb, minimized_coords_file=self.minimized_coords, restraint_energy_constant=restraint_energy_constant, + force_tolerance=force_tolerance, ) elif self.simulation_platform == "gmx": - if "energy_tolerance" in self.target_parameters: - energy_tolerance = self.target_parameters["energy_tolerance"] - - else: - energy_tolerance = ENERGY_TOLERANCE - minimize_gmx( parametrized_system=self.parametrized_system, solvated_pdb_file=solvated_pdb, minimized_coords_file=self.minimized_coords, setup_prefix=self.setup_prefix, gmx_executable=self.gmx_executable, - energy_tolerance=energy_tolerance, + nonbonded_cutoff=self.nonbonded_cutoff, + vdw_switch_width=self.vdw_switch_width, + restraint_energy_constant=restraint_energy_constant, + force_tolerance=force_tolerance, ) print(f"Setup complete for system {self.system_name}") @@ -283,7 +274,7 @@ def run_simulations(self, replica: int = 1): if self.simulation_platform == "openmm": equilibrated_state = f"{equil_prefix}-1.xml" elif self.simulation_platform == "gmx": - equilibrated_state = f"{equil_prefix}.gro" + equilibrated_state = f"{equil_prefix}.cpt" # Equilibrate at constant pressure and temperature if not exists_and_not_empty(equilibrated_state): @@ -305,25 +296,23 @@ def run_simulations(self, replica: int = 1): else: equil_frame_length = EQUIL_FRAME_LENGTH + if "equil_langevin_friction" in self.target_parameters: + equil_langevin_friction = self.target_parameters[ + "equil_langevin_friction" + ] + else: + equil_langevin_friction = EQUIL_LANGEVIN_FRICTION + + # Initialize the equilibration simulation if self.simulation_platform == "openmm": # Get OpenMM Specific Parameters - if "equil_langevin_friction" in self.target_parameters: - equil_langevin_friction = self.target_parameters[ - "equil_langevin_friction" - ] - - else: - equil_langevin_friction = EQUIL_OPENMM_LANGEVIN_FRICTION - if "equil_barostat_frequency" in self.target_parameters: equil_barostat_frequency = self.target_parameters[ "equil_barostat_frequency" ] - else: equil_barostat_frequency = EQUIL_OPENMM_BAROSTAT_FREQUENCY - # Initialize the equilibration simulation equilibration_dcd = f"{equil_prefix}.dcd" equilibration_state_data = f"{equil_prefix}.out" equilibration_checkpoint = f"{equil_prefix}.chk" @@ -346,9 +335,6 @@ def run_simulations(self, replica: int = 1): save_state_length=equil_traj_length.to_openmm(), ) - # Run equilibration - equilibration_simulation.start_from_pdb() - else: # Get GROMACS Specific Parameters if "equil_barostat_constant" in self.target_parameters: @@ -358,30 +344,24 @@ def run_simulations(self, replica: int = 1): else: equil_barostat_constant = EQUIL_GMX_BAROSTAT_CONSTANT - if "equil_thermostat_constant" in self.target_parameters: - equil_thermostat_constant = self.target_parameters[ - "equil_thermostat_constant" - ] - else: - equil_thermostat_constant = EQUIL_GMX_THERMOSTAT_CONSTANT - - NPT_simulation = GMXSimulation( + equilibration_simulation = GMXSimulation( gmx_executable=self.gmx_executable, - initial_coords_file=self.minimized_coords, parametrized_system=self.parametrized_system, - save_state_prefix=equil_prefix, - setup_prefix=setup_prefix, + initial_coords_file=self.minimized_coords, + output_prefix=equil_prefix, + nonbonded_cutoff=self.nonbonded_cutoff, + vdw_switch_width=self.vdw_switch_width, temperature=self.target_parameters["temperature"], pressure=self.target_parameters["pressure"], + langevin_friction=equil_langevin_friction, barostat_constant=equil_barostat_constant, - thermostat_constant=equil_thermostat_constant, timestep=equil_timestep, traj_length=equil_traj_length, frame_length=equil_frame_length, - restraints_present="NPT", ) - NPT_simulation.run() + # Run equilibration + equilibration_simulation.start_from_pdb() print(f"Running NPT production for system {self.system_name}") @@ -411,6 +391,12 @@ def run_simulations(self, replica: int = 1): else: frame_length = FRAME_LENGTH + if "langevin_friction" in self.target_parameters: + langevin_friction = self.target_parameters["langevin_friction"] + else: + langevin_friction = LANGEVIN_FRICTION + + # Initialize the production simulation if self.simulation_platform == "openmm": # Get OpenMM Specific Parameters if "checkpoint_length" in self.target_parameters: @@ -423,17 +409,11 @@ def run_simulations(self, replica: int = 1): else: save_state_length = SAVE_STATE_LENGTH - if "langevin_friction" in self.target_parameters: - langevin_friction = self.target_parameters["langevin_friction"] - else: - langevin_friction = OPENMM_LANGEVIN_FRICTION - if "barostat_frequency" in self.target_parameters: barostat_frequency = self.target_parameters["barostat_frequency"] else: barostat_frequency = OPENMM_BAROSTAT_FREQUENCY - # Initialize the production simulation production_dcd = f"{prod_prefix}.dcd" production_state_data = f"{prod_prefix}.out" production_checkpoint = f"{prod_prefix}.chk" @@ -456,16 +436,6 @@ def run_simulations(self, replica: int = 1): save_state_length=save_state_length.to_openmm(), ) - # Run production - if not exists_and_not_empty(production_checkpoint): - # Start production simulation, initializing positions and - # velocities to the final state from the equilibration simulation - production_simulation.start_from_save_state(equilibrated_state) - - else: - # Resume from a previous production checkpoint - production_simulation.resume_from_checkpoint() - else: # Get GROMACS Specific Parameters if "barostat_constant" in self.target_parameters: @@ -473,37 +443,33 @@ def run_simulations(self, replica: int = 1): else: barostat_constant = GMX_BAROSTAT_CONSTANT - if "thermostat_constant" in self.target_parameters: - thermostat_constant = self.target_parameters["thermostat_constant"] - else: - thermostat_constant = GMX_THERMOSTAT_CONSTANT - production_checkpoint = f"{prod_prefix}.cpt" production_simulation = GMXSimulation( gmx_executable=self.gmx_executable, parametrized_system=self.parametrized_system, initial_coords_file=equil_prefix, - setup_prefix=setup_prefix, - save_state_prefix=prod_prefix, + output_prefix=prod_prefix, + nonbonded_cutoff=self.nonbonded_cutoff, + vdw_switch_width=self.vdw_switch_width, temperature=self.target_parameters["temperature"], pressure=self.target_parameters["pressure"], + langevin_friction=langevin_friction, barostat_constant=barostat_constant, - thermostat_constant=thermostat_constant, timestep=timestep, traj_length=traj_length, frame_length=frame_length, - restraints_present=False, - load_state_prefix=equil_prefix, ) - # Run Production - if not exists_and_not_empty(production_checkpoint): - # Start production simulation, initializing positions and - # velocities to the final state from the equilibration simulation - production_simulation.run() - else: - production_simulation.start_from_save_state(production_checkpoint) + # Run production + if not exists_and_not_empty(production_checkpoint): + # Start production simulation, initializing positions and + # velocities to the final state from the equilibration simulation + production_simulation.start_from_save_state(equilibrated_state) + + else: + # Resume from a previous production checkpoint + production_simulation.resume_from_checkpoint() def analyze_observables(self, replica: int = 1): """Process trajectories and estimate observables.""" diff --git a/proteinbenchmark/simulation_parameters.py b/proteinbenchmark/simulation_parameters.py index 07751a9..0f845c4 100644 --- a/proteinbenchmark/simulation_parameters.py +++ b/proteinbenchmark/simulation_parameters.py @@ -11,26 +11,24 @@ # Default minimization parameters RESTRAINT_ENERGY_CONSTANT = 1.0 * unit.kilocalorie_per_mole / unit.angstrom**2 -ENERGY_TOLERANCE = 100 * unit.kilojoules_per_mole / unit.nanometer +FORCE_TOLERANCE = 10 * unit.kilojoules_per_mole / unit.nanometer # Default equilibration simulation parameters EQUIL_TIMESTEP = 1.0 * unit.femtosecond EQUIL_TRAJ_LENGTH = 1.0 * unit.nanosecond EQUIL_FRAME_LENGTH = 10.0 * unit.picosecond -EQUIL_OPENMM_LANGEVIN_FRICTION = 5.0 / unit.picosecond +EQUIL_LANGEVIN_FRICTION = 5.0 / unit.picosecond EQUIL_OPENMM_BAROSTAT_FREQUENCY = 5 EQUIL_GMX_BAROSTAT_CONSTANT = 10.0 * unit.picosecond -EQUIL_GMX_THERMOSTAT_CONSTANT = 2.0 * unit.picosecond # Default production simulation parameters TIMESTEP = 4.0 * unit.femtosecond FRAME_LENGTH = 100.0 * unit.picosecond CHECKPOINT_LENGTH = FRAME_LENGTH * 100 SAVE_STATE_LENGTH = CHECKPOINT_LENGTH * 10 -OPENMM_LANGEVIN_FRICTION = 1.0 / unit.picosecond +LANGEVIN_FRICTION = 1.0 / unit.picosecond OPENMM_BAROSTAT_FREQUENCY = 25 GMX_BAROSTAT_CONSTANT = 10.0 * unit.picosecond -GMX_THERMOSTAT_CONSTANT = 2.0 * unit.picosecond # Production trajectory length for peptides, folded proteins, and disordered # proteins diff --git a/proteinbenchmark/system_setup.py b/proteinbenchmark/system_setup.py index 5b5fb78..334d385 100644 --- a/proteinbenchmark/system_setup.py +++ b/proteinbenchmark/system_setup.py @@ -627,7 +627,7 @@ def assign_parameters( # Manually change hydrogen masses in OpenMM system. Taken from # https://github.com/openmm/openmm/blob/f30d716ace8331003c5115bdfa9e03341a757878/wrappers/python/openmm/app/forcefield.py#L1249 _hydrogen = app.element.hydrogen - for atom1, atom2 in modeller.topology.bonds(): + for atom1, atom2 in openmm_topology.bonds(): if atom1.element == _hydrogen: (atom1, atom2) = (atom2, atom1) if ( @@ -767,11 +767,12 @@ def assign_parameters( write_xml(parametrized_system, openmm_system) elif simulation_platform == "gmx": + # TODO write GMX files with Interchange import parmed as pmd # Write GROMACS files struct = pmd.openmm.load_topology( - openmm_topology, openmm_system, xyz=modeller.positions + openmm_topology, openmm_system, xyz=solvated_pdb.positions ) hmass = pmd.tools.HMassRepartition(struct, hydrogen_mass.m_as(unit.dalton)) hmass.execute() @@ -809,6 +810,7 @@ def minimize_openmm( solvated_pdb_file: str, minimized_coords_file: str, restraint_energy_constant: unit.Quantity, + force_tolerance: unit.Quantity, ): """ Minimize energy of solvated system with Cartesian restraints on non-hydrogen @@ -823,8 +825,9 @@ def minimize_openmm( minimized_coords_file The path to write the minimized coords (OpenMM PDB or GMX GRO). restraint_energy_constant - Energy constant for Cartesian restraints in OpenMM (units Energy - Length^-2). + Energy constant for Cartesian restraints (units Energy Length^-2). + force_tolerance + Force tolerance for minimization (units Energy Length^-1). """ # Check units of arguments @@ -834,6 +837,8 @@ def minimize_openmm( raise ValueError( "restraint_energy_constant does not have units of Energy Length^-2" ) + if not force_tolerance.is_compatible_with(unit.kilojoule_per_mole / unit.nanometer): + raise ValueError("force_tolerance does not have units of Energy Length^-1") k = restraint_energy_constant.m_as(unit.kilocalories_per_mole / unit.angstrom**2) @@ -906,7 +911,10 @@ def minimize_gmx( minimized_coords_file: str, setup_prefix: str, gmx_executable: str, - energy_tolerance: unit.Quantity, + nonbonded_cutoff: unit.Quantity, + vdw_switch_width: unit.Quantity, + restraint_energy_constant: unit.Quantity, + force_tolerance: unit.Quantity, ): """ Minimize energy of solvated system with Cartesian restraints on non-hydrogen @@ -922,30 +930,71 @@ def minimize_gmx( The path to write the minimized coords (OpenMM PDB or GMX GRO). gmx_executable Name of GROMACS executable to pass to subprocess. - energy_tolerance - Energy tolerance for minimization in GROMACS (units Energy Length^-1). + nonbonded_cutoff + The cutoff for the Lennard-Jones potential and PME direct space + summation. + vdw_switch_width + The distance from the nonbonded cutoff at which to apply the + switching function. + restraint_energy_constant + Energy constant for Cartesian restraints (units Energy Length^-2). + force_tolerance + Force tolerance for minimization (units Energy Length^-1). """ import subprocess # Check units of arguments - if not energy_tolerance.is_compatible_with( - unit.kilojoule_per_mole / unit.nanometer + if not nonbonded_cutoff.is_compatible_with(unit.nanometer): + raise ValueError("nonbonded_cutoff does not have units of Length") + if not vdw_switch_width.is_compatible_with(unit.nanometer): + raise ValueError("vdw_switch_width does not have units of Length") + if not restraint_energy_constant.is_compatible_with( + unit.kilojoule_per_mole / unit.nanometer**2 ): - raise ValueError("energy_tolerance does not have units of Energy Length^-1") + raise ValueError( + "restraint_energy_constant does not have units of Energy Length^-2" + ) + if not force_tolerance.is_compatible_with(unit.kilojoule_per_mole / unit.nanometer): + raise ValueError("force_tolerance does not have units of Energy Length^-1") - k = energy_tolerance.m_as(unit.kilojoule_per_mole / unit.nanometer) + nonbonded_cutoff = nonbonded_cutoff.m_as(unit.nanometer) + vdw_switch_width = vdw_switch_width.m_as(unit.nanometer) + switch_distance = nonbonded_cutoff - vdw_switch_width + restraint_energy_constant = restraint_energy_constant.m_as( + unit.kilojoule_per_mole / unit.nanometer**2 + ) + force_tolerance = force_tolerance.m_as(unit.kilojoule_per_mole / unit.nanometer) # Create MDP file mdp_file = str(setup_prefix) + "-minimized.mdp" with open(mdp_file, "w") as mdp_file_w: mdp_file_w.write( - f"integrator = steep\nemtol = {k}\nemstep = 0.01\nnsteps=50000\n" - "nstlist = 1\ncutoff-scheme = Verlet\nns_type = grid\nrlist = 0.9\n" - "coulombtype = PME\nrcoulomb = 1.0\nrvdw = 0.9\npbc = xyz\n" + "define=-DPOSRES\n" + "integrator = steep\n" + f"emtol = {10*force_tolerance}\n" + "emstep = 0.01\n" + "nsteps = -1\n" + "constraint_algorithm = lincs\n" + "lincs-warnangle = 45\n" + "constraints = h-bonds\n" + "lincs_iter = 2\n" + "lincs_order = 4\n" + "cutoff-scheme = Verlet\n" + "nstlist = 40\n" + "vdwtype = cut-off\n" + "vdw-modifier = potential-switch\n" + f"rvdw = {nonbonded_cutoff}\n" + f"rvdw-switch = {switch_distance}\n" + "coulombtype = PME\n" + f"rcoulomb = {nonbonded_cutoff}\n" + "pme_order = 4\n" + "fourierspacing = 0.16\n" + "pbc = xyz\n" + "DispCorr = EnerPres\n" ) - # Create position restrints file for backbone atoms + # Create position restraints file for backbone atoms restr = subprocess.Popen( [ gmx_executable, @@ -954,6 +1003,8 @@ def minimize_gmx( solvated_pdb_file, "-o", f"{setup_prefix}_posre.itp", + "-fc", + str(restraint_energy_constant), ], stdin=subprocess.PIPE, ) @@ -972,6 +1023,10 @@ def minimize_gmx( parametrized_system, "-c", solvated_pdb_file, + "-r", + solvated_pdb_file, + "-maxwarn", + f"{1}", "-o", out_prefix, ]