diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index c2c65ebb..6d047c0a 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -16,6 +16,7 @@ dependencies: - pyyaml - plumed - intermol + - openff-interchange >=0.3.7 - openff-units >=0.2.0 - openff-utilities diff --git a/docs/tutorials/08-generating-taproom-systems.ipynb b/docs/tutorials/08-generating-taproom-systems.ipynb new file mode 100644 index 00000000..f4ec791e --- /dev/null +++ b/docs/tutorials/08-generating-taproom-systems.ipynb @@ -0,0 +1,891 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4fa4e04a-c418-46d6-b4c0-92c4f20a484a", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "LICENSE: Could not open license file \"oe_license.txt\" in local directory\n", + "LICENSE: N.B. OE_LICENSE environment variable is not set\n", + "LICENSE: N.B. OE_DIR environment variable is not set\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "The OpenEye Toolkits are found to be installed but not licensed and therefore will not be used.\n", + "The OpenEye Toolkits require a (free for academics) license, see https://docs.eyesopen.com/toolkits/python/quickstart-python/license.html\n", + "LICENSE: No product keys!\n", + "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", + "Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "61ce7d2e2022466cb18b85a279cb6b52", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import os\n", + "import warnings\n", + "\n", + "warnings.filterwarnings(action=\"ignore\")\n", + "\n", + "from openff.toolkit.typing.engines.smirnoff import ForceField\n", + "from pkg_resources import resource_filename\n", + "from tqdm.auto import tqdm\n", + "\n", + "from paprika.build.system.taproom import BuildTaproomAPR" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "e9e0c9bd-1771-4c52-9c46-b5ece3d8c7c2", + "metadata": {}, + "outputs": [], + "source": [ + "os.environ[\"INTERCHANGE_EXPERIMENTAL\"] = \"1\"" + ] + }, + { + "cell_type": "markdown", + "id": "e5140dbf", + "metadata": {}, + "source": [ + "* The `BuildTaproomAPR` class requires the OpenFF-Interchange modules. Install in your conda environment with:\n", + "\n", + "`conda install -c conda-forge openff-interchange`" + ] + }, + { + "cell_type": "markdown", + "id": "ac82a11a", + "metadata": {}, + "source": [ + "## 01) Build System in explicit solvent" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b8c3ab15", + "metadata": {}, + "outputs": [], + "source": [ + "# Define force field (OpenFF)\n", + "force_field = ForceField(\"openff-2.0.0.offxml\")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "50986520", + "metadata": {}, + "outputs": [], + "source": [ + "system = BuildTaproomAPR(\n", + " host_code=\"bcd\",\n", + " guest_code=\"hex\",\n", + " n_water=3000,\n", + " force_field=force_field,\n", + " working_folder=\"explicit_solvent\",\n", + " disable_progress=False,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "2c92a6bc", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "a34b780a4b4243da857665d4bcd0ecec", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/46 [00:00:228: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject\n", + "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", + ":228: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject\n", + "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", + "Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.\n", + ":228: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject\n", + "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", + "Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.\n", + ":228: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject\n", + "LICENSE: Could not open license file \"oe_license.txt\" in local directory\n", + "LICENSE: N.B. OE_LICENSE environment variable is not set\n", + "LICENSE: N.B. OE_DIR environment variable is not set\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "The OpenEye Toolkits are found to be installed but not licensed and therefore will not be used.\n", + "The OpenEye Toolkits require a (free for academics) license, see https://docs.eyesopen.com/toolkits/python/quickstart-python/license.html\n", + "Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.\n", + "LICENSE: Could not open license file \"oe_license.txt\" in local directory\n", + "LICENSE: N.B. OE_LICENSE environment variable is not set\n", + "LICENSE: N.B. OE_DIR environment variable is not set\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.\n", + "LICENSE: No product keys!\n", + "The OpenEye Toolkits are found to be installed but not licensed and therefore will not be used.\n", + "The OpenEye Toolkits require a (free for academics) license, see https://docs.eyesopen.com/toolkits/python/quickstart-python/license.html\n", + "LICENSE: Could not open license file \"oe_license.txt\" in local directory\n", + "LICENSE: N.B. OE_LICENSE environment variable is not set\n", + "LICENSE: N.B. OE_DIR environment variable is not set\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "The OpenEye Toolkits are found to be installed but not licensed and therefore will not be used.\n", + "The OpenEye Toolkits require a (free for academics) license, see https://docs.eyesopen.com/toolkits/python/quickstart-python/license.html\n", + "LICENSE: Could not open license file \"oe_license.txt\" in local directory\n", + "LICENSE: N.B. OE_LICENSE environment variable is not set\n", + "LICENSE: N.B. OE_DIR environment variable is not set\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "LICENSE: No product keys!\n", + "The OpenEye Toolkits are found to be installed but not licensed and therefore will not be used.\n", + "The OpenEye Toolkits require a (free for academics) license, see https://docs.eyesopen.com/toolkits/python/quickstart-python/license.html\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/_experimental.py:35: UserWarning: Interchange object combination is experimental and likely to produce strange results. Any workflow using this method is not guaranteed to be suitable for production. Use with extreme caution and thoroughly validate results!\n", + " return func(*args, **kwargs)\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/_experimental.py:35: UserWarning: Interchange object combination is experimental and likely to produce strange results. Any workflow using this method is not guaranteed to be suitable for production. Use with extreme caution and thoroughly validate results!\n", + " return func(*args, **kwargs)\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/_experimental.py:35: UserWarning: Interchange object combination is experimental and likely to produce strange results. Any workflow using this method is not guaranteed to be suitable for production. Use with extreme caution and thoroughly validate results!\n", + " return func(*args, **kwargs)\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/_experimental.py:35: UserWarning: Interchange object combination is experimental and likely to produce strange results. Any workflow using this method is not guaranteed to be suitable for production. Use with extreme caution and thoroughly validate results!\n", + " return func(*args, **kwargs)\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/components/interchange.py:844: UserWarning: Setting positions to None because one or both objects added together were missing positions.\n", + " warnings.warn(\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/components/interchange.py:844: UserWarning: Setting positions to None because one or both objects added together were missing positions.\n", + " warnings.warn(\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/components/interchange.py:844: UserWarning: Setting positions to None because one or both objects added together were missing positions.\n", + " warnings.warn(\n", + "/Users/jsetiadi/opt/anaconda3/envs/paprika-dev/lib/python3.9/site-packages/openff/interchange/components/interchange.py:844: UserWarning: Setting positions to None because one or both objects added together were missing positions.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "e60193dde26d4eee918037f8300d7367", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/15 [00:00= 4 @@ -850,11 +839,15 @@ def run_mbar(self, phase, prepared_data, method, verbose=False): # Estimate uncertainty from decorrelated samples # Create subsampled indices and count their lengths - frac_N_ss = np.array([int(fraction * n) for n in N_ss], dtype=np.int32) + frac_N_ss = numpy.array( + [int(fraction * n) for n in N_ss], dtype=numpy.int32 + ) # Create a new potential array for the uncertainty calculation # (are we using too much memory?) - u_kln_err = np.zeros([num_win, num_win, np.max(frac_N_ss)], np.float64) + u_kln_err = numpy.zeros( + [num_win, num_win, numpy.max(frac_N_ss)], numpy.float64 + ) # Populate the subsampled array, drawing the appropriate # fraction of subsamples from the original @@ -940,7 +933,7 @@ def run_ti(self, phase, prepared_data, method): ---------- phase: str The phase of the calculation to analyze. - prepared_data: :class:`np.array` + prepared_data: :class:`numpy.array` The list of "prepared data" including the number of windows, data points, which restraints are changing, their force constants and targets, and well as the order of the windows. This probably ought to be redesigned. @@ -968,39 +961,44 @@ def run_ti(self, phase, prepared_data, method): ) = prepared_data # Number of data points in each restraint value array - N_k = np.array(data_points) + N_k = numpy.array(data_points) # The dU array to store the partial derivative of the potential with respect lambda or target, # depending on the whether attach/release or pull. Data stored for each frame. This just a # temporary storage space. - dU = np.zeros([num_win, max_data_points], np.float64) + dU = numpy.zeros([num_win, max_data_points], numpy.float64) # The mean, SEM, standard deviation, and number of uncorrelated dU values for each window. - dU_avgs = np.zeros([num_win], np.float64) - dU_sems = np.zeros([num_win], np.float64) - dU_stdv = np.zeros([num_win], np.float64) - dU_Nunc = np.zeros([num_win], np.float64) + dU_avgs = numpy.zeros([num_win], numpy.float64) + dU_sems = numpy.zeros([num_win], numpy.float64) + dU_stdv = numpy.zeros([num_win], numpy.float64) + dU_Nunc = numpy.zeros([num_win], numpy.float64) # The statistical inefficiency - g = np.zeros([num_win], np.float64) + g = numpy.zeros([num_win], numpy.float64) # Array for values of the changing coordinate (x-axis), either lambda or target. # I'll name them dl_vals for dlambda values. - dl_vals = np.zeros([num_win], np.float64) + dl_vals = numpy.zeros([num_win], numpy.float64) # Setup interpolation array for the dLambda (dl) coordinate. We're gonna create # this progressively by appending ... - dl_intp = np.zeros([0], np.float64) + dl_intp = numpy.zeros([0], numpy.float64) # Get units - target_units = np.array([targets[r][0].units for r in range(len(active_rest))]) - force_units = np.array( + target_units = numpy.array( + [targets[r][0].units for r in range(len(active_rest))] + ) + force_units = numpy.array( [force_constants[r][0].units for r in range(len(active_rest))] ) # Store the max force constant value for each restraint. max_force_constants = ( - np.array( - [np.max(force_constants[r]).magnitude for r in range(len(active_rest))] + numpy.array( + [ + numpy.max(force_constants[r]).magnitude + for r in range(len(active_rest)) + ] ) * force_units ) @@ -1008,8 +1006,10 @@ def run_ti(self, phase, prepared_data, method): # Transpose force_constants and targets into "per window" format, instead of # the "per restraint" format. # print(targets) - force_constants_T = np.asarray(force_constants).T * force_units - targets_T = np.asarray(targets).T * target_units + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + force_constants_T = numpy.asarray(force_constants).T * force_units + targets_T = numpy.asarray(targets).T * target_units # For each window: do dihedral wrapping, compute forces, append dl_intp for k in range(num_win): # Coordinate windows @@ -1070,28 +1070,28 @@ def run_ti(self, phase, prepared_data, method): # Compute standard deviations and SEMs, unless we're going to do # exact_sem_each_ti_fraction - dU_avgs[k] = np.mean(dU[k, 0 : N_k[k]]) - dU_stdv[k] = np.std(dU[k, 0 : N_k[k]]) + dU_avgs[k] = numpy.mean(dU[k, 0 : N_k[k]]) + dU_stdv[k] = numpy.std(dU[k, 0 : N_k[k]]) if method == "ti-block": nearest_max = get_nearest_max(N_k[k]) dU_sems[k] = get_block_sem(dU[k, 0:nearest_max]) # Rearrange SEM = StdDev/sqrt(N) to get N_uncorrelated dU_Nunc[k] = (dU_stdv[k] / dU_sems[k]) ** 2 elif method == "ti-nocor": - dU_sems[k] = dU_stdv[k] / np.sqrt(N_k[k]) + dU_sems[k] = dU_stdv[k] / numpy.sqrt(N_k[k]) dU_Nunc[k] = N_k[k] g[k] = N_k[k] / dU_Nunc[k] # Create the interpolation by appending 100 points between each window. - # Start with k=1 so we don't double count. + # Start with k=1, so we don't double count. if k > 0: - dl_intp = np.append( + dl_intp = numpy.append( dl_intp, - np.linspace(dl_vals[k - 1], dl_vals[k], num=100, endpoint=False), + numpy.linspace(dl_vals[k - 1], dl_vals[k], num=100, endpoint=False), ) # Tack on the final value to the dl interpolation - dl_intp = np.append(dl_intp, dl_vals[-1]) + dl_intp = numpy.append(dl_intp, dl_vals[-1]) logger.debug("Running bootstrap calculations...") @@ -1104,28 +1104,28 @@ def run_ti(self, phase, prepared_data, method): logger.debug("Working on fraction ... {}".format(fraction)) # Compute means for this fraction. - frac_dU_avgs = np.array( - [np.mean(dU[k, 0 : int(fraction * n)]) for k, n in enumerate(N_k)] + frac_dU_avgs = numpy.array( + [numpy.mean(dU[k, 0 : int(fraction * n)]) for k, n in enumerate(N_k)] ) # If self.exact_sem_each_ti_fraction, we're gonna recompute the SEM for each fraction # rather than estimating it from the standard deviation (dU_stdv) and number of # uncorrelated data points (dU_Nunc) from the total data set. if method == "ti-block" and self.exact_sem_each_ti_fraction: - frac_dU_sems = np.zero([k], np.float64) + frac_dU_sems = numpy.zeros([num_win], numpy.float64) for k in range(num_win): nearest_max = get_nearest_max(int(fraction * N_k[k])) frac_dU_sems[k] = get_block_sem(dU[k, 0:nearest_max]) elif method == "ti-nocor" and self.exact_sem_each_ti_fraction: - frac_dU_sems = np.zero([k], np.float64) + frac_dU_sems = numpy.zeros([num_win], numpy.float64) for k in range(num_win): - frac_dU_sems[k] = np.std( + frac_dU_sems[k] = numpy.std( dU[k, 0 : int(fraction * N_k[k])] - ) / np.sqrt(int(fraction * N_k[k])) + ) / numpy.sqrt(int(fraction * N_k[k])) else: - frac_dU_sems = dU_stdv / np.sqrt(fraction * dU_Nunc) + frac_dU_sems = dU_stdv / numpy.sqrt(fraction * dU_Nunc) - dU_samples = np.random.normal( + dU_samples = numpy.random.normal( frac_dU_avgs, frac_dU_sems, size=(self.boot_cycles, frac_dU_avgs.size) ) @@ -1153,12 +1153,12 @@ def run_ti(self, phase, prepared_data, method): if self.compute_roi: logger.info(phase + ": computing ROI for " + method) # Do ROI calc - max_fraction = np.max(self.fractions) + max_fraction = numpy.max(self.fractions) # If we didn't compute fe/sem for fraction 1.0 already, do it now - dU_samples = np.random.normal( + dU_samples = numpy.random.normal( dU_avgs, dU_sems, size=(self.boot_cycles, dU_avgs.size) ) - if not np.isclose(max_fraction, 1.0): + if not numpy.isclose(max_fraction, 1.0): junk_fe, total_sem_matrix = integrate_bootstraps( dl_vals, dU_samples, x_intp=dl_intp, matrix=self.ti_matrix ) @@ -1166,12 +1166,12 @@ def run_ti(self, phase, prepared_data, method): total_sem_matrix = self.results[phase][method]["fraction_sem_matrix"][ max_fraction ].magnitude - self.results[phase][method]["roi"] = np.zeros([num_win], np.float64) + self.results[phase][method]["roi"] = numpy.zeros([num_win], numpy.float64) for k in range(num_win): # Compute overall integrated SEM with 10% smaller SEM for dU[k] - cnvg_dU_samples = np.array(dU_samples) - cnvg_dU_samples[:, k] = np.random.normal( + cnvg_dU_samples = numpy.array(dU_samples) + cnvg_dU_samples[:, k] = numpy.random.normal( dU_avgs[k], 0.9 * dU_sems[k], self.boot_cycles ) junk_fe, cnvg_sem_matrix = integrate_bootstraps( @@ -1229,7 +1229,7 @@ def compute_free_energy(self, phases=["attach", "pull", "release"], seed=None): for method in self.methods: if seed is not None: - np.random.seed(seed) + numpy.random.seed(seed) logger.debug(f"Setting random number seed = {seed}") self.results[phase][method] = {} @@ -1239,7 +1239,7 @@ def compute_free_energy(self, phases=["attach", "pull", "release"], seed=None): logger.debug("Skipping free energy calculation for %s" % phase) continue prepared_data = self.prepare_data(phase) - self.results[phase][method]["n_frames"] = np.sum(prepared_data[1]) + self.results[phase][method]["n_frames"] = numpy.sum(prepared_data[1]) logger.debug( "Running {} analysis on {} phase ...".format(method, phase) @@ -1251,7 +1251,7 @@ def compute_free_energy(self, phases=["attach", "pull", "release"], seed=None): or method == "mbar-boot" ): self.run_mbar(phase, prepared_data, method) - elif method == "ti-block": + elif method == "ti-block" or method == "ti-nocor": self.run_ti(phase, prepared_data, method) else: raise NotImplementedError( @@ -1280,7 +1280,7 @@ def compute_free_energy(self, phases=["attach", "pull", "release"], seed=None): # Set these higher level (total) values, which will be slightly # easier to access - max_fraction = np.max(self.fractions) + max_fraction = numpy.max(self.fractions) self.results[phase][method]["fe_matrix"] = self.results[phase][method][ "fraction_fe_matrix" ][max_fraction] @@ -1301,7 +1301,7 @@ def compute_free_energy(self, phases=["attach", "pull", "release"], seed=None): self.results[phase][method][ "largest_neighbor" ] = openff_unit.Quantity( - np.ones([windows], np.float64) * -1.0, + numpy.ones([windows], numpy.float64) * -1.0, units=self.energy_unit, ) logger.info(f"{phase}: computing largest_neighbor for {method}...") @@ -1405,12 +1405,12 @@ def compute_ref_state_work(self, restraints, state="final"): force_index = 0 fcs.append( - np.sort(restraint.phase[phase]["force_constants"])[ + numpy.sort(restraint.phase[phase]["force_constants"])[ force_index ] ) targs.append( - np.sort(restraint.phase[phase]["targets"])[target_index] + numpy.sort(restraint.phase[phase]["targets"])[target_index] ) target_and_force_exist = True @@ -1473,7 +1473,9 @@ def load_results(self, filepath): self.results = json.loads(data, cls=PaprikaDecoder) - def save_data(self, filepath="simulation_data.json", overwrite=False): + def save_simulation_data_to_json( + self, filepath="simulation_data.json", overwrite=False + ): """ Save the simulation data (DAT values) to a JSON file. @@ -1498,23 +1500,22 @@ def save_data(self, filepath="simulation_data.json", overwrite=False): ) f.write(dumped) - def load_data(self, filepath): + def load_simulation_data_from_json(self, filepath): """ - Load the simulation data (DAT values) from a JSON file. + Read in simulation data from a JSON file. Parameters ---------- - filepath: str - The name of the JSON file to read. + filepath: os.PathLike + The name of the JSON file. """ with open(filepath, "r") as f: - data = f.read() - - simulation_data = json.loads(data, cls=PaprikaDecoder) + json_data = f.read() + data = json.loads(json_data, cls=PaprikaDecoder) - self.simulation_data = simulation_data["simulation_data"] - self.changing_restraints = simulation_data["changing_restraints"] - self.orders = simulation_data["orders"] + self.changing_restraints = data["changing_restraints"] + self.orders = data["orders"] + self.simulation_data = data["simulation_data"] def ref_state_work( @@ -1599,7 +1600,7 @@ def ref_state_work( Returns ------- - RT * np.log(trans * orient): openff.units.unit.Quantity + RT * numpy.log(trans * orient): openff.units.unit.Quantity The free energy associated with releasing the restraints (in kcal/mol openff units). """ @@ -1613,41 +1614,45 @@ def ref_state_work( # Distance Integration Function def dist_int(RT, fc, targ): def potential(arange, RT, fc, targ): - return (arange**2) * np.exp((-1.0 / RT) * fc * (arange - targ) ** 2) + return (arange**2) * numpy.exp((-1.0 / RT) * fc * (arange - targ) ** 2) targ = targ.to(distance_unit) fc = fc.to(energy_unit / distance_unit**2) - arange = (np.arange(0.0, 100.0, 0.0001) * openff_unit.angstrom).to( + arange = (numpy.arange(0.0, 100.0, 0.0001) * openff_unit.angstrom).to( distance_unit ) - return np.trapz(potential(arange, RT, fc, targ), arange) + return numpy.trapz(potential(arange, RT, fc, targ), arange) # Angle Integration Function def ang_int(RT, fc, targ): def potential(arange, RT, fc, targ): - return np.sin(arange) * np.exp((-1.0 / RT) * fc * (arange - targ) ** 2) + return numpy.sin(arange) * numpy.exp( + (-1.0 / RT) * fc * (arange - targ) ** 2 + ) targ = targ.to(angle_unit) fc = fc.to(energy_unit / angle_unit**2) - arange = (np.arange(0.0, np.pi, 0.00005) * openff_unit.radians).to(angle_unit) + arange = (numpy.arange(0.0, numpy.pi, 0.00005) * openff_unit.radians).to( + angle_unit + ) - return np.trapz(potential(arange, RT, fc, targ), arange) + return numpy.trapz(potential(arange, RT, fc, targ), arange) # Torsion Integration Function def tors_int(RT, fc, targ): def potential(arange, RT, fc, targ): - return np.exp((-1.0 / RT) * fc * (arange - targ) ** 2) + return numpy.exp((-1.0 / RT) * fc * (arange - targ) ** 2) # Note, because of periodicity, I'm gonna wrap +/- pi around target for integration. targ = targ.to(angle_unit) fc = fc.to(energy_unit / angle_unit**2) arange = ( - np.arange(targ.magnitude - np.pi, targ.magnitude + np.pi, 0.00005) + numpy.arange(targ.magnitude - numpy.pi, targ.magnitude + numpy.pi, 0.00005) * openff_unit.radians ).to(angle_unit) - return np.trapz(potential(arange, RT, fc, targ), arange) + return numpy.trapz(potential(arange, RT, fc, targ), arange) # Distance restraint, r if None in [r_fc, r_tg]: @@ -1663,13 +1668,13 @@ def potential(arange, RT, fc, targ): # Torsion restraint, phi if None in [ph_fc, ph_tg]: - ph_int = 2.0 * np.pi * openff_unit.radians + ph_int = 2.0 * numpy.pi * openff_unit.radians else: ph_int = tors_int(RT, ph_fc, ph_tg) # Torsion restraint, alpha if None in [a_fc, a_tg]: - a_int = 2.0 * np.pi * openff_unit.radians + a_int = 2.0 * numpy.pi * openff_unit.radians else: a_int = tors_int(RT, a_fc, a_tg) @@ -1681,7 +1686,7 @@ def potential(arange, RT, fc, targ): # Torsion restraint, gamma if None in [g_fc, g_tg]: - g_int = 2.0 * np.pi * openff_unit.radians + g_int = 2.0 * numpy.pi * openff_unit.radians else: g_int = tors_int(RT, g_fc, g_tg) @@ -1690,8 +1695,8 @@ def potential(arange, RT, fc, targ): translational = r_int * th_int * ph_int * (1.0 / V0) # C^o = 1/V^o # Orientational term - rotational_volume = 8.0 * np.pi**2 + rotational_volume = 8.0 * numpy.pi**2 orientational = a_int * b_int * g_int / rotational_volume # Return the free energy - return RT * np.log(translational * orientational) + return RT * numpy.log(translational * orientational) diff --git a/paprika/analysis/bootstrap.py b/paprika/analysis/bootstrap.py index adbe1c7a..584b124f 100644 --- a/paprika/analysis/bootstrap.py +++ b/paprika/analysis/bootstrap.py @@ -1,8 +1,8 @@ from typing import Dict, List, Union -import numpy as np +import numpy import openmm.unit as openmm_unit -from openff.units import unit +from openff.units import unit as openff_unit from openff.units.openmm import from_openmm from scipy import stats as statistics from scipy.interpolate import Akima1DInterpolator @@ -20,13 +20,13 @@ def integrate_bootstraps(x, ys, x_intp=None, matrix="full"): Parameters ---------- - x: :class:`np.array` + x: :class:`numpy.array` The x coordinate of the curve to be integrated. - ys: :class:`np.array` + ys: :class:`numpy.array` Two dimensional array in which the first dimension is boot_cycles and the second dimension contains the arrays of y values which correspond to the x values and will be used for integration. The shape of this is :code:`(boot_cycles, len(x))`. - x_intp: :class:`np.array`, optional, default=None + x_intp: :class:`numpy.array`, optional, default=None An array which finely interpolates the x values. If not provided, it will be generated by adding 100 evenly spaced points between each x value. Default: None. matrix: str, optional, default='full` @@ -37,9 +37,9 @@ def integrate_bootstraps(x, ys, x_intp=None, matrix="full"): Returns ------- - avg_matrix: :class:`np.array` + avg_matrix: :class:`numpy.array` Matrix of the integration mean between each x value (as specified by 'matrix') - sem_matrix: :class:`np.array` + sem_matrix: :class:`numpy.array` Matrix of the uncertainty (SEM) between each x value (as specified by 'matrix') """ @@ -47,25 +47,25 @@ def integrate_bootstraps(x, ys, x_intp=None, matrix="full"): num_x = len(x) # Prepare to store the index location of the x values in the x_intp array - x_idxs = np.zeros([num_x], np.int32) + x_idxs = numpy.zeros([num_x], numpy.int32) # If not provided, generate x interpolation with 100 inpolated points between # each x value. Store the index locations of the x values in the x_intp # array. if x_intp is None: - x_intp = np.zeros([0], np.float64) + x_intp = numpy.zeros([0], numpy.float64) for i in range(1, num_x): - x_intp = np.append( - x_intp, np.linspace(x[i - 1], x[i], num=100, endpoint=False) + x_intp = numpy.append( + x_intp, numpy.linspace(x[i - 1], x[i], num=100, endpoint=False) ) x_idxs = len(x_intp) # Tack on the final value onto the interpolation - x_intp = np.append(x_intp, x[-1]) + x_intp = numpy.append(x_intp, x[-1]) # If x_intp is provided, find the locations of x values in x_intp else: i = 0 for j in range(len(x_intp)): - if np.isclose(x[i], x_intp[j]): + if numpy.isclose(x[i], x_intp[j]): x_idxs[i] = j i += 1 if i != num_x: @@ -77,7 +77,7 @@ def integrate_bootstraps(x, ys, x_intp=None, matrix="full"): cycles = len(ys) # Setup array to store integration bootstraps - int_matrix = np.zeros([num_x, num_x, cycles], np.float64) + int_matrix = numpy.zeros([num_x, num_x, cycles], numpy.float64) # Do the integration bootstraps. Originally, I had matrix=endpoints in the loop # below with everything else, but I'll split it out here in case that's faster @@ -88,8 +88,8 @@ def integrate_bootstraps(x, ys, x_intp=None, matrix="full"): y_intp = intp_func(x_intp) # for i in range(0, num_x): # for j in range(i+1, num_x): - # int_matrix[i, j, cycle] = np.trapz( y_intp, x_intp ) - int_matrix[0, num_x - 1, cycle] = np.trapz(y_intp, x_intp) + # int_matrix[i, j, cycle] = numpy.trapz( y_intp, x_intp ) + int_matrix[0, num_x - 1, cycle] = numpy.trapz(y_intp, x_intp) else: for cycle in range(cycles): intp_func = Akima1DInterpolator(x, ys[cycle]) @@ -100,12 +100,14 @@ def integrate_bootstraps(x, ys, x_intp=None, matrix="full"): continue beg = x_idxs[i] end = x_idxs[j] - int_matrix[i, j, cycle] = np.trapz(y_intp[beg:end], x_intp[beg:end]) + int_matrix[i, j, cycle] = numpy.trapz( + y_intp[beg:end], x_intp[beg:end] + ) # Setup matrices to store the average/sem values. # Is it bad that the default is 0.0 rather than None? - avg_matrix = np.zeros([num_x, num_x], np.float64) - sem_matrix = np.zeros([num_x, num_x], np.float64) + avg_matrix = numpy.zeros([num_x, num_x], numpy.float64) + sem_matrix = numpy.zeros([num_x, num_x], numpy.float64) # Second pass to compute the mean and standard deviation. for i in range(0, num_x): @@ -116,9 +118,9 @@ def integrate_bootstraps(x, ys, x_intp=None, matrix="full"): continue if matrix == "endpoints" and i != 0 and j != num_x - 1: continue - avg_matrix[i, j] = np.mean(int_matrix[i, j]) + avg_matrix[i, j] = numpy.mean(int_matrix[i, j]) avg_matrix[j, i] = -1.0 * avg_matrix[i, j] - sem_matrix[i, j] = np.std(int_matrix[i, j]) + sem_matrix[i, j] = numpy.std(int_matrix[i, j]) sem_matrix[j, i] = sem_matrix[i, j] return avg_matrix, sem_matrix @@ -155,56 +157,56 @@ def regression_bootstrap( Returns the mean, sem, ci_low, ci_high for [slope, intercept, R, R^2, RMSE, MSE, MUE, Kendall's Tau] """ - summary_statistics = np.empty((cycles, 8)) + summary_statistics = numpy.empty((cycles, 8)) # Bootstrapping for cycle in range(cycles): - new_x = np.empty_like(x_data) - new_y = np.empty_like(y_data) + new_x = numpy.empty_like(x_data) + new_y = numpy.empty_like(y_data) for index in range(len(x_data)): if with_replacement: - j = np.random.randint(len(x_data)) + j = numpy.random.randint(len(x_data)) else: j = index if with_uncertainty and x_sem is not None: - new_x[index] = np.random.normal(x_data[j], x_sem[j]) + new_x[index] = numpy.random.normal(x_data[j], x_sem[j]) elif with_uncertainty and x_sem is None: new_x[index] = x_data[j] if with_uncertainty and y_sem is not None: - new_y[index] = np.random.normal(y_data[j], y_sem[j]) + new_y[index] = numpy.random.normal(y_data[j], y_sem[j]) elif with_uncertainty and y_sem is None: new_y[index] = y_data[j] summary_statistics[cycle] = summarize_statistics(new_x, new_y) # Confidence interval - ci = np.empty((8, 2)) + ci = numpy.empty((8, 2)) for statistic in range(8): - sorted_statistic = np.sort(summary_statistics[:, statistic]) + sorted_statistic = numpy.sort(summary_statistics[:, statistic]) ci[statistic][0] = sorted_statistic[int(0.025 * cycles)] ci[statistic][1] = sorted_statistic[int(0.975 * cycles)] # Summarize results results = { "mean": { - "slope": np.mean(summary_statistics[:, 0]), - "intercept": np.mean(summary_statistics[:, 1]), - "R": np.mean(summary_statistics[:, 2]), - "R**2": np.mean(summary_statistics[:, 3]), - "RMSE": np.mean(summary_statistics[:, 4]), - "MSE": np.mean(summary_statistics[:, 5]), - "MUE": np.mean(summary_statistics[:, 6]), - "Tau": np.mean(summary_statistics[:, 7]), + "slope": numpy.mean(summary_statistics[:, 0]), + "intercept": numpy.mean(summary_statistics[:, 1]), + "R": numpy.mean(summary_statistics[:, 2]), + "R**2": numpy.mean(summary_statistics[:, 3]), + "RMSE": numpy.mean(summary_statistics[:, 4]), + "MSE": numpy.mean(summary_statistics[:, 5]), + "MUE": numpy.mean(summary_statistics[:, 6]), + "Tau": numpy.mean(summary_statistics[:, 7]), }, "sem": { - "slope": np.std(summary_statistics[:, 0]), - "intercept": np.std(summary_statistics[:, 1]), - "R": np.std(summary_statistics[:, 2]), - "R**2": np.std(summary_statistics[:, 3]), - "RMSE": np.std(summary_statistics[:, 4]), - "MSE": np.std(summary_statistics[:, 5]), - "MUE": np.std(summary_statistics[:, 6]), - "Tau": np.std(summary_statistics[:, 7]), + "slope": numpy.std(summary_statistics[:, 0]), + "intercept": numpy.std(summary_statistics[:, 1]), + "R": numpy.std(summary_statistics[:, 2]), + "R**2": numpy.std(summary_statistics[:, 3]), + "RMSE": numpy.std(summary_statistics[:, 4]), + "MSE": numpy.std(summary_statistics[:, 5]), + "MUE": numpy.std(summary_statistics[:, 6]), + "Tau": numpy.std(summary_statistics[:, 7]), }, "ci_low": { "slope": ci[0][0], @@ -232,14 +234,14 @@ def regression_bootstrap( def dG_bootstrap( - x_data: Union[List, np.array], - x_sem: Union[List, np.array], - y_data: Union[List, np.array], - y_sem: Union[List, np.array], + x_data: Union[List, numpy.array], + x_sem: Union[List, numpy.array], + y_data: Union[List, numpy.array], + y_sem: Union[List, numpy.array], cycles: int = 1000, - temperature: Union[float, unit.Quantity] = 298.15 * unit.kelvin, + temperature: Union[float, openff_unit.Quantity] = 298.15 * openff_unit.kelvin, with_uncertainty: bool = True, - energy_units: unit.Quantity = None, + energy_units: openff_unit.Quantity = None, ): """ Combine dG when for multiple binding poses. This is from Eq(A12) and Eq(A13) from @@ -261,7 +263,7 @@ def dG_bootstrap( Temperature of the simulation with_uncertainty: bool If true, generate samples from normal distribution based on mean=data, mu=sem - energy_units: unit.Quantity + energy_units: openff.units.unit.Quantity If true, return values as openff.units.unit.Quantity. Default is kcal/mol. Return @@ -270,41 +272,41 @@ def dG_bootstrap( Returns the mean, standard deviation and confidence interval from bootstrapping. """ - x_data = check_unit(x_data, base_unit=unit.kilocalorie_per_mole).magnitude - x_sem = check_unit(x_sem, base_unit=unit.kilocalorie_per_mole).magnitude - y_data = check_unit(y_data, base_unit=unit.kilocalorie_per_mole).magnitude - y_sem = check_unit(y_sem, base_unit=unit.kilocalorie_per_mole).magnitude - temperature = check_unit(temperature, base_unit=unit.kelvin) + x_data = check_unit(x_data, base_unit=openff_unit.kilocalorie_per_mole).magnitude + x_sem = check_unit(x_sem, base_unit=openff_unit.kilocalorie_per_mole).magnitude + y_data = check_unit(y_data, base_unit=openff_unit.kilocalorie_per_mole).magnitude + y_sem = check_unit(y_sem, base_unit=openff_unit.kilocalorie_per_mole).magnitude + temperature = check_unit(temperature, base_unit=openff_unit.kelvin) - summary_statistics = np.empty((cycles)) - RT = (R_gas * temperature).to(unit.kilocalorie_per_mole).magnitude + summary_statistics = numpy.empty((cycles)) + RT = (R_gas * temperature).to(openff_unit.kilocalorie_per_mole).magnitude beta = 1.0 / RT - ci = np.empty((2)) + ci = numpy.empty((2)) for cycle in range(cycles): - new_x = np.empty_like(x_data) - new_y = np.empty_like(y_data) + new_x = numpy.empty_like(x_data) + new_y = numpy.empty_like(y_data) if with_uncertainty and x_sem is not None: - new_x = np.random.normal(x_data, x_sem) + new_x = numpy.random.normal(x_data, x_sem) elif with_uncertainty and x_sem is None: new_x = x_data if with_uncertainty and y_sem is not None: - new_y = np.random.normal(y_data, y_sem) + new_y = numpy.random.normal(y_data, y_sem) elif with_uncertainty and y_sem is None: new_y = y_data - summary_statistics[cycle] = -RT * np.log( - np.exp(-beta * new_x) + np.exp(-beta * new_y) + summary_statistics[cycle] = -RT * numpy.log( + numpy.exp(-beta * new_x) + numpy.exp(-beta * new_y) ) # Get confidence interval - sorted_statistic = np.sort(summary_statistics) + sorted_statistic = numpy.sort(summary_statistics) ci[0] = sorted_statistic[int(0.025 * cycles)] ci[1] = sorted_statistic[int(0.975 * cycles)] results = { - "mean": np.mean(summary_statistics), - "sem": np.std(summary_statistics), + "mean": numpy.mean(summary_statistics), + "sem": numpy.std(summary_statistics), "ci": ci, } @@ -325,9 +327,9 @@ def dH_bootstrap( dG_y_data, dG_y_sem, cycles=1000, - temperature=298.15 * unit.kelvin, + temperature=298.15 * openff_unit.kelvin, with_uncertainty=True, - energy_units: unit.Quantity = None, + energy_units: openff_unit.Quantity = None, ): """ Combine dH when for multiple binding poses. This is from Eq(A16) and Eq(A17) from @@ -357,7 +359,7 @@ def dH_bootstrap( Temperature of the simulation with_uncertainty: bool If true, generate samples from normal distribution based on mean=data, mu=sem - energy_units: unit.Quantity + energy_units: openff.units.unit.Quantity If true, return values as openff.units.unit.Quantity. Default is kcal/mol. Return @@ -366,60 +368,77 @@ def dH_bootstrap( Returns the mean, standard deviation and confidence interval from bootstrapping. """ - dH_x_data = check_unit(dH_x_data, base_unit=unit.kilocalorie_per_mole).magnitude - dH_x_sem = check_unit(dH_x_sem, base_unit=unit.kilocalorie_per_mole).magnitude - dH_y_data = check_unit(dH_y_data, base_unit=unit.kilocalorie_per_mole).magnitude - dH_y_sem = check_unit(dH_y_sem, base_unit=unit.kilocalorie_per_mole).magnitude - dG_x_data = check_unit(dG_x_data, base_unit=unit.kilocalorie_per_mole).magnitude - dG_x_sem = check_unit(dG_x_sem, base_unit=unit.kilocalorie_per_mole).magnitude - dG_y_data = check_unit(dG_y_data, base_unit=unit.kilocalorie_per_mole).magnitude - dG_y_sem = check_unit(dG_y_sem, base_unit=unit.kilocalorie_per_mole).magnitude - temperature = check_unit(temperature, base_unit=unit.kelvin) - - summary_statistics = np.empty((cycles)) - RT = (R_gas * temperature).to(unit.kcal / unit.mole).magnitude + dH_x_data = check_unit( + dH_x_data, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + dH_x_sem = check_unit( + dH_x_sem, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + dH_y_data = check_unit( + dH_y_data, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + dH_y_sem = check_unit( + dH_y_sem, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + dG_x_data = check_unit( + dG_x_data, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + dG_x_sem = check_unit( + dG_x_sem, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + dG_y_data = check_unit( + dG_y_data, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + dG_y_sem = check_unit( + dG_y_sem, base_unit=openff_unit.kilocalorie_per_mole + ).magnitude + temperature = check_unit(temperature, base_unit=openff_unit.kelvin) + + summary_statistics = numpy.empty((cycles)) + RT = (R_gas * temperature).to(openff_unit.kcal / openff_unit.mole).magnitude beta = 1.0 / RT - ci = np.empty((2)) + ci = numpy.empty((2)) for cycle in range(cycles): - new_dH_x = np.empty_like(dH_x_data) - new_dH_y = np.empty_like(dH_y_data) + new_dH_x = numpy.empty_like(dH_x_data) + new_dH_y = numpy.empty_like(dH_y_data) - new_dG_x = np.empty_like(dG_x_data) - new_dG_y = np.empty_like(dG_y_data) + new_dG_x = numpy.empty_like(dG_x_data) + new_dG_y = numpy.empty_like(dG_y_data) # Resample dH if with_uncertainty and dH_x_sem is not None: - new_dH_x = np.random.normal(dH_x_data, dH_x_sem) + new_dH_x = numpy.random.normal(dH_x_data, dH_x_sem) elif with_uncertainty and dH_x_sem is None: new_dH_x = dH_x_data if with_uncertainty and dH_y_sem is not None: - new_dH_y = np.random.normal(dH_y_data, dH_y_sem) + new_dH_y = numpy.random.normal(dH_y_data, dH_y_sem) elif with_uncertainty and dH_y_sem is None: new_dH_y = dH_y_data # Resample dG if with_uncertainty and dG_x_sem is not None: - new_dG_x = np.random.normal(dG_x_data, dG_x_sem) + new_dG_x = numpy.random.normal(dG_x_data, dG_x_sem) elif with_uncertainty and dG_x_sem is None: new_dG_x = dG_x_data if with_uncertainty and dG_y_sem is not None: - new_dG_y = np.random.normal(dG_y_data, dG_y_sem) + new_dG_y = numpy.random.normal(dG_y_data, dG_y_sem) elif with_uncertainty and dG_y_sem is None: new_dG_y = dG_y_data summary_statistics[cycle] = ( - new_dH_x * np.exp(-beta * new_dG_x) + new_dH_y * np.exp(-beta * new_dG_y) - ) / (np.exp(-beta * new_dG_x) + np.exp(-beta * new_dG_y)) + new_dH_x * numpy.exp(-beta * new_dG_x) + + new_dH_y * numpy.exp(-beta * new_dG_y) + ) / (numpy.exp(-beta * new_dG_x) + numpy.exp(-beta * new_dG_y)) # Confidence interval - sorted_statistic = np.sort(summary_statistics) + sorted_statistic = numpy.sort(summary_statistics) ci[0] = sorted_statistic[int(0.025 * cycles)] ci[1] = sorted_statistic[int(0.975 * cycles)] results = { - "mean": np.mean(summary_statistics), - "sem": np.std(summary_statistics), + "mean": numpy.mean(summary_statistics), + "sem": numpy.std(summary_statistics), "ci": ci, } @@ -436,9 +455,9 @@ def summarize_statistics(x, y): Parameters ---------- - x: np.array + x: numpy.array X data - y: np.array + y: numpy.array Y data Returns @@ -453,7 +472,7 @@ def summarize_statistics(x, y): * MUE - Mean Unsigned Error or Mean Absolute Error * Tau - Kendall's Tau """ - summary_statistics = np.empty(8) + summary_statistics = numpy.empty(8) # Slope, intercept, R - Pearson correlation coefficient ( summary_statistics[0], @@ -467,13 +486,13 @@ def summarize_statistics(x, y): summary_statistics[3] = summary_statistics[2] ** 2 # RMSE - Root-Mean-Squared-Error - summary_statistics[4] = np.sqrt(np.mean((y - x) ** 2)) + summary_statistics[4] = numpy.sqrt(numpy.mean((y - x) ** 2)) # MSE - Mean Signed Error - summary_statistics[5] = np.mean(y - x) + summary_statistics[5] = numpy.mean(y - x) # MUE - Mean Unsigned Error - summary_statistics[6] = np.mean(np.absolute(y - x)) + summary_statistics[6] = numpy.mean(numpy.absolute(y - x)) # Tau - Kendall's Tau summary_statistics[7], prob = statistics.kendalltau(x, y) diff --git a/paprika/analysis/utils.py b/paprika/analysis/utils.py index 989792d9..a2afca86 100644 --- a/paprika/analysis/utils.py +++ b/paprika/analysis/utils.py @@ -1,4 +1,4 @@ -import numpy as np +import numpy def get_factors(n): @@ -17,7 +17,7 @@ def get_factors(n): """ factors = [] - sqrt_n = int(round(np.sqrt(n) + 0.5)) + sqrt_n = int(round(numpy.sqrt(n) + 0.5)) i = 1 while i <= sqrt_n: if n % i == 0: @@ -26,6 +26,7 @@ def get_factors(n): if j != i: factors.append(int(j)) i += 1 + return sorted(factors, key=int) @@ -45,19 +46,24 @@ def get_nearest_max(n): """ max_factors = 0 + if n % 2 == 0: beg = n - 100 end = n else: beg = n - 101 end = n - 1 + if beg < 0: beg = 0 + + most_factors = 0 for i in range(beg, end + 2, 2): num_factors = len(get_factors(i)) if num_factors >= max_factors: max_factors = num_factors most_factors = i + return most_factors @@ -72,12 +78,12 @@ def get_block_sem(data_array): Parameters ---------- - data_array: :class:`np.array` + data_array: :class:`numpy.array` Array containing data values. Returns ------- - np.max(sems): float + numpy.max(sems): float The maximum SEM obtained from te blocking curve. """ @@ -86,11 +92,11 @@ def get_block_sem(data_array): block_sizes = get_factors(len(data_array)) # An array to store means for each block ... make it bigger than we need. - block_means = np.zeros([block_sizes[-1]], np.float64) + block_means = numpy.zeros([block_sizes[-1]], numpy.float64) # Store the SEM for each block size, except the last two size for which # there will only be two or one blocks total and thus very noisy. - sems = np.zeros([len(block_sizes) - 2], np.float64) + sems = numpy.zeros([len(block_sizes) - 2], numpy.float64) # Check each block size except the last two. for size_idx in range(len(block_sizes) - 2): @@ -102,15 +108,15 @@ def get_block_sem(data_array): data_beg_idx = blk_idx * block_sizes[size_idx] data_end_idx = (blk_idx + 1) * block_sizes[size_idx] # Compute the mean of this block and store in array - block_means[blk_idx] = np.mean(data_array[data_beg_idx:data_end_idx]) + block_means[blk_idx] = numpy.mean(data_array[data_beg_idx:data_end_idx]) # Compute the standard deviation across all blocks, devide by # num_blocks-1 for SEM - sems[size_idx] = np.std(block_means[0:num_blocks], ddof=0) / np.sqrt( + sems[size_idx] = numpy.std(block_means[0:num_blocks], ddof=0) / numpy.sqrt( num_blocks - 1 ) # Hmm or should ddof=1? I think 0, see Flyvbjerg -----^ - return np.max(sems) + return numpy.max(sems) def get_subsampled_indices(N, g, conservative=False): @@ -138,16 +144,16 @@ def get_subsampled_indices(N, g, conservative=False): # if conservative, assume integer g and round up if conservative: - g = np.ceil(g) + g = numpy.ceil(g) # initialize indices = [0] g_idx = 1.0 - int_step = int(np.round(g_idx * g)) + int_step = int(numpy.round(g_idx * g)) while int_step < N: indices.append(int_step) g_idx += 1.0 - int_step = int(np.round(g_idx * g)) + int_step = int(numpy.round(g_idx * g)) return indices diff --git a/paprika/build/system/__init__.py b/paprika/build/system/__init__.py index 148d7b7d..448b125d 100644 --- a/paprika/build/system/__init__.py +++ b/paprika/build/system/__init__.py @@ -1,8 +1,10 @@ from paprika.build.system.tleap import TLeap from paprika.build.system.utils import ConversionToolkit, PBCBox +from paprika.build.system.taproom import BuildTaproomAPR __all__ = [ "ConversionToolkit", "PBCBox", "TLeap", + "BuildTaproomAPR", ] diff --git a/paprika/build/system/taproom.py b/paprika/build/system/taproom.py new file mode 100644 index 00000000..74079122 --- /dev/null +++ b/paprika/build/system/taproom.py @@ -0,0 +1,931 @@ +import json +import os +import shutil +from typing import Any, Dict, List, Union + +import numpy +import openmm +import openmm.app as app +import parmed +from joblib import Parallel, delayed +from openff.interchange import Interchange +from openff.interchange.components._packmol import pack_box +from openff.toolkit import ForceField, Molecule, Topology +from openff.units import unit +from tqdm.auto import tqdm + +from paprika.build import align +from paprika.evaluator import Setup +from paprika.io import PaprikaEncoder, save_restraints +from paprika.restraints import DAT_restraint, create_window_list, parse_window +from paprika.restraints.openmm import apply_dat_restraint, apply_positional_restraints +from paprika.taproom import get_benchmarks, read_yaml_schema +from paprika.utils import index_from_mask + + +class BuildTaproomAPR: + """A class to generate APR files from Taproom database in pAPRika. + + As of now this class can generate files for a single host-guest pair with + explicit and implicit solvent based on OpenFF specifications. + + TODO: Implement an option to build an array of host-guest pairs from Taproom. + TODO: Implement an option to use the GAFF force field (possible through the OpenMMForceFields package). + + Parameters + ---------- + host_code: str + The 3-letter Taproom code for the host molecule. + guest_code: str + The 3-letter Taproom code for the guest molecule. + host_guest_codes: dict + Selection of the host-guest pairs as a dictionary. + n_water: int + Number of Water molecules. If set as `0` or `None` then the system will be built without water. + build_folder: str + Temporary folder to save intermediate files. + working_folder: str + The main folder to write the APR structure files. + + Examples + -------- + >>> from paprika.build.system import BuildTaproomAPR + >>> from openff.toolkit import ForceField + >>> + >>> # Select OpenFF force field version 2.0.0 + >>> force_field = ForceField("openff-2.0.0.offxml") + >>> + >>> # Initiate system object + >>> system = BuildTaproomAPR(host_code="bcd", guest_code="hex", n_water=3000, force_field=force_field) + >>> + >>> # Build APR files + >>> system.build_system() + >>> + >>> # We can also extend the default `r_final` specified in Taproom if it's not far enough + >>> from openff.units import unit + >>> system.extend_pull_distance(extend_by=6*unit.angstrom) + >>> system.build_system() + >>> + >>> # We can create the APR system in a vacuum (useful if you want to add your own custom implicit solvent later on) + >>> system = BuildTaproomAPR(host_code="bcd", guest_code="hex", n_water=None, force_field=force_field) + >>> system.build_system() + >>> + >>> # We can create the APR system with the OBC2 implicit solvent + >>> from pkg_resources import resource_filename + >>> GBSA = resource_filename( + >>> " openff.toolkit", + >>> os.path.join("data", "test_forcefields", "GBSA_OBC2-1.0.offxml"), + >>> ) + >>> force_field = ForceField("openff-2.0.0.offxml", GBSA) + >>> + >>> system = BuildTaproomAPR(host_code="bcd", guest_code="hex", n_water=None, force_field=force_field) + >>> system.build_system() + >>> + >>> # Creating these files can take 10-20 minutes on one core (especially with explicit solvent). + >>> # We can speed things up running these in parallel. + >>> system.build_system(n_cpus=4) + """ + + def __init__( + self, + host_code: Union[str, None] = None, + guest_code: Union[str, None] = None, + host_guest_codes: Union[Dict[str, List[str]], None] = None, + n_water: Union[int, None] = None, + force_field: Union[ForceField, None] = None, + use_taproom_mol2: bool = True, + build_folder: str = "build_files", + working_folder: str = "simulations", + disable_progress: bool = True, + ): + if force_field is None: + raise ValueError( + "The option `force_field` cannot be a None. Please specify an OpenFF `ForceField` object." + ) + + self._host_code = host_code + self._guest_code = guest_code + self._host_guest_code = host_guest_codes + self._n_water = n_water + if n_water == 0 or n_water is None: + self._n_water = None + self._use_taproom_mol2 = use_taproom_mol2 + self._force_field = force_field + self._build_folder = build_folder + self._working_folder = working_folder + self._disable_progress = disable_progress + self._host_metadata = None + self._guest_metadata = None + self._orientations = None + self._water_mol = None + self._water_intrcg = None + self._restraints = None + self._n_cpus = 1 + self._build_array = False + + if self._host_code is not None and self._guest_code is not None: + self._initialize_single() + else: + self._build_array = True + raise NotImplementedError( + "Creating APR files for an array of host-guest pairs is not implemented yet." + ) + + def _initialize_single(self): + # Create folder + os.makedirs(self._build_folder, exist_ok=True) + os.makedirs(self._working_folder, exist_ok=True) + + # Load Host-Guest system from Taproom + taproom = get_benchmarks() + self._host_metadata = taproom["host_guest_systems"][self._host_code] + self._guest_metadata = taproom["host_guest_systems"][self._host_code][ + self._guest_code + ] + self._orientations = list(self._host_metadata["yaml"].keys()) + self._host_yaml_schema = read_yaml_schema( + self._host_metadata["yaml"][self._orientations[0]] + ) + self._guest_yaml_schema = read_yaml_schema(self._guest_metadata["yaml"]) + + # Water Molecule + if self._n_water is not None: + self._water_mol = Molecule.from_smiles("O") + self._water_intrcg = Interchange.from_smirnoff( + force_field=self._force_field, + topology=[self._water_mol] * self._n_water, + ) + + # Load restraints + self._restraints = { + "static_restraints": self._unnest_restraint_specs( + self._host_yaml_schema["restraints"]["static"] + ), + "conformational_restraints": self._unnest_restraint_specs( + self._host_yaml_schema["restraints"]["conformational"] + ), + "guest_restraints": self._unnest_restraint_specs( + self._guest_yaml_schema["restraints"]["guest"] + ), + "wall_restraints": self._unnest_restraint_specs( + self._guest_yaml_schema["restraints"]["wall_restraints"] + ), + "symmetry_restraints": self._unnest_restraint_specs( + self._guest_yaml_schema["symmetry_correction"]["restraints"] + ), + } + + def _initialize_array(self): + """Not implemented yet.""" + pass + + @staticmethod + def _unnest_restraint_specs( + restraint_specs: List[Dict[str, Any]] + ) -> List[Dict[str, Any]]: + """A helper method to un-nest restraint lists parsed from a taproom + yaml file. + + Parameters + ---------- + restraint_specs + The restraint specs to un-nest. + """ + return [ + value["restraint"] + for value in restraint_specs + if value["restraint"] is not None + ] + + @staticmethod + def _restraints_to_dict(restraints: List[DAT_restraint]): + """Converts a list of ``paprika`` restraint objects to + a list of JSON compatible dictionary representations + """ + + return [ + json.loads(json.dumps(restraint.__dict__, cls=PaprikaEncoder)) + for restraint in restraints + ] + + def _solvate_and_add_dummy( + self, + complex_path: str, + solvated_path: str, + host_resname: str, + unique_molecules: List[Molecule], + initial_distance: float, + offset_mask: Union[str, None] = None, + ): + """Solvate a PDB file with PackMol through OpenFF-Interchange. + + Parameters + ---------- + complex_path: str + The file path of the complex PDB. + solvated_path: str + The output file path for the solvated complex. + unique_molecules: List[Molecule] + List of unique molecules (to generate the OpenFF Topology) + initial_distance: float + An offset to place the dummy atoms. + offset_mask: str + The AMBER atom mask for G1 - used for offsetting Dummy atoms + + Returns + ------- + system_intrcg: openff.interchange.Interchange + The solvated system as OpenFF-Interchange object. + """ + rectangular_box = numpy.asarray( + [ + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 2.0], + ] + ) + + # 01 - Solvate structure + pdbfile = app.PDBFile(complex_path) + solute_topology = Topology.from_openmm( + pdbfile.topology, unique_molecules=unique_molecules + ) + solute_intrcg = Interchange.from_smirnoff( + force_field=self._force_field, + topology=solute_topology, + charge_from_molecules=unique_molecules if self._use_taproom_mol2 else None, + ) + if self._n_water is not None: + solvated_topology = pack_box( + molecules=[self._water_mol], + number_of_copies=[self._n_water], + solute=solute_intrcg.topology, + box_shape=rectangular_box, + mass_density=0.95 * unit.grams / unit.milliliters, + center_solute="ORIGIN", + ) + solute_intrcg.box = solvated_topology.box_vectors + self._water_intrcg.box = solvated_topology.box_vectors + solvated_topology.to_file(solvated_path) + else: + solute_intrcg.box = None + solute_intrcg.topology.to_file(solvated_path) + + # 02 - Add Dummy Atoms to PDB + input_structure = parmed.load_file( + solvated_path if self._n_water is not None else complex_path, + structure=True, + ) + if self._n_water is None: + input_structure = align.translate_to_origin( + input_structure, atom_mask=f":{host_resname}", weight="geo" + ) + + offset_array = numpy.array([0.0, 0.0, 0.0]) + if offset_mask is not None: + coordinates_z = input_structure[offset_mask].coordinates[-1][-1] + offset_array = numpy.array([0.0, 0.0, coordinates_z]) + + Setup.add_dummy_atoms_to_structure( + input_structure, + dummy_atom_offsets=[ + numpy.array([0, 0, -initial_distance]), + numpy.array([0, 0, -3.0 - initial_distance]), + numpy.array([0, 2.2, -5.2 - initial_distance]), + ], + offset_coordinates=offset_array, + ) + + # 03 - Shift structure to avoid issues with PBC + if self._n_water is not None: + input_structure.coordinates += numpy.array( + [ + input_structure.box[0] * 0.5, + input_structure.box[1] * 0.5, + -input_structure.coordinates[-1, 2] + 5.0, + ] + ) + + # 04 - Write PDB for solvated system + with open(solvated_path, "w") as f: + app.PDBFile.writeFile( + input_structure.topology, + input_structure.positions, + f, + keepIds=True, + ) + + # 05 - Combine interchange objects + if self._n_water is not None: + system_intrcg = solute_intrcg + self._water_intrcg + system_intrcg.box = solvated_topology.box_vectors + else: + system_intrcg = solute_intrcg + system_intrcg.box = None + + return system_intrcg + + @staticmethod + def _create_system_and_add_dummy( + system_intrcg: Interchange, system_output_path: str + ): + """Convert `Interchange` object to OpenMM System and add dummy atoms. + + Parameters + ---------- + system_intrcg: Interchange + The Interchange object to convert. + system_output_path: str + The file path to save the OpenMM System to XML file. + """ + openmm_system = system_intrcg.to_openmm() + + for _ in range(3): + openmm_system.addParticle(mass=207) + + for i, force in enumerate(openmm_system.getForces()): + if isinstance(force, openmm.NonbondedForce): + force.addParticle(0.0, 1.0, 0.0) + force.addParticle(0.0, 1.0, 0.0) + force.addParticle(0.0, 1.0, 0.0) + + with open(system_output_path, "w") as f: + f.write(openmm.XmlSerializer.serialize(openmm_system)) + + def _build_pull_structures( + self, + i, + orient, + complex_path, + guest_atom_indices, + guest_orientation_mask, + pulling_distance, + host_resname, + initial_distance, + n_windows, + unique_molecules, + ): + """Function that translates guest molecules from host that is to be wrapped in `delayed` for parallelism.""" + + folder = f"{self._working_folder}/pull-{orient}/p{i:03}" + os.makedirs(folder, exist_ok=True) + + # 01 - Prepare complex structure + structure = Setup.prepare_complex_structure( + complex_path, + guest_atom_indices, + guest_orientation_mask, + pull_distance=pulling_distance, + pull_window_index=i, + n_pull_windows=n_windows["pull"], + ) + complex_prepared_path = ( + f"{folder}/{self._host_code}-{self._guest_code}-{orient}.pdb" + ) + with open(complex_prepared_path, "w") as f: + app.PDBFile.writeFile( + structure.topology, + structure.positions, + f, + keepIds=True, + ) + + # 02 - Solvate structure + r_i = numpy.linspace(0.0, pulling_distance, n_windows["pull"])[i] + complex_solvated_path = f"{folder}/restrained.pdb" + offset_mask = guest_orientation_mask.split(" ")[0] + host_guest_system_intrcg = self._solvate_and_add_dummy( + complex_prepared_path, + complex_solvated_path, + host_resname=host_resname, + unique_molecules=unique_molecules, + initial_distance=initial_distance + r_i, + offset_mask=offset_mask, + ) + + if i == 0: + # 03 - Create Host-Guest OpenMM System with Dummy Atoms + system_output_path = f"{self._build_folder}/{self._host_code}-{self._guest_code}-dum-solv.xml" + self._create_system_and_add_dummy( + host_guest_system_intrcg, system_output_path + ) + + # 04 - Clean up + os.remove(complex_prepared_path) + + def _build_apr_structures(self): + """Build and prepare the APR structures and windows.""" + host_resname = self._host_yaml_schema["resname"] + n_windows = self._host_yaml_schema["calculation"]["windows"] + + # Create OpenFF Molecule instances of molecules + guest_mol = Molecule.from_file( + str( + self._host_metadata["path"] + .joinpath(self._guest_yaml_schema["name"]) + .joinpath(self._guest_yaml_schema["structure"]["sdf"]) + ) + ) + host_mol = Molecule.from_file( + str( + self._host_metadata["path"].joinpath( + self._host_yaml_schema["structure"]["sdf"] + ) + ) + ) + + # --------------------------------------------------------------------- # + # Prepare Host-Guest Complex + # --------------------------------------------------------------------- # + # print("Generating files for the `pull` phase.") + for orient in self._orientations: + # 01 - Load complex structure + complex_path = str( + self._host_metadata["path"] + .joinpath(self._guest_yaml_schema["name"]) + .joinpath(self._guest_yaml_schema["complex"]) + ).replace(".pdb", f"-{orient}.pdb") + structure = parmed.load_file(complex_path, structure=True) + + # 02 - Get Guest indices and mask + guest_atom_indices = index_from_mask( + structure, f":{self._guest_yaml_schema['name'].upper()}" + ) + G1 = self._guest_yaml_schema["aliases"][3]["G1"] + G2 = self._guest_yaml_schema["aliases"][4]["G2"] + guest_orientation_mask = f"{G1} {G2}" + + # 03 - Initial `r` values + r_initial = self._guest_yaml_schema["restraints"]["guest"][0]["restraint"][ + "attach" + ]["target"] + r_final = self._guest_yaml_schema["restraints"]["guest"][0]["restraint"][ + "pull" + ]["target"] + pulling_distance = (r_final - r_initial).m_as(unit.angstrom) + initial_distance = self._guest_yaml_schema["restraints"]["guest"][0][ + "restraint" + ]["attach"]["target"].m_as(unit.angstrom) + + # --------------------------------------------------------------------- # + # Prepare `pull` windows + # --------------------------------------------------------------------- # + Parallel(n_jobs=self._n_cpus)( + delayed(self._build_pull_structures)( + i, + orient, + complex_path, + guest_atom_indices, + guest_orientation_mask, + pulling_distance, + host_resname, + initial_distance, + n_windows, + [host_mol, guest_mol], + ) + for i in tqdm(range(n_windows["pull"]), disable=self._disable_progress) + ) + + # --------------------------------------------------------------------- # + # Prepare `attach` windows - Copy PDB from p000 + # --------------------------------------------------------------------- # + # print("Generating files for the `attach` phase.") + for i in tqdm(range(n_windows["attach"]), disable=self._disable_progress): + folder = f"{self._working_folder}/attach-{orient}/a{i:03}" + os.makedirs(folder, exist_ok=True) + shutil.copy( + f"{self._working_folder}/pull-{orient}/p000/restrained.pdb", + f"{folder}/restrained.pdb", + ) + + # --------------------------------------------------------------------- # + # Prepare Host-only Structure + # --------------------------------------------------------------------- # + # print("Generating files for the `release` phase.") + # 01 - Remove Guest molecule and Dummy atoms from complex structure + complex_solvate_path = f"{self._working_folder}/pull-p/p000/restrained.pdb" + host_pdb = app.PDBFile( + str( + self._host_metadata["path"].joinpath( + self._host_yaml_schema["structure"]["pdb"] + ) + ) + ) + pdbfile = app.PDBFile(complex_solvate_path) + modeller = app.Modeller(pdbfile.topology, pdbfile.positions) + guest_atoms = [ + atom + for atom in pdbfile.topology.atoms() + if atom.residue.name != host_resname + ] + modeller.delete(guest_atoms) + + with open(f"{self._build_folder}/host_input.pdb", "w") as f: + app.PDBFile.writeFile( + host_pdb.topology, + modeller.positions, + f, + keepIds=False, + ) + + # 02 - Align host molecule + host_structure = Setup.prepare_host_structure( + f"{self._build_folder}/host_input.pdb" + ) + output_coordinate_path = f"{self._build_folder}/host_input_aligned.pdb" + with open(output_coordinate_path, "w") as file: + app.PDBFile.writeFile( + host_structure.topology, host_structure.positions, file, True + ) + + # 03 - Solvate host molecule and add dummy atoms + host_solvated_path = f"{self._build_folder}/{self._host_code}-dum-solv.pdb" + host_system_intrcg = self._solvate_and_add_dummy( + output_coordinate_path, + host_solvated_path, + host_resname=host_resname, + unique_molecules=[host_mol], + initial_distance=initial_distance, + offset_mask=None, + ) + + # 04 - Create Host-only OpenMM System with Dummy Atoms + system_output_path = f"{self._build_folder}/{self._host_code}-dum-solv.xml" + self._create_system_and_add_dummy(host_system_intrcg, system_output_path) + + # --------------------------------------------------------------------- # + # Prepare `release` windows - Copy PDB + # --------------------------------------------------------------------- # + for i in tqdm(range(n_windows["release"]), disable=self._disable_progress): + folder = f"{self._working_folder}/release/r{i:03}" + os.makedirs(folder, exist_ok=True) + shutil.copy( + host_solvated_path, + f"{folder}/restrained.pdb", + ) + + def _apply_attach_restraints(self): + """Apply restraints for the `attach` phase.""" + + # print("Applying restraints for `attach` phase...") + attach_lambdas = self._host_yaml_schema["calculation"]["lambda"]["attach"] + n_windows = self._host_yaml_schema["calculation"]["windows"] + + for orient in self._orientations: + attach_folder = f"{self._working_folder}/attach-{orient}" + complex_path = f"{self._working_folder}/pull-{orient}/p000/restrained.pdb" + + static_restraints = Setup.build_static_restraints( + complex_path, + n_attach_windows=n_windows["attach"], + n_pull_windows=None, + n_release_windows=None, + restraint_schemas=self._restraints["static_restraints"], + ) + conformational_restraints = Setup.build_conformational_restraints( + complex_path, + attach_lambdas=attach_lambdas, + n_pull_windows=None, + release_lambdas=None, + restraint_schemas=self._restraints["conformational_restraints"], + ) + guest_restraints = Setup.build_guest_restraints( + complex_path, + attach_lambdas=attach_lambdas, + n_pull_windows=None, + restraint_schemas=self._restraints["guest_restraints"], + ) + symmetry_restraints = Setup.build_symmetry_restraints( + complex_path, + n_attach_windows=n_windows["attach"], + restraint_schemas=self._restraints["symmetry_restraints"], + ) + wall_restraints = Setup.build_wall_restraints( + complex_path, + n_attach_windows=n_windows["attach"], + restraint_schemas=self._restraints["wall_restraints"], + ) + + symmetry_restraints = ( + [] if symmetry_restraints is None else symmetry_restraints + ) + wall_restraints = [] if wall_restraints is None else wall_restraints + guest_restraints = [] if guest_restraints is None else guest_restraints + + restraints_dictionary = { + "static": self._restraints_to_dict(static_restraints), + "conformational": self._restraints_to_dict(conformational_restraints), + "symmetry": self._restraints_to_dict(symmetry_restraints), + "wall": self._restraints_to_dict(wall_restraints), + "guest": self._restraints_to_dict(guest_restraints), + } + + with open(f"{attach_folder}/restraints.json", "w") as file: + json.dump(restraints_dictionary, file) + + save_restraints( + conformational_restraints + guest_restraints, + filepath=f"{attach_folder}/apr_restraints.json", + ) + + # Apply restraints + system_path = f"{self._build_folder}/{self._host_code}-{self._guest_code}-dum-solv.xml" + attach_windows = create_window_list(guest_restraints) + Parallel(n_jobs=self._n_cpus)( + delayed(self._apply_attach_to_system)( + window, + complex_path, + system_path, + attach_folder, + static_restraints, + conformational_restraints, + guest_restraints, + symmetry_restraints, + wall_restraints, + ) + for window in tqdm(attach_windows, disable=self._disable_progress) + ) + + def _apply_pull_restraints(self): + """Apply restraints for the `pull` phase.""" + + # print("Applying restraints for `pull` phase...") + attach_lambdas = self._host_yaml_schema["calculation"]["lambda"]["attach"] + n_windows = self._host_yaml_schema["calculation"]["windows"] + + for orient in self._orientations: + pull_folder = f"{self._working_folder}/pull-{orient}" + complex_path = f"{self._working_folder}/pull-{orient}/p000/restrained.pdb" + + static_restraints = Setup.build_static_restraints( + complex_path, + n_attach_windows=n_windows["attach"], + n_pull_windows=n_windows["pull"], + n_release_windows=None, + restraint_schemas=self._restraints["static_restraints"], + ) + conformational_restraints = Setup.build_conformational_restraints( + complex_path, + attach_lambdas=attach_lambdas, + n_pull_windows=n_windows["pull"], + release_lambdas=None, + restraint_schemas=self._restraints["conformational_restraints"], + ) + guest_restraints = Setup.build_guest_restraints( + complex_path, + attach_lambdas=attach_lambdas, + n_pull_windows=n_windows["pull"], + restraint_schemas=self._restraints["guest_restraints"], + ) + guest_restraints = [] if guest_restraints is None else guest_restraints + + # Remove the `attach` phases from the restraints as these restraints are + # only being used for the pull phase. + for restraint in ( + static_restraints + conformational_restraints + guest_restraints + ): + for key in restraint.phase["attach"]: + restraint.phase["attach"][key] = None + + restraints_dictionary = { + "static": self._restraints_to_dict(static_restraints), + "conformational": self._restraints_to_dict(conformational_restraints), + "symmetry": None, + "wall": None, + "guest": self._restraints_to_dict(guest_restraints), + } + + with open(f"{pull_folder}/restraints.json", "w") as file: + json.dump(restraints_dictionary, file) + + save_restraints( + conformational_restraints + guest_restraints, + filepath=f"{pull_folder}/apr_restraints.json", + ) + + # Apply restraints + system_path = f"{self._build_folder}/{self._host_code}-{self._guest_code}-dum-solv.xml" + pull_windows = create_window_list(guest_restraints) + Parallel(n_jobs=self._n_cpus)( + delayed(self._apply_pull_to_system)( + window, + complex_path, + system_path, + pull_folder, + static_restraints, + conformational_restraints, + guest_restraints, + ) + for window in tqdm(pull_windows, disable=self._disable_progress) + ) + + def _apply_release_restraints(self): + """Apply restraints for the `release` phase.""" + + # print("Applying restraints for `release` phase...") + release_lambdas = self._host_yaml_schema["calculation"]["lambda"]["release"] + n_windows = self._host_yaml_schema["calculation"]["windows"] + + release_folder = f"{self._working_folder}/release" + host_solvated_path = f"{release_folder}/r000/restrained.pdb" + + static_restraints = Setup.build_static_restraints( + host_solvated_path, + n_attach_windows=None, + n_pull_windows=None, + n_release_windows=n_windows["release"], + restraint_schemas=self._restraints["static_restraints"], + ) + conformational_restraints = Setup.build_conformational_restraints( + host_solvated_path, + attach_lambdas=None, + n_pull_windows=None, + release_lambdas=release_lambdas, + restraint_schemas=self._restraints["conformational_restraints"], + ) + + restraints_dictionary = { + "static": self._restraints_to_dict(static_restraints), + "conformational": self._restraints_to_dict(conformational_restraints), + "symmetry": None, + "wall": None, + "guest": None, + } + + with open(f"{release_folder}/restraints.json", "w") as file: + json.dump(restraints_dictionary, file) + + save_restraints( + conformational_restraints, + filepath=f"{release_folder}/apr_restraints.json", + ) + + # Apply restraints + system_path = f"{self._build_folder}/{self._host_code}-dum-solv.xml" + release_windows = create_window_list(conformational_restraints) + Parallel(n_jobs=self._n_cpus)( + delayed(self._apply_release_to_system)( + window, + host_solvated_path, + system_path, + release_folder, + static_restraints, + conformational_restraints, + ) + for window in tqdm(release_windows, disable=self._disable_progress) + ) + + @staticmethod + def _apply_attach_to_system( + window, + complex_path, + system_path, + attach_folder, + static_restraints, + conformational_restraints, + guest_restraints, + symmetry_restraints, + wall_restraints, + ): + """Function that applies `attach` restraints that is to be wrapped in `delayed` for parallelism.""" + + window_number, phase = parse_window(window) + folder = f"{attach_folder}/{window}" + os.makedirs(folder, exist_ok=True) + + with open(system_path, "r") as file: + system = openmm.XmlSerializer.deserialize(file.read()) + + for restraint in static_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=10) + for restraint in conformational_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=11) + for restraint in guest_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=12) + for restraint in symmetry_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=13) + for restraint in wall_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=14) + + apply_positional_restraints(complex_path, system, force_group=15) + + new_system_path = f"{folder}/restrained.xml" + with open(new_system_path, "w") as file: + file.write(openmm.XmlSerializer.serialize(system)) + + @staticmethod + def _apply_pull_to_system( + window, + complex_path, + system_path, + pull_folder, + static_restraints, + conformational_restraints, + guest_restraints, + ): + """Function that applies `pull` restraints that is to be wrapped in `delayed` for parallelism.""" + + window_number, phase = parse_window(window) + folder = f"{pull_folder}/{window}" + os.makedirs(folder, exist_ok=True) + + with open(system_path, "r") as file: + system = openmm.XmlSerializer.deserialize(file.read()) + + for restraint in static_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=10) + for restraint in conformational_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=11) + for restraint in guest_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=12) + + apply_positional_restraints(complex_path, system, force_group=15) + + new_system_path = f"{folder}/restrained.xml" + with open(new_system_path, "w") as file: + file.write(openmm.XmlSerializer.serialize(system)) + + @staticmethod + def _apply_release_to_system( + window, + host_solvated_path, + system_path, + release_folder, + static_restraints, + conformational_restraints, + ): + """Function that applies `release` restraints that is to be wrapped in `delayed` for parallelism.""" + + window_number, phase = parse_window(window) + folder = f"{release_folder}/{window}" + os.makedirs(folder, exist_ok=True) + + with open(system_path, "r") as file: + system = openmm.XmlSerializer.deserialize(file.read()) + + for restraint in static_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=10) + for restraint in conformational_restraints: + apply_dat_restraint(system, restraint, phase, window_number, force_group=11) + + apply_positional_restraints(host_solvated_path, system, force_group=15) + + new_system_path = f"{folder}/restrained.xml" + with open(new_system_path, "w") as file: + file.write(openmm.XmlSerializer.serialize(system)) + + def extend_pull_distance(self, extend_by: unit.Quantity): + """Extend the pull distance further than what is in the original Taproom Metadata. + + Parameters + ---------- + extend_by: openff.unit.units.Quantity + The extended distance to pull the guest molecule to. + """ + + # Determine dr between windows + pull_distance = ( + self._restraints["guest_restraints"][0]["pull"]["target"] + - self._restraints["guest_restraints"][0]["attach"]["target"] + ).m_as(unit.angstrom) + n_pull_windows = self._host_yaml_schema["calculation"]["windows"]["pull"] + dr = pull_distance / n_pull_windows + + # Update distance + self._restraints["guest_restraints"][0]["pull"]["target"] = extend_by + new_pull_distance = ( + self._restraints["guest_restraints"][0]["pull"]["target"] + - self._restraints["guest_restraints"][0]["attach"]["target"] + ).m_as(unit.angstrom) + updated_n_windows = int(new_pull_distance / dr) + + self._host_yaml_schema["calculation"]["windows"]["pull"] = updated_n_windows + + def build_system(self, n_cpus: int = 1, clean_files: bool = False): + """Build the APR files in the order: + + (1) Generate Structures and add dummy atoms + (2) Apply `attach` restraints + (3) Apply `pull` restraints + (4) Apply `release` restraints + + Parameters + ---------- + n_cpus: int + Number of CPUs to spread the workload. + clean_files: bool + Option to delete temporary files. + """ + + self._n_cpus = n_cpus + + if not self._build_array: + self._build_apr_structures() + self._apply_attach_restraints() + self._apply_pull_restraints() + self._apply_release_restraints() + else: + self._build_apr_structures() + + if clean_files: + shutil.rmtree(self._build_folder) diff --git a/paprika/evaluator/amber.py b/paprika/evaluator/amber.py index e2e4ef23..a633948c 100644 --- a/paprika/evaluator/amber.py +++ b/paprika/evaluator/amber.py @@ -7,11 +7,11 @@ import subprocess from typing import Optional -import numpy as np -import parmed as pmd +import numpy +import parmed logger = logging.getLogger(__name__) -_PI_ = np.pi +_PI_ = numpy.pi def generate_gaff( @@ -59,7 +59,7 @@ def generate_gaff( "Checking to see if we have a multi-residue MOL2 file that should be converted " "to single-residue..." ) - structure = pmd.load_file( + structure = parmed.load_file( os.path.join(directory_path, f"{output_name}.{gaff_version}.mol2"), structure=True, ) @@ -129,7 +129,6 @@ def _generate_gaff_atom_types( cwd=directory_path, ) p.communicate() - print(p) remove_files = [ "ANTECHAMBER_AC.AC", diff --git a/paprika/evaluator/analyze.py b/paprika/evaluator/analyze.py index 95371359..7dd879f6 100644 --- a/paprika/evaluator/analyze.py +++ b/paprika/evaluator/analyze.py @@ -3,12 +3,14 @@ """ import logging -from typing import Any, Dict, List +from typing import Any, Dict, List, Union -import numpy as np +import numpy +from openff.units import unit as openff_unit from paprika.analysis import fe_calc from paprika.restraints import DAT_restraint +from paprika.utils import check_unit logger = logging.getLogger(__name__) @@ -75,10 +77,10 @@ def compute_phase_free_energy( @classmethod def compute_ref_state_work( cls, - temperature: float, + temperature: Union[float, openff_unit.Quantity], guest_restraints: List[DAT_restraint], - ) -> float: - """Computes the reference state work of the attach phase. + ) -> openff_unit.Quantity: + """Computes the reference state work of the 'attach' phase. Parameters ---------- @@ -95,50 +97,49 @@ def compute_ref_state_work( analysis = fe_calc() analysis.temperature = temperature - distance_restraint = next( - ( - restraint - for restraint in guest_restraints - if "DM" in restraint.mask1 - and restraint.mask2 is not None - and restraint.mask3 is None - and restraint.mask4 is None - ), - None, - ) - - theta_restraint = next( - ( - restraint - for restraint in guest_restraints - if "DM" in restraint.mask1 - and "DM" in restraint.mask2 - and restraint.mask3 is not None - and restraint.mask4 is None - ), - None, - ) - beta_restraint = next( - ( - restraint - for restraint in guest_restraints - if "DM" in restraint.mask1 - and "DM" not in restraint.mask2 - and restraint.mask3 is not None - and restraint.mask4 is None - ), - None, - ) - - if not distance_restraint or not theta_restraint or not beta_restraint: + boresch_restraints = { + "r": None, + "theta": None, + "phi": None, + "alpha": None, + "beta": None, + "gamma": None, + } + + for restraint in guest_restraints: + # Distance + if not restraint.mask3 and not restraint.mask4: + boresch_restraints["r"] = restraint + + # Angle + elif restraint.mask3 and not restraint.mask4: + if "DM2" in restraint.mask1 and "DM1" in restraint.mask2: + boresch_restraints["theta"] = restraint + elif "DM1" in restraint.mask1 and "DM" not in restraint.mask2: + boresch_restraints["beta"] = restraint + + # Dihedral restraints + if restraint.mask4: + if ( + "DM3" in restraint.mask1 + and "DM2" in restraint.mask2 + and "DM1" in restraint.mask3 + ): + boresch_restraints["phi"] = restraint + + elif "DM2" in restraint.mask1 and "DM1" in restraint.mask2: + boresch_restraints["beta"] = restraint + + elif "DM1" in restraint.mask1: + boresch_restraints["gamma"] = restraint + + if boresch_restraints["r"] is None: raise RuntimeError( - "Could not determine the r, θ, or β restraint for computing the " - "analytic release step." + "Could not determine the `r` restraint! Need at least `r` restraint " + "for computing the analytic release step." ) - analysis.compute_ref_state_work( - [distance_restraint, theta_restraint, None, None, beta_restraint, None] - ) + analysis.compute_ref_state_work(boresch_restraints) return analysis.results["ref_state_work"] @@ -146,19 +147,19 @@ def compute_ref_state_work( def symmetry_correction( cls, n_microstates: int, - temperature: float, - ) -> float: + temperature: Union[float, openff_unit.Quantity], + ) -> openff_unit.Quantity: """Computes the free energy corrections to apply to symmetrical guest when the guest is restrained to just one of several possible symmetrical configurations (e.g. butane restrained in a cyclic host). Parameters ---------- - temperature - The temperature at which the calculation was performed in units of kelvin. n_microstates The number of different symmetrical microstates that the guest can exist in (e.g. for butane this is two). + temperature + The temperature at which the calculation was performed in units of kelvin. Returns ------- @@ -168,6 +169,9 @@ def symmetry_correction( assert n_microstates > 0 if n_microstates == 1: - return 0.0 + return 0.0 * openff_unit.kcal / openff_unit.mol + + k_B = 1.987204118e-3 * openff_unit.kcal / openff_unit.mol / openff_unit.kelvin + temperature = check_unit(temperature, base_unit=openff_unit.kelvin) - return -temperature * 0.001987204258640832 * np.log(n_microstates) + return -temperature * k_B * numpy.log(n_microstates) diff --git a/paprika/evaluator/setup.py b/paprika/evaluator/setup.py index cfcbd30a..10a4a8e0 100644 --- a/paprika/evaluator/setup.py +++ b/paprika/evaluator/setup.py @@ -5,15 +5,14 @@ import logging from typing import Any, Dict, List, Optional -import numpy as np -import parmed as pmd -from openff.units import unit as openff_unit +import numpy +import parmed from paprika.build import align from paprika.restraints import DAT_restraint, static_DAT_restraint logger = logging.getLogger(__name__) -_PI_ = np.pi +_PI_ = numpy.pi class Setup: @@ -25,7 +24,7 @@ class Setup: @classmethod def prepare_host_structure( cls, coordinate_path: str, host_atom_indices: Optional[List[int]] = None - ) -> pmd.Structure: + ) -> parmed.Structure: """Prepares the coordinates of a host molecule ready for the release phase of an APR calculation. This currently involves aligning the cavity of the host along the z-axis, and @@ -53,7 +52,7 @@ def prepare_host_structure( """ # noinspection PyTypeChecker - structure = pmd.load_file(coordinate_path, structure=True) + structure = parmed.load_file(coordinate_path, structure=True) # Extract the host from the full structure. if not host_atom_indices: @@ -64,8 +63,9 @@ def prepare_host_structure( ] # noinspection PyTypeChecker - center_of_mass: np.ndarray = pmd.geometry.center_of_mass( - host_structure.coordinates, masses=np.ones(len(host_structure.coordinates)) + center_of_mass: numpy.ndarray = parmed.geometry.center_of_mass( + host_structure.coordinates, + masses=numpy.ones(len(host_structure.coordinates)), ) # Remove the COM from the host coordinates to make alignment easier. @@ -75,19 +75,21 @@ def prepare_host_structure( # Find the principal components of the host, take the two largest, and find # the vector orthogonal to that. Use that vector to align with the z-axis. # This may not generalize to non-radially-symmetric host molecules. - inertia_tensor = np.dot( + inertia_tensor = numpy.dot( host_structure.coordinates.transpose(), host_structure.coordinates ) - eigenvalues, eigenvectors = np.linalg.eig(inertia_tensor) - order = np.argsort(eigenvalues) + eigenvalues, eigenvectors = numpy.linalg.eig(inertia_tensor) + order = numpy.argsort(eigenvalues) _, axis_2, axis_1 = eigenvectors[:, order].transpose() - cavity_axis = np.cross(axis_1, axis_2) + cavity_axis = numpy.cross(axis_1, axis_2) # Add dummy atoms which will be used to align the structure. - cls.add_dummy_atoms_to_structure(structure, [np.array([0, 0, 0]), cavity_axis]) + cls.add_dummy_atoms_to_structure( + structure, [numpy.array([0, 0, 0]), cavity_axis] + ) # Give atoms uniform mass so that the align code uses the center # of geometry rather than the center of mass. @@ -100,7 +102,7 @@ def prepare_host_structure( # have not been changed and dummy atoms not added. # noinspection PyTypeChecker - structure: pmd.Structure = pmd.load_file(coordinate_path, structure=True) + structure: parmed.Structure = parmed.load_file(coordinate_path, structure=True) structure.coordinates = aligned_structure["!:DM1&!:DM2"].coordinates return structure @@ -114,7 +116,7 @@ def prepare_complex_structure( pull_distance: float, pull_window_index: int, n_pull_windows: int, - ) -> pmd.Structure: + ) -> parmed.Structure: """Prepares the coordinates of a host molecule ready for the pull (+ attach) phase of an APR calculation. @@ -158,7 +160,7 @@ def prepare_complex_structure( # Align the host-guest complex so the first guest atom is at (0, 0, 0) and the # second guest atom lies along the positive z-axis. # noinspection PyTypeChecker - structure: pmd.Structure = pmd.load_file(coordinate_path, structure=True) + structure: parmed.Structure = parmed.load_file(coordinate_path, structure=True) ( guest_orientation_mask_0, @@ -169,7 +171,7 @@ def prepare_complex_structure( structure, guest_orientation_mask_0, guest_orientation_mask_1 ) - target_distance = np.linspace(0.0, pull_distance, n_pull_windows)[ + target_distance = numpy.linspace(0.0, pull_distance, n_pull_windows)[ pull_window_index ] target_difference = target_distance @@ -181,9 +183,9 @@ def prepare_complex_structure( @staticmethod def add_dummy_atoms_to_structure( - structure: pmd.Structure, - dummy_atom_offsets: List[np.ndarray], - offset_coordinates: Optional[np.ndarray] = None, + structure: parmed.Structure, + dummy_atom_offsets: List[numpy.ndarray], + offset_coordinates: Optional[numpy.ndarray] = None, ): """A convenience method to add a number of dummy atoms to an existing ParmEd structure, and to position those atoms at a specified set of positions. @@ -201,9 +203,9 @@ def add_dummy_atoms_to_structure( """ if offset_coordinates is None: - offset_coordinates = np.zeros(3) + offset_coordinates = numpy.zeros(3) - full_coordinates = np.vstack( + full_coordinates = numpy.vstack( [ structure.coordinates, *[ @@ -214,7 +216,7 @@ def add_dummy_atoms_to_structure( ) for index in range(len(dummy_atom_offsets)): - structure.add_atom(pmd.Atom(name="DUM"), f"DM{index + 1}", 1) + structure.add_atom(parmed.Atom(name="DUM"), f"DM{index + 1}", 1) structure.positions = full_coordinates @@ -314,7 +316,7 @@ def build_conformational_restraints( The path to the coordinate file which the restraints will be applied to. This should contain either the host or the complex, the dummy atoms and solvent. attach_lambdas - The values 'lambda' being used during the attach phase of the APR + The values 'lambda' being used during the 'attach' phase of the APR calculation. n_pull_windows The total number of pull windows being used in the APR calculation. @@ -389,7 +391,7 @@ def build_symmetry_restraints( * a ``force_constant`` entry which specifies the force constant of the restraint. - These 'schemas` map directly to the 'restraints -> symmetry_correction + These `schemas` map directly to the 'restraints -> symmetry_correction -> restraint' dictionaries specified in the `taproom` guest YAML files. Parameters @@ -397,7 +399,7 @@ def build_symmetry_restraints( coordinate_path The path to the coordinate file which the restraints will be applied to. This should contain either the host or the complex, the dummy atoms and - and solvent. + solvent. n_attach_windows The total number of attach windows being used in the APR calculation. restraint_schemas @@ -426,16 +428,24 @@ def build_symmetry_restraints( restraint.attach["fc_final"] = restraint_schema["force_constant"] restraint.attach["fraction_list"] = [1.0] * n_attach_windows + restraint.attach["target"] = restraint_schema["target"] - # This target will be overridden by the custom values. - restraint.attach["target"] = 91 * openff_unit.degrees - restraint.custom_restraint_values["r2"] = 91 * openff_unit.degrees - restraint.custom_restraint_values["r3"] = 91 * openff_unit.degrees - - # 0 force constant between 91 degrees and 180 degrees. - restraint.custom_restraint_values["rk3"] = ( - 0.0 * openff_unit.kcal / openff_unit.mole / openff_unit.radians**2 + # Set upper bounds to zero + restraint.custom_restraint_values["r3"] = ( + 0 * restraint.attach["target"].units ) + restraint.custom_restraint_values["r4"] = ( + 0 * restraint.attach["target"].units + ) + + # Harmonic force constant beyond target distance. + restraint.custom_restraint_values["rk2"] = restraint_schema[ + "force_constant" + ] + restraint.custom_restraint_values["rk3"] = restraint_schema[ + "force_constant" + ] + restraint.initialize() restraints.append(restraint) @@ -461,7 +471,7 @@ def build_wall_restraints( restraint. * a ``target`` entry which specifies the target value of the restraint. - These 'schemas` map directly to the 'restraints -> wall_restraints -> restraint' + These `schemas` map directly to the 'restraints -> wall_restraints -> restraint' dictionaries specified in the `taproom` guest YAML files. Parameters @@ -469,7 +479,7 @@ def build_wall_restraints( coordinate_path The path to the coordinate file which the restraints will be applied to. This should contain either the host or the complex, the dummy atoms and - and solvent. + solvent. n_attach_windows The total number of attach windows being used in the APR calculation. restraint_schemas @@ -487,21 +497,29 @@ def build_wall_restraints( restraints = [] for restraint_schema in restraint_schemas: + mask = restraint_schema["atoms"].split() + restraint = DAT_restraint() restraint.auto_apr = True restraint.continuous_apr = False restraint.amber_index = use_amber_indices restraint.topology = coordinate_path - restraint.mask1 = restraint_schema["atoms"].split()[0] - restraint.mask2 = restraint_schema["atoms"].split()[1] + restraint.mask1 = mask[0] + restraint.mask2 = mask[1] + restraint.mask3 = mask[2] if len(mask) > 2 else None + restraint.mask4 = mask[3] if len(mask) > 3 else None restraint.attach["fc_final"] = restraint_schema["force_constant"] restraint.attach["fraction_list"] = [1.0] * n_attach_windows restraint.attach["target"] = restraint_schema["target"] - # Minimum distance is 0 Angstrom - restraint.custom_restraint_values["r1"] = 0 * openff_unit.degrees - restraint.custom_restraint_values["r2"] = 0 * openff_unit.degrees + # Set lower bounds to zero + restraint.custom_restraint_values["r1"] = ( + 0 * restraint.attach["target"].units + ) + restraint.custom_restraint_values["r2"] = ( + 0 * restraint.attach["target"].units + ) # Harmonic force constant beyond target distance. restraint.custom_restraint_values["rk2"] = restraint_schema[ @@ -539,7 +557,7 @@ def build_guest_restraints( restraint. * a ``target`` entry which specifies the target value of the restraint. - These 'schemas` map directly to the 'restraints -> guest -> restraint' + These `schemas` map directly to the 'restraints -> guest -> restraint' dictionaries specified in the `taproom` guest YAML files. Parameters @@ -547,9 +565,9 @@ def build_guest_restraints( coordinate_path The path to the coordinate file which the restraints will be applied to. This should contain either the host or the complex, the dummy atoms and - and solvent. + solvent. attach_lambdas - The values 'lambda' being used during the attach phase of the APR + The values 'lambda' being used during the 'attach' phase of the APR calculation. n_pull_windows The total number of pull windows being used in the APR calculation. diff --git a/paprika/evaluator/utils.py b/paprika/evaluator/utils.py deleted file mode 100644 index a602ab2b..00000000 --- a/paprika/evaluator/utils.py +++ /dev/null @@ -1,13 +0,0 @@ -import pkg_resources - - -def get_benchmarks(): - """ - Determine the installed ``taproom`` benchmarks. - """ - installed_benchmarks = {} - - for entry_point in pkg_resources.iter_entry_points(group="taproom.benchmarks"): - installed_benchmarks[entry_point.name] = entry_point.load() - - return installed_benchmarks diff --git a/paprika/io.py b/paprika/io.py index d5eaba0f..47a83b85 100644 --- a/paprika/io.py +++ b/paprika/io.py @@ -6,13 +6,13 @@ import traceback from enum import Enum -import numpy as np -import pytraj as pt +import numpy +import pytraj from openff.units import unit as openff_unit from parmed import Structure from parmed.amber import AmberParm -from paprika.restraints import BiasPotentialType, DAT_restraint, RestraintType +from paprika.restraints import DAT_restraint # https://stackoverflow.com/questions/27909658/json-encoder-and-decoder-for-complex-numpy-arrays # https://stackoverflow.com/a/24375113/901925 @@ -80,11 +80,11 @@ def default(self, obj): logging.warning("Encountered Structure, which does not store filename.") return "" - if isinstance(obj, np.ndarray): + if isinstance(obj, numpy.ndarray): if obj.flags["C_CONTIGUOUS"]: obj_data = obj.data else: - cont_obj = np.ascontiguousarray(obj) + cont_obj = numpy.ascontiguousarray(obj) assert cont_obj.flags["C_CONTIGUOUS"] obj_data = cont_obj.data data_b64 = base64.b64encode(obj_data) @@ -97,23 +97,25 @@ def default(self, obj): elif isinstance( obj, ( - np.int_, - np.intc, - np.intp, - np.int8, - np.int16, - np.int32, - np.int64, - np.uint8, - np.uint16, - np.uint32, - np.uint64, + numpy.int_, + numpy.intc, + numpy.intp, + numpy.int8, + numpy.int16, + numpy.int32, + numpy.int64, + numpy.uint8, + numpy.uint16, + numpy.uint32, + numpy.uint64, ), ): return int(obj) - elif isinstance(obj, (np.float_, np.float16, np.float32, np.float64)): + elif isinstance( + obj, (numpy.float_, numpy.float16, numpy.float32, numpy.float64) + ): return float(obj) - elif isinstance(obj, (np.ndarray,)): + elif isinstance(obj, (numpy.ndarray,)): return obj.tolist() elif isinstance(obj, openff_unit.Quantity): return serialize_quantity(obj) @@ -139,7 +141,7 @@ def __init__(self, *args, **kwargs): def custom_object_hook(self, obj): if "__ndarray__" in obj: data = base64.b64decode(obj["__ndarray__"]) - return np.frombuffer(data, obj["dtype"]).reshape(obj["shape"]) + return numpy.frombuffer(data, obj["dtype"]).reshape(obj["shape"]) if "@type" in obj: if obj["@type"] == "openff.units.unit.Quantity": @@ -343,7 +345,7 @@ def load_trajectory(window, trajectory, topology, single_topology=False): ) logger.debug(f"Loading {os.path.join(window, topology)} and {trajectory_path}") try: - traj = pt.iterload(trajectory_path, os.path.join(window, topology)) + traj = pytraj.iterload(trajectory_path, os.path.join(window, topology)) except ValueError as e: formatted_exception = traceback.format_exception(None, e, e.__traceback__) logger.info( @@ -351,10 +353,10 @@ def load_trajectory(window, trajectory, topology, single_topology=False): f"{formatted_exception}" ) elif isinstance(topology, str) and single_topology: - traj = pt.iterload(trajectory_path, os.path.join(topology)) + traj = pytraj.iterload(trajectory_path, os.path.join(topology)) else: try: - traj = pt.iterload(trajectory_path, topology) + traj = pytraj.iterload(trajectory_path, topology) except BaseException: raise Exception("Tried to load `topology` object directly and failed.") @@ -384,7 +386,7 @@ def read_restraint_data( Returns ------- - data: :class:`np.array` + data: :class:`numpy.array` The values for this restraint in this window """ @@ -397,7 +399,7 @@ def read_restraint_data( and not restraint.mask4 ): data = openff_unit.Quantity( - pt.distance( + pytraj.distance( trajectory, " ".join([restraint.mask1, restraint.mask2]), image=True ), units=openff_unit.angstrom, @@ -407,7 +409,7 @@ def read_restraint_data( restraint.mask1 and restraint.mask2 and restraint.mask3 and not restraint.mask4 ): data = openff_unit.Quantity( - pt.angle( + pytraj.angle( trajectory, " ".join([restraint.mask1, restraint.mask2, restraint.mask3]), ), @@ -416,7 +418,7 @@ def read_restraint_data( elif restraint.mask1 and restraint.mask2 and restraint.mask3 and restraint.mask4: data = openff_unit.Quantity( - pt.dihedral( + pytraj.dihedral( trajectory, " ".join( [restraint.mask1, restraint.mask2, restraint.mask3, restraint.mask4] diff --git a/paprika/restraints/__init__.py b/paprika/restraints/__init__.py index 40c473b7..645264e0 100644 --- a/paprika/restraints/__init__.py +++ b/paprika/restraints/__init__.py @@ -7,7 +7,7 @@ RestraintType, static_DAT_restraint, ) -from .utils import create_window_list +from .utils import create_window_list, parse_window __all__ = [ "BiasPotentialType", @@ -18,4 +18,5 @@ "Plumed", "Colvars", "create_window_list", + "parse_window", ] diff --git a/paprika/restraints/colvars.py b/paprika/restraints/colvars.py index 5149d904..dcdf5ced 100644 --- a/paprika/restraints/colvars.py +++ b/paprika/restraints/colvars.py @@ -1,7 +1,7 @@ import logging import os -import numpy as np +import numpy from openff.units import unit as openff_unit from paprika.restraints.plumed import Plumed @@ -10,7 +10,7 @@ logger = logging.getLogger(__name__) -_PI_ = np.pi +_PI_ = numpy.pi class Colvars(Plumed): diff --git a/paprika/restraints/openmm.py b/paprika/restraints/openmm.py index 84684900..415534b7 100644 --- a/paprika/restraints/openmm.py +++ b/paprika/restraints/openmm.py @@ -1,7 +1,7 @@ """A module aimed at applying restraints directly to OpenMM systems.""" import logging -import numpy as np +import numpy try: import openmm @@ -12,14 +12,14 @@ from typing import Optional, Union -import parmed as pmd +import parmed from openff.units import unit as openff_unit from openff.units.openmm import to_openmm from paprika.restraints import DAT_restraint logger = logging.getLogger(__name__) -_PI_ = np.pi +_PI_ = numpy.pi def apply_positional_restraints( @@ -50,7 +50,7 @@ def apply_positional_restraints( """ # noinspection PyTypeChecker - structure: pmd.Structure = pmd.load_file(coordinate_path, structure=True) + structure: parmed.Structure = parmed.load_file(coordinate_path, structure=True) for atom in structure.atoms: if atom.name == atom_name: diff --git a/paprika/restraints/plumed.py b/paprika/restraints/plumed.py index 5943ea62..32dfed43 100644 --- a/paprika/restraints/plumed.py +++ b/paprika/restraints/plumed.py @@ -1,7 +1,7 @@ import logging import os -import numpy as np +import numpy from openff.units import unit as openff_unit from parmed.structure import Structure as ParmedStructureClass @@ -11,7 +11,7 @@ logger = logging.getLogger(__name__) -_PI_ = np.pi +_PI_ = numpy.pi _plumed_unit_dict = { openff_unit.kcal / openff_unit.mole: "kcal/mol", diff --git a/paprika/restraints/read_yaml.py b/paprika/restraints/read_yaml.py deleted file mode 100644 index 59c34dd0..00000000 --- a/paprika/restraints/read_yaml.py +++ /dev/null @@ -1,75 +0,0 @@ -import logging -import re - -import yaml - -logger = logging.getLogger(__name__) - - -def read_yaml(file): - """ - Read `Taproom `_ -style YAML-formatted instructions for - preparing host-guest systems. - - Parameters - ---------- - file: os.PathLike - A YAML-formatted file. - - Returns - ------- - yaml_data: dict - Dictionary containing simulation setup parameters. - - """ - - with open(file, "r") as f: - yaml_data = yaml.safe_load(f) - logger.debug(yaml_data) - - if "aliases" in yaml_data.keys(): - logger.debug("Dealiasing atom masks...") - yaml_data = de_alias(yaml_data) - - return yaml_data - - -def multiple_replace(dct, text): - """ - Create a regular expression to do multiple find and replace. - """ - - # Create a regular expression from the dictionary keys - regex = re.compile("(%s)" % "|".join(map(re.escape, dct.keys()))) - - # For each match, look-up corresponding value in dictionary - return regex.sub(lambda mo: dct[mo.string[mo.start() : mo.end()]], text) - - -def de_alias(yaml_data): - """ - Replace aliased atoms in a ``taproom`` recipe. - """ - - mapping_list = yaml_data["aliases"] - mapping_dictionary = {} - - for atom_pair in mapping_list: - mapping_dictionary.update(atom_pair) - logger.debug(f"Found mapping: {mapping_dictionary}") - - for restraint_type, restraint_type_list in yaml_data["restraints"].items(): - for restraint in restraint_type_list: - atoms = restraint["restraint"]["atoms"] - mapped_atoms = multiple_replace(mapping_dictionary, atoms) - logger.info(f"{atoms} → {mapped_atoms}") - restraint["restraint"]["atoms"] = mapped_atoms - - if "symmetry_correction" in yaml_data.keys(): - for restraint in yaml_data["symmetry_correction"]["restraints"]: - atoms = restraint["restraint"]["atoms"] - mapped_atoms = multiple_replace(mapping_dictionary, atoms) - logger.info(f"{atoms} → {mapped_atoms}") - restraint["restraint"]["atoms"] = mapped_atoms - - return yaml_data diff --git a/paprika/restraints/utils.py b/paprika/restraints/utils.py index 89380b7e..158eda6e 100644 --- a/paprika/restraints/utils.py +++ b/paprika/restraints/utils.py @@ -1,6 +1,6 @@ import logging -import numpy as np +import numpy from openff.units import unit as openff_unit from paprika.restraints.restraints import ( @@ -12,7 +12,7 @@ logger = logging.getLogger(__name__) -_PI_ = np.pi +_PI_ = numpy.pi def parse_window(window): diff --git a/paprika/simulate/amber.py b/paprika/simulate/amber.py index 08ac89f1..6faed277 100644 --- a/paprika/simulate/amber.py +++ b/paprika/simulate/amber.py @@ -1,7 +1,7 @@ import abc import logging import os -import subprocess as sp +import subprocess from collections import OrderedDict from enum import Enum @@ -633,11 +633,11 @@ def run(self, soft_minimize=False, overwrite=False, fail_ok=False): logger.debug("Exec line: " + " ".join(exec_list)) # Execute - amber_output = sp.Popen( + amber_output = subprocess.Popen( exec_list, cwd=self.path, - stdout=sp.PIPE, - stderr=sp.PIPE, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, env=os.environ, ) diff --git a/paprika/simulate/gromacs.py b/paprika/simulate/gromacs.py index 09ca00b1..2f665754 100644 --- a/paprika/simulate/gromacs.py +++ b/paprika/simulate/gromacs.py @@ -2,7 +2,7 @@ import glob import logging import os -import subprocess as sp +import subprocess from collections import OrderedDict from enum import Enum @@ -623,11 +623,11 @@ def run(self, run_grompp=True, overwrite=False, fail_ok=False): grompp_list += ["-n", self.index_file] # Run GROMPP - grompp_output = sp.Popen( + grompp_output = subprocess.Popen( grompp_list, cwd=self.path, - stdout=sp.PIPE, - stderr=sp.PIPE, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, env=os.environ, ) grompp_stdout = grompp_output.stdout.read().splitlines() @@ -707,11 +707,11 @@ def run(self, run_grompp=True, overwrite=False, fail_ok=False): mdrun_list += ["-plumed", self.plumed_file] # Run MDRUN - mdrun_output = sp.Popen( + mdrun_output = subprocess.Popen( mdrun_list, cwd=self.path, - stdout=sp.PIPE, - stderr=sp.PIPE, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, env=os.environ, ) mdrun_out = mdrun_output.stdout.read().splitlines() diff --git a/paprika/simulate/namd.py b/paprika/simulate/namd.py index 2f832c8d..03278984 100644 --- a/paprika/simulate/namd.py +++ b/paprika/simulate/namd.py @@ -1,12 +1,12 @@ import abc import logging import os -import subprocess as sp +import subprocess from collections import OrderedDict from enum import Enum -import numpy as np -import parmed as pmd +import numpy +import parmed from paprika.utils import get_dict_without_keys @@ -584,14 +584,14 @@ def _get_cell_basis_vectors(self): """ Function to calculate the PBC cell basis vectors (needed when running a simulation for the first time). """ - structure = pmd.load_file( + structure = parmed.load_file( os.path.join(self.path, self.topology), os.path.join(self.path, self.coordinates), structure=True, ) coordinates = structure.coordinates - masses = np.ones(len(coordinates)) - center = pmd.geometry.center_of_mass(coordinates, masses) + masses = numpy.ones(len(coordinates)) + center = parmed.geometry.center_of_mass(coordinates, masses) self.cell_basis_vectors["cellOrigin"] = list(center) @@ -878,11 +878,11 @@ def run(self, overwrite=True, fail_ok=False): logger.debug("Exec line: " + " ".join(exec_list)) # Execute - namd_output = sp.Popen( + namd_output = subprocess.Popen( exec_list, cwd=self.path, stdout=open(os.path.join(self.path, self.logfile), "w"), - stderr=sp.PIPE, + stderr=subprocess.PIPE, env=os.environ, ) namd_stderr = namd_output.stderr.read().splitlines() diff --git a/paprika/taproom/__init__.py b/paprika/taproom/__init__.py new file mode 100644 index 00000000..2c3e1e4e --- /dev/null +++ b/paprika/taproom/__init__.py @@ -0,0 +1,6 @@ +from .taproom import get_benchmarks, read_yaml_schema + +__all__ = [ + "get_benchmarks", + "read_yaml_schema", +] diff --git a/paprika/taproom/taproom.py b/paprika/taproom/taproom.py new file mode 100644 index 00000000..6a520457 --- /dev/null +++ b/paprika/taproom/taproom.py @@ -0,0 +1,60 @@ +import logging +import sys + +import yaml + +from paprika.taproom.utils import convert_string_to_quantity, de_alias + +if sys.version_info.minor < 10: + from pkg_resources import iter_entry_points as entry_points +else: + from importlib.metadata import entry_points + + +logger = logging.getLogger(__name__) + + +def get_benchmarks(): + """ + Determine the installed ``taproom`` benchmarks. + """ + installed_benchmarks = {} + + for entry_point in entry_points(group="taproom.benchmarks"): + installed_benchmarks[entry_point.name] = entry_point.load() + + return installed_benchmarks + + +def read_yaml_schema(file): + """ + Read `Taproom `_ -style YAML-formatted instructions for + preparing host-guest systems. + + Parameters + ---------- + file: os.PathLike + A YAML-formatted file. + + Returns + ------- + yaml_data: dict + Dictionary containing simulation setup parameters. + + """ + + # Read YAML file + with open(file, "r") as f: + yaml_data = yaml.safe_load(f) + logger.debug(yaml_data) + + # Convert aliases to atom masks + if "aliases" in yaml_data.keys(): + logger.debug("Dealiasing atom masks...") + yaml_data = de_alias(yaml_data) + + # Convert all string to OpenFF Quantity + logger.debug("Converting string to unit.Quantity...") + convert_string_to_quantity(yaml_data) + + return yaml_data diff --git a/paprika/taproom/utils.py b/paprika/taproom/utils.py new file mode 100644 index 00000000..229e1ece --- /dev/null +++ b/paprika/taproom/utils.py @@ -0,0 +1,63 @@ +import logging + +from openff.units import unit as openff_unit + +from paprika.utils import multiple_replace + +logger = logging.getLogger(__name__) + + +def de_alias(yaml_data): + """ + Replace aliased atoms in a ``taproom`` recipe. + """ + + mapping_list = yaml_data["aliases"] + mapping_dictionary = {} + + for atom_pair in mapping_list: + mapping_dictionary.update(atom_pair) + logger.debug(f"Found mapping: {mapping_dictionary}") + + for restraint_type, restraint_type_list in yaml_data["restraints"].items(): + for restraint in restraint_type_list: + atoms = restraint["restraint"]["atoms"] + mapped_atoms = multiple_replace(mapping_dictionary, atoms) + logger.info(f"{atoms} → {mapped_atoms}") + restraint["restraint"]["atoms"] = mapped_atoms + + if "symmetry_correction" in yaml_data.keys(): + for restraint in yaml_data["symmetry_correction"]["restraints"]: + atoms = restraint["restraint"]["atoms"] + mapped_atoms = multiple_replace(mapping_dictionary, atoms) + logger.info(f"{atoms} → {mapped_atoms}") + restraint["restraint"]["atoms"] = mapped_atoms + + return yaml_data + + +def convert_string_to_quantity(yaml_data): + """ + Convert strings for 'force_constant' and 'targets' to `unit.Quantity` + """ + + def _to_quantity(string, key): + if string in key: + value = key[string] + key[string] = openff_unit.Quantity(value) + + for restraint_type, restraint_type_list in yaml_data["restraints"].items(): + for restraint in restraint_type_list: + _to_quantity("force_constant", restraint["restraint"]) + _to_quantity("target", restraint["restraint"]) + if "attach" in restraint["restraint"]: + _to_quantity("force_constant", restraint["restraint"]["attach"]) + _to_quantity("target", restraint["restraint"]["attach"]) + if "pull" in restraint["restraint"]: + _to_quantity("force_constant", restraint["restraint"]["pull"]) + _to_quantity("target", restraint["restraint"]["pull"]) + + if "symmetry_correction" in yaml_data.keys(): + for restraint in yaml_data["symmetry_correction"]["restraints"]: + _to_quantity("force_constant", restraint["restraint"]) + _to_quantity("target", restraint["restraint"]) diff --git a/paprika/tests/test_analysis.py b/paprika/tests/test_analysis.py index b0a30e71..5787605d 100644 --- a/paprika/tests/test_analysis.py +++ b/paprika/tests/test_analysis.py @@ -1,10 +1,12 @@ import logging import os import shutil +from copy import deepcopy -import numpy as np -import parmed as pmd +import numpy +import parmed import pytest +from openff.units import unit as openff_unit from pytest import approx from paprika import analysis, log, restraints @@ -15,6 +17,9 @@ logger = logging.getLogger(__name__) +random_seed = 12345 + + @pytest.fixture(scope="module", autouse=True) def clean_files(directory="tmp"): # This happens before the test function call @@ -28,7 +33,7 @@ def clean_files(directory="tmp"): @pytest.fixture(scope="module", autouse=True) def setup_free_energy_calculation(): - input_pdb = pmd.load_file( + input_pdb = parmed.load_file( os.path.join(os.path.dirname(__file__), "../data/cb6-but/vac.pdb") ) @@ -83,8 +88,7 @@ def setup_free_energy_calculation(): rest3.pull["num_windows"] = 19 rest3.initialize() - seed = 12345 - + # Create Analysis Instance fecalc = analysis.fe_calc() fecalc.topology = os.path.join( os.path.dirname(__file__), "../data/cb6-but-apr/vac.prmtop" @@ -92,74 +96,123 @@ def setup_free_energy_calculation(): fecalc.trajectory = "*.nc" fecalc.path = os.path.join(os.path.dirname(__file__), "../data/cb6-but-apr/") fecalc.restraint_list = [rest1, rest2, rest3] - fecalc.methods = ["ti-block", "mbar-block", "mbar-autoc", "mbar-boot"] fecalc.boot_cycles = 100 fecalc.ti_matrix = "diagonal" fecalc.compute_largest_neighbor = True fecalc.compute_roi = True + fecalc.conservative_subsample = False + fecalc.exact_sem_each_ti_fraction = False + fecalc.fractions = [0.2, 0.4, 0.6, 0.8, 1.0] + fecalc.energy_unit = openff_unit.kcal / openff_unit.mole + fecalc.distance_unit = openff_unit.angstrom + fecalc.angle_unit = openff_unit.degrees + fecalc.temperature_unit = openff_unit.kelvin fecalc.collect_data(single_topology=True) - fecalc.compute_free_energy(seed=seed) fecalc.compute_ref_state_work([rest1, rest2, rest3, None, None, None]) return fecalc def test_setup(clean_files, setup_free_energy_calculation): - pass + assert setup_free_energy_calculation.temperature_unit == openff_unit.kelvin + assert setup_free_energy_calculation.distance_unit == openff_unit.angstrom + assert setup_free_energy_calculation.angle_unit == openff_unit.degrees + assert ( + setup_free_energy_calculation.energy_unit == openff_unit.kcal / openff_unit.mole + ) + assert setup_free_energy_calculation.fractions == [0.2, 0.4, 0.6, 0.8, 1.0] + results = deepcopy(setup_free_energy_calculation.results) + setup_free_energy_calculation.results = results + assert setup_free_energy_calculation.results == results + assert setup_free_energy_calculation.exact_sem_each_ti_fraction is False + assert setup_free_energy_calculation.conservative_subsample is False def test_mbar_block(clean_files, setup_free_energy_calculation): - results = setup_free_energy_calculation.results method = "mbar-block" - # Test mbar-block free energies and uncertainties + # Estimate FE with `mbar-block` + setup_free_energy_calculation.methods = [method] + setup_free_energy_calculation.compute_free_energy(seed=random_seed) + results = setup_free_energy_calculation.results + + # Test `mbar-block` free energies and uncertainties test_vals = [ results["attach"][method]["fe"].magnitude, results["attach"][method]["sem"].magnitude, results["pull"][method]["fe"].magnitude, results["pull"][method]["sem"].magnitude, ] - reference_values = [13.267731176, 0.16892084090, -2.1791430735, 0.93638948302] + reference_values = [13.2677, 0.1689, -2.1791, 0.9364] assert reference_values == approx(test_vals, abs=0.01) - # Test attach mbar-block largest_neighbor values + # Test attach `mbar-block` largest_neighbor values test_vals = results["attach"][method]["largest_neighbor"].magnitude - reference_values = np.array([0.0198918, 0.0451676, 0.0564517, 0.1079282, 0.1079282]) + reference_values = numpy.array([0.0199, 0.0452, 0.0565, 0.1079, 0.1079]) assert reference_values == approx(test_vals, abs=0.01) - # Test pull mbar-block largest_neighbor values + # Test pull `mbar-block` largest_neighbor values test_vals = results["pull"][method]["largest_neighbor"].magnitude - reference_values = np.array( + reference_values = numpy.array( [ - 0.2053769, - 0.2053769, - 0.1617423, - 0.1747668, - 0.5255023, - 0.5255023, - 0.1149945, - 0.1707901, - 0.2129136, - 0.2129136, - 0.1942189, - 0.1768906, - 0.1997338, - 0.1997338, - 0.2014766, - 0.2014766, - 0.1470727, - 0.1442517, - 0.1434395, + # fmt: off + 0.20538, 0.20538, 0.16174, 0.17477, 0.52550, + 0.52550, 0.11499, 0.17079, 0.21291, 0.21291, + 0.19422, 0.17689, 0.19973, 0.19973, 0.20148, + 0.20148, 0.14707, 0.14425, 0.14344 + # fmt: on ] ) assert reference_values == approx(test_vals, abs=0.01) -def test_ti_block(clean_files, setup_free_energy_calculation): +def test_mbar_autoc(clean_files, setup_free_energy_calculation): + method = "mbar-autoc" + + setup_free_energy_calculation.methods = [method] + setup_free_energy_calculation.compute_free_energy(seed=random_seed) results = setup_free_energy_calculation.results + # Test `mbar-autoc` free energies and uncertainties + test_vals = [ + results["attach"][method]["fe"].magnitude, + results["attach"][method]["sem"].magnitude, + results["pull"][method]["fe"].magnitude, + results["pull"][method]["sem"].magnitude, + ] + reference_values = [13.26773, 0.0808, -2.1791, 0.6969] + assert reference_values == approx(test_vals, abs=0.01) + + # Test attach `mbar-autoc` largest_neighbor values + test_vals = results["attach"][method]["largest_neighbor"].magnitude + reference_values = numpy.array([0.0199, 0.0259, 0.0336, 0.0383, 0.0383]) + assert reference_values == approx(test_vals, abs=0.01) + + # Test pull `mbar-autoc` largest_neighbor values + test_vals = results["pull"][method]["largest_neighbor"].magnitude + numpy.savetxt("tmp/test.txt", test_vals, fmt="%.5f") + reference_values = numpy.array( + [ + # fmt: off + 0.10274, 0.11361, 0.13074, 0.14136, 0.38928, + 0.38928, 0.11215, 0.12951, 0.13145, 0.13145, + 0.13113, 0.13113, 0.13751, 0.14186, 0.14186, + 0.12847, 0.13550, 0.13550, 0.13404 + # fmt: on + ] + ) + assert reference_values == approx(test_vals, abs=0.01) + + +def test_ti_block(clean_files, setup_free_energy_calculation): method = "ti-block" + # Estimate FE with exact sem + setup_free_energy_calculation.methods = [method] + setup_free_energy_calculation.exact_sem_each_ti_fraction = True + setup_free_energy_calculation.compute_free_energy(seed=random_seed) + results = setup_free_energy_calculation.results + # Test ti-block free energies and uncertainties test_vals = [ results["attach"][method]["fe"].magnitude, @@ -167,77 +220,178 @@ def test_ti_block(clean_files, setup_free_energy_calculation): results["pull"][method]["fe"].magnitude, results["pull"][method]["sem"].magnitude, ] - reference_values = np.array([13.35, 0.26, -1.85, 0.78]) + reference_values = numpy.array([13.31, 0.25, -1.62, 0.87]) assert reference_values == approx(test_vals, abs=0.01) # ROI only runs during TI. # Test attach ti-block largest_neighbor values test_vals = results["attach"][method]["largest_neighbor"].magnitude - reference_values = np.array([0.03, 0.07, 0.10, 0.18, 0.18]) + reference_values = numpy.array([0.03, 0.07, 0.10, 0.18, 0.18]) assert reference_values == approx(test_vals, abs=0.01) # Test pull ti-block largest_neighbor values test_vals = results["pull"][method]["largest_neighbor"].magnitude - reference_values = np.array( + reference_values = numpy.array( [ - 0.33156402, - 0.33156402, - 0.22515133, - 0.2219127, - 0.2219127, - 0.1311959, - 0.13514015, - 0.15078472, - 0.15078472, - 0.12448228, - 0.10678047, - 0.10678047, - 0.10157904, - 0.14122943, - 0.16608568, - 0.16608568, - 0.14718857, - 0.14090383, - 0.11005729, + # fmt: off + 0.33156, 0.33156, 0.21509, 0.22191, 0.22191, + 0.10746, 0.13514, 0.15078, 0.15078, 0.15518, + 0.10678, 0.10678, 0.10158, 0.14123, 0.16609, + 0.16609, 0.14719, 0.14090, 0.11006 + # fmt: on ] ) assert reference_values == approx(test_vals, abs=0.01) + # Estimate FE without exact sem + setup_free_energy_calculation.exact_sem_each_ti_fraction = False + setup_free_energy_calculation.compute_free_energy(seed=random_seed) + results = setup_free_energy_calculation.results -def test_reference_state_work(clean_files, setup_free_energy_calculation): + # Test ti-block free energies and uncertainties + test_vals = [ + results["attach"][method]["fe"].magnitude, + results["attach"][method]["sem"].magnitude, + results["pull"][method]["fe"].magnitude, + results["pull"][method]["sem"].magnitude, + ] + reference_values = numpy.array([13.31, 0.25, -1.62, 0.87]) + assert reference_values == approx(test_vals, abs=0.01) + + # ROI only runs during TI. + + # Test attach ti-block largest_neighbor values + test_vals = results["attach"][method]["largest_neighbor"].magnitude + reference_values = numpy.array([0.03, 0.07, 0.10, 0.18, 0.18]) + assert reference_values == approx(test_vals, abs=0.01) + + # Test pull ti-block largest_neighbor values + test_vals = results["pull"][method]["largest_neighbor"].magnitude + + reference_values = numpy.array( + [ + # fmt: off + 0.33156, 0.33156, 0.21509, 0.22191, 0.22191, + 0.10746, 0.13514, 0.15078, 0.15078, 0.15518, + 0.10678, 0.10678, 0.10157, 0.14122, 0.16608, + 0.16608, 0.14718, 0.14090, 0.11005 + # fmt: on + ] + ) + assert reference_values == approx(test_vals, abs=0.01) + + +def test_ti_nocor(clean_files, setup_free_energy_calculation): + method = "ti-nocor" + + # Estimate FE with exact sem + setup_free_energy_calculation.methods = [method] + setup_free_energy_calculation.exact_sem_each_ti_fraction = True + setup_free_energy_calculation.compute_free_energy(seed=random_seed) results = setup_free_energy_calculation.results - assert np.isclose(-4.34372240, results["ref_state_work"].magnitude) + # Test `ti-nocor` free energies and uncertainties + test_vals = [ + results["attach"][method]["fe"].magnitude, + results["attach"][method]["sem"].magnitude, + results["pull"][method]["fe"].magnitude, + results["pull"][method]["sem"].magnitude, + ] + reference_values = numpy.array([13.34, 0.09, -1.71, 0.56]) + assert reference_values == approx(test_vals, abs=0.01) -def test_save_and_loading(clean_files, setup_free_energy_calculation): - # Save FE results to file - setup_free_energy_calculation.save_results("tmp/results.json") - assert is_file_and_not_empty("tmp/results.json") is True - setup_free_energy_calculation.save_results("tmp/results.json", overwrite=True) - assert is_file_and_not_empty("tmp/results.json") is True + # ROI only runs during TI. - # Load results - fe_calc = analysis.fe_calc() - fe_calc.load_results("tmp/results.json") - assert np.isclose(-4.34372240, fe_calc.results["ref_state_work"].magnitude) + # Test attach `ti-nocor` largest_neighbor values + test_vals = results["attach"][method]["largest_neighbor"].magnitude + reference_values = numpy.array([0.014, 0.036, 0.055, 0.066, 0.066]) + assert reference_values == approx(test_vals, abs=0.01) - # Save Simulation data to file - setup_free_energy_calculation.save_data("tmp/simulation_data.json") - assert is_file_and_not_empty("tmp/simulation_data.json") is True - setup_free_energy_calculation.save_data("tmp/simulation_data.json", overwrite=True) - assert is_file_and_not_empty("tmp/simulation_data.json") is True + # Test pull `ti-nocor` largest_neighbor values + test_vals = results["pull"][method]["largest_neighbor"].magnitude - # Load simulation data - fe_calc.load_data("tmp/simulation_data.json") - assert fe_calc.simulation_data is not None - assert fe_calc.changing_restraints is not None - assert fe_calc.orders is not None + reference_values = numpy.array( + [ + # fmt: off + 0.09992, 0.09992, 0.08598, 0.10866, 0.10866, + 0.08704, 0.10761, 0.10761, 0.11872, 0.11872, + 0.10678, 0.10678, 0.10158, 0.09112, 0.10342, + 0.10342, 0.10319, 0.10608, 0.11005 + # fmt: on + ] + ) + assert reference_values == approx(test_vals, abs=0.01) + + # No-exact sem + setup_free_energy_calculation.exact_sem_each_ti_fraction = False + setup_free_energy_calculation.compute_free_energy(seed=random_seed) + results = setup_free_energy_calculation.results + + # Test `ti-nocor` free energies and uncertainties + test_vals = [ + results["attach"][method]["fe"].magnitude, + results["attach"][method]["sem"].magnitude, + results["pull"][method]["fe"].magnitude, + results["pull"][method]["sem"].magnitude, + ] + reference_values = numpy.array([13.34, 0.098, -1.71, 0.56]) + assert reference_values == approx(test_vals, abs=0.01) + + # ROI only runs during TI. + + # Test attach `ti-nocor` largest_neighbor values + test_vals = results["attach"][method]["largest_neighbor"].magnitude + reference_values = numpy.array([0.0138, 0.0362, 0.055, 0.0657, 0.0657]) + assert reference_values == approx(test_vals, abs=0.01) + + # Test pull `ti-nocor` largest_neighbor values + test_vals = results["pull"][method]["largest_neighbor"].magnitude + + reference_values = numpy.array( + [ + # fmt: off + 0.09992, 0.09992, 0.08598, 0.10867, 0.10866, + 0.08704, 0.10761, 0.10761, 0.11872, 0.11872, + 0.10678, 0.10678, 0.10158, 0.09112, 0.10342, + 0.10342, 0.10319, 0.10608, 0.11006 + # fmt: on + ] + ) + assert reference_values == approx(test_vals, abs=0.01) + + +def test_reference_state_work(clean_files, setup_free_energy_calculation): + results = setup_free_energy_calculation.results + assert numpy.isclose(-4.34372, results["ref_state_work"].magnitude) + + +def test_save_and_loading(clean_files, setup_free_energy_calculation): + # Test save and load results -- JSON + results = deepcopy(setup_free_energy_calculation.results) + setup_free_energy_calculation.save_results("tmp/results.json", overwrite=True) + assert is_file_and_not_empty("tmp/results.json") + setup_free_energy_calculation.results = {} + setup_free_energy_calculation.load_results("tmp/results.json") + assert len(setup_free_energy_calculation.results) == len(results) + + # Test save and load simulation data -- JSON + setup_free_energy_calculation.save_simulation_data_to_json( + "tmp/simulation.json", overwrite=True + ) + assert is_file_and_not_empty("tmp/simulation.json") + setup_free_energy_calculation.changing_restraints = None + setup_free_energy_calculation.orders = None + setup_free_energy_calculation.simulation_data = None + setup_free_energy_calculation.load_simulation_data_from_json("tmp/simulation.json") + assert setup_free_energy_calculation.changing_restraints is not None + assert setup_free_energy_calculation.orders is not None + assert setup_free_energy_calculation.simulation_data is not None def test_temperature(clean_files): - input_pdb = pmd.load_file( + input_pdb = parmed.load_file( os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-dum.pdb") ) @@ -320,71 +474,69 @@ def test_bootstrap(): """Test the utility modules in `analysis`""" # Test regression statistics - x = np.linspace(0, 10, 11) - y = np.linspace(0, 10, 11) + x = numpy.linspace(0, 10, 11) + y = numpy.linspace(0, 10, 11) stats = analysis.summarize_statistics(x, y) - assert all(stats == np.array([1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0])) + assert all(stats == numpy.array([1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0])) # Test Regression bootstrap - x_sem = np.ones_like(x) - y_sem = np.ones_like(y) + x_sem = numpy.ones_like(x) + y_sem = numpy.ones_like(y) - np.random.seed(0) + numpy.random.seed(0) results = analysis.regression_bootstrap(x, x_sem, y, y_sem, cycles=1) + # fmt: off compare = { - "slope": 0.7703030693742714, - "intercept": 0.7777706430286471, - "R": 0.9184956779750857, - "R**2": 0.8436343104589124, - "RMSE": 1.637675285465762, - "MSE": -0.40918918963339473, - "MUE": 1.2814477376354574, - "Tau": 0.8545454545454545, + "slope": 0.77030, + "intercept": 0.77778, + "R": 0.91850, + "R**2": 0.84363, + "RMSE": 1.63768, + "MSE": -0.40919, + "MUE": 1.28145, + "Tau": 0.85454, } + # fmt: on for stat in results["mean"]: assert pytest.approx(results["mean"][stat], abs=1e-3) == compare[stat] # Test dG Bootstrap - np.random.seed(0) + numpy.random.seed(0) results = analysis.dG_bootstrap(-12, 2, -12, 2, cycles=1, with_uncertainty=True) - assert pytest.approx(results["mean"], abs=1e-3) == -11.205587976469898 + assert pytest.approx(results["mean"], abs=1e-3) == -11.20559 assert pytest.approx(results["sem"], abs=1e-3) == 0.0 - assert pytest.approx(results["ci"][0], abs=1e-3) == -11.20558798 - assert pytest.approx(results["ci"][1], abs=1e-3) == -11.20558798 + assert pytest.approx(results["ci"][0], abs=1e-3) == -11.20559 + assert pytest.approx(results["ci"][1], abs=1e-3) == -11.20559 results = analysis.dG_bootstrap(-12, 2, -12, 2, cycles=1, with_uncertainty=False) - assert pytest.approx(results["mean"], abs=1e-3) == -0.4106792724182964 + assert pytest.approx(results["mean"], abs=1e-3) == -0.41068 assert pytest.approx(results["sem"], abs=1e-3) == 0.0 - assert pytest.approx(results["ci"][0], abs=1e-3) == -0.41067927 - assert pytest.approx(results["ci"][1], abs=1e-3) == -0.41067927 + assert pytest.approx(results["ci"][0], abs=1e-3) == -0.41068 + assert pytest.approx(results["ci"][1], abs=1e-3) == -0.41068 # Test dH Bootstrap - np.random.seed(0) + numpy.random.seed(0) + # fmt: off results = analysis.dH_bootstrap( - -15, - 2, - -15, - 2, - -10, - 2, - -10, - 2, + -15, 2, -15, 2, -10, 2, -10, 2, cycles=1, with_uncertainty=True, ) - assert pytest.approx(results["mean"], abs=1e-3) == -11.50986102184404 + # fmt: on + assert pytest.approx(results["mean"], abs=1e-3) == -11.50986 assert pytest.approx(results["sem"], abs=1e-3) == 0.0 - assert pytest.approx(results["ci"][0], abs=1e-3) == -11.5098610 - assert pytest.approx(results["ci"][1], abs=1e-3) == -11.5098610 + assert pytest.approx(results["ci"][0], abs=1e-3) == -11.50986 + assert pytest.approx(results["ci"][1], abs=1e-3) == -11.50986 def test_utils(): """Test the utility modules in `analysis`""" assert utils.get_factors(10) == [1, 2, 5, 10] assert utils.get_nearest_max(100) == 90 - np.random.seed(0) - results = utils.get_block_sem(np.random.normal(10.0, 2.0, 100)) - assert pytest.approx(results, abs=1e-3) == 0.3916368835724714 + + numpy.random.seed(0) + results = utils.get_block_sem(numpy.random.normal(10.0, 2.0, 100)) + assert pytest.approx(results, abs=1e-3) == 0.39164 assert utils.get_subsampled_indices(10, 2.0) == [0, 2, 4, 6, 8] diff --git a/paprika/tests/test_evaluator.py b/paprika/tests/test_evaluator.py index 11fd4fa2..d4cadc36 100644 --- a/paprika/tests/test_evaluator.py +++ b/paprika/tests/test_evaluator.py @@ -5,15 +5,18 @@ import os import shutil -import numpy as np -import parmed as pmd +import numpy +import parmed import pytest -import pytraj as pt +import pytraj +import yaml from openff.units import unit as openff_unit from paprika.evaluator import Analyze, Setup from paprika.evaluator.amber import generate_gaff from paprika.restraints import DAT_restraint +from paprika.taproom.taproom import read_yaml_schema +from paprika.taproom.utils import convert_string_to_quantity, de_alias logger = logging.getLogger(__name__) @@ -29,6 +32,300 @@ def clean_files(directory=os.path.join(os.path.dirname(__file__), "tmp")): shutil.rmtree(directory) +@pytest.fixture() +def complex_file(): + complex_pdb = os.path.join(os.path.dirname(__file__), "../data/cb6-but/vac.pdb") + + butane_molecule = [] + structure = parmed.load_file(complex_pdb, structure=True) + for atom in structure.topology.atoms(): + if atom.residue.name == "BUT": + butane_molecule.append(atom.index) + + G1 = ":BUT@C" + G2 = ":BUT@C3" + + host_guest = Setup.prepare_complex_structure( + complex_pdb, + butane_molecule, + f"{G1} {G2}", + 24.0, + 0, + 46, + ) + + Setup.add_dummy_atoms_to_structure( + host_guest, + [ + numpy.array([0, 0, 0]), + numpy.array([0, 0, -3.0]), + numpy.array([0, 2.2, -5.2]), + ], + numpy.zeros(3), + ) + + return host_guest + + +@pytest.fixture(scope="module") +def yaml_restraint_schema(): + yaml_file = """name: bam +structure: bam.mol2 +complex: a-bam.pdb +net_charge: +1e +aliases: + - D1: :DM1 + - D2: :DM2 + - D3: :DM3 + - G1: :BAM@C4 + - G2: :BAM@N1 +restraints: + guest: + - restraint: + atoms: D1 G1 + attach: + # During the 'attach' phase, the `force_constant` argument is the + # final force constant. + force_constant: 5.0 * kilocalorie / mole / angstrom**2 + target: 6.0 * angstrom + pull: + # During the 'pull' phase, the `target` argument is the final value of + # the restraint. + force_constant: 5.0 * kilocalorie / mole / angstrom**2 + target: 24.0 * angstrom + - restraint: + atoms: D2 D1 G1 + attach: + force_constant: 100.0 * kilocalorie / mole / radians**2 + target: 180.0 * degrees + pull: + force_constant: 100.0 * kilocalorie / mole / radians**2 + target: 180.0 * degrees + - restraint: + atoms: D1 G1 G2 + attach: + force_constant: 100.0 * kilocalorie / mole / radians**2 + target: 180.0 * degrees + pull: + force_constant: 100.0 * kilocalorie / mole / radians**2 + target: 180.0 * degrees + + wall_restraints: + - restraint: + atoms: ":1@O2 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 9.3 * angstrom + - restraint: + atoms: ":2@O2 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 9.3 * angstrom + - restraint: + atoms: ":3@O2 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 9.3 * angstrom + - restraint: + atoms: ":4@O2 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 9.3 * angstrom + - restraint: + atoms: ":5@O2 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 9.3 * angstrom + - restraint: + atoms: ":6@O2 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 9.3 * angstrom + - restraint: + atoms: ":1@O6 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 11.3 * angstrom + - restraint: + atoms: ":2@O6 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 11.3 * angstrom + - restraint: + atoms: ":3@O6 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 11.3 * angstrom + - restraint: + atoms: ":4@O6 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 11.3 * angstrom + - restraint: + atoms: ":5@O6 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 11.3 * angstrom + - restraint: + atoms: ":6@O6 G1" + force_constant: 50.0 * kilocalorie / mole / angstrom**2 + target: 11.3 * angstrom + +symmetry_correction: + restraints: + - restraint: + atoms: D2 G1 G2 + force_constant: 200.0 * kilocalorie / mole / radian**2 + target: 91 * degrees + # Do not attempt to automatically correct for the symmetry restraint by adding -RT \ln (microstates). + # Instead, we will apply the symmetry restraint, which locks in a particular binding orientation, and then + # perform separate calculations. + microstates: 1""" + return yaml_file + + +@pytest.fixture(scope="module") +def restraints_schema(): + schema = { + "static": [ + { + "atoms": ":DM1 :CB6@O", + "force_constant": 5.0 + * openff_unit.kcal + / openff_unit.mol + / openff_unit.angstrom**2, + } + ], + "conformational": [ + { + "atoms": ":CB6@O :CB6@O2 :CB6@O4 :CB6@O6", + "force_constant": 6.0 + * openff_unit.kcal + / openff_unit.mol + / openff_unit.radians**2, + "target": 104.3 * openff_unit.degrees, + } + ], + "symmetry": [ + { + "atoms": ":DM2 :BUT@C :BUT@C3", + "force_constant": 50.0 + * openff_unit.kcal + / openff_unit.mol + / openff_unit.radians**2, + "target": 11.0 * openff_unit.degrees, + } + ], + "wall": [ + { + "atoms": ":CB6@O :BUT@C", + "force_constant": 50.0 + * openff_unit.kcal + / openff_unit.mol + / openff_unit.angstrom**2, + "target": 11.0 * openff_unit.angstrom, + } + ], + "guest": [ + { + "atoms": ":DM1 :BUT@C", + "attach": { + "force_constant": 5.0 + * openff_unit.kcal + / openff_unit.mol + / openff_unit.angstrom**2, + "target": 6.0 * openff_unit.angstrom, + }, + "pull": { + "force_constant": 5.0 + * openff_unit.kcal + / openff_unit.mol + / openff_unit.angstrom**2, + "target": 24.0 * openff_unit.angstrom, + }, + } + ], + } + return schema + + +def test_taproom_yaml(clean_files, yaml_restraint_schema): + temporary_directory = os.path.join(os.path.dirname(__file__), "tmp") + + k_dist = 5.0 * openff_unit.kcal / openff_unit.mole / openff_unit.angstrom**2 + k_wall = 50.0 * openff_unit.kcal / openff_unit.mole / openff_unit.angstrom**2 + k_angle = 100.0 * openff_unit.kcal / openff_unit.mole / openff_unit.radian**2 + r_initial = 6.0 * openff_unit.angstrom + r_final = 24.0 * openff_unit.angstrom + angle = 180.0 * openff_unit.degrees + + # Write yaml string to file + with open(f"{temporary_directory}/guest.yaml", "w") as f: + f.write(yaml_restraint_schema) + + # Check de_alias + with open(f"{temporary_directory}/guest.yaml", "r") as f: + yaml_data = yaml.safe_load(f) + + if "aliases" in yaml_data.keys(): + yaml_data = de_alias(yaml_data) + + distance_restraint = yaml_data["restraints"]["guest"][0]["restraint"] + theta_restraint = yaml_data["restraints"]["guest"][1]["restraint"] + beta_restraint = yaml_data["restraints"]["guest"][2]["restraint"] + symmetry_restraint = yaml_data["symmetry_correction"]["restraints"][0]["restraint"] + + assert distance_restraint["atoms"] == ":DM1 :BAM@C4" + assert theta_restraint["atoms"] == ":DM2 :DM1 :BAM@C4" + assert beta_restraint["atoms"] == ":DM1 :BAM@C4 :BAM@N1" + assert symmetry_restraint["atoms"] == ":DM2 :BAM@C4 :BAM@N1" + + # Check convert string to quantity + convert_string_to_quantity(yaml_data) + + assert distance_restraint["attach"]["force_constant"] == k_dist + assert theta_restraint["attach"]["force_constant"] == k_angle + assert beta_restraint["attach"]["force_constant"] == k_angle + + assert distance_restraint["pull"]["force_constant"] == k_dist + assert theta_restraint["pull"]["force_constant"] == k_angle + assert beta_restraint["pull"]["force_constant"] == k_angle + + assert distance_restraint["attach"]["target"] == r_initial + assert distance_restraint["pull"]["target"] == r_final + + assert theta_restraint["attach"]["target"] == angle + assert theta_restraint["pull"]["target"] == angle + assert beta_restraint["attach"]["target"] == angle + assert beta_restraint["pull"]["target"] == angle + + for i, restraint in enumerate(yaml_data["restraints"]["wall_restraints"]): + assert restraint["restraint"]["force_constant"] == k_wall + if i < 6: + assert restraint["restraint"]["target"] == 9.3 * openff_unit.angstrom + else: + assert restraint["restraint"]["target"] == 11.3 * openff_unit.angstrom + + # Check full conversion + guest_spec = read_yaml_schema(f"{temporary_directory}/guest.yaml") + + distance_restraint = guest_spec["restraints"]["guest"][0]["restraint"] + theta_restraint = guest_spec["restraints"]["guest"][1]["restraint"] + beta_restraint = guest_spec["restraints"]["guest"][2]["restraint"] + + assert distance_restraint["attach"]["force_constant"] == k_dist + assert theta_restraint["attach"]["force_constant"] == k_angle + assert beta_restraint["attach"]["force_constant"] == k_angle + + assert distance_restraint["pull"]["force_constant"] == k_dist + assert theta_restraint["pull"]["force_constant"] == k_angle + assert beta_restraint["pull"]["force_constant"] == k_angle + + assert distance_restraint["attach"]["target"] == r_initial + assert distance_restraint["pull"]["target"] == r_final + + assert theta_restraint["attach"]["target"] == angle + assert theta_restraint["pull"]["target"] == angle + assert beta_restraint["attach"]["target"] == angle + assert beta_restraint["pull"]["target"] == angle + + for i, restraint in enumerate(guest_spec["restraints"]["wall_restraints"]): + assert restraint["restraint"]["force_constant"] == k_wall + if i < 6: + assert restraint["restraint"]["target"] == 9.3 * openff_unit.angstrom + else: + assert restraint["restraint"]["target"] == 11.3 * openff_unit.angstrom + + def test_evaluator_setup_structure(clean_files): temporary_directory = os.path.join(os.path.dirname(__file__), "tmp") @@ -38,7 +335,7 @@ def test_evaluator_setup_structure(clean_files): # Test prepare_complex_structure host_guest_pdb = os.path.join(os.path.dirname(__file__), "../data/cb6-but/vac.pdb") guest_atom_indices = [] - structure = pmd.load_file(host_guest_pdb, structure=True) + structure = parmed.load_file(host_guest_pdb, structure=True) for atom in structure.topology.atoms(): if atom.residue.name == "BUT": guest_atom_indices.append(atom.index) @@ -67,12 +364,14 @@ def test_evaluator_setup_structure(clean_files): cG1 = host_guest_structure_final[G1].coordinates[0] cG2 = host_guest_structure_final[G2].coordinates[0] vec = cG2 - cG1 - axis = np.array([0, 0, 1]) - theta = np.arccos(np.dot(vec, axis) / (np.linalg.norm(vec) * np.linalg.norm(axis))) + axis = numpy.array([0, 0, 1]) + theta = numpy.arccos( + numpy.dot(vec, axis) / (numpy.linalg.norm(vec) * numpy.linalg.norm(axis)) + ) assert theta == 0.0 # Test prepare_host_structure - structure = pmd.load_file(host_guest_pdb, structure=True) + structure = parmed.load_file(host_guest_pdb, structure=True) structure[":CB6"].save(os.path.join(temporary_directory, "cb6.pdb")) host_pdb = os.path.join(temporary_directory, "cb6.pdb") host_atom_indices = [] @@ -84,17 +383,28 @@ def test_evaluator_setup_structure(clean_files): host_pdb, host_atom_indices, ) - center_of_mass = pmd.geometry.center_of_mass( - host_structure.coordinates, masses=np.ones(len(host_structure.coordinates)) + center_of_mass = parmed.geometry.center_of_mass( + host_structure.coordinates, masses=numpy.ones(len(host_structure.coordinates)) + ) + assert pytest.approx(center_of_mass[0], abs=1e-3) == 0.0 + assert pytest.approx(center_of_mass[1], abs=1e-3) == 0.0 + assert pytest.approx(center_of_mass[2], abs=1e-3) == 0.0 + + host_structure = Setup.prepare_host_structure( + host_pdb, + host_atom_indices=None, + ) + center_of_mass = parmed.geometry.center_of_mass( + host_structure.coordinates, masses=numpy.ones(len(host_structure.coordinates)) ) assert pytest.approx(center_of_mass[0], abs=1e-3) == 0.0 assert pytest.approx(center_of_mass[1], abs=1e-3) == 0.0 assert pytest.approx(center_of_mass[2], abs=1e-3) == 0.0 - inertia_tensor = np.dot( + inertia_tensor = numpy.dot( host_structure.coordinates.transpose(), host_structure.coordinates ) - eig_val, eig_vec = np.linalg.eig(inertia_tensor) + eig_val, eig_vec = numpy.linalg.eig(inertia_tensor) assert pytest.approx(eig_vec[0, -1], abs=1e-3) == 0.0 assert pytest.approx(eig_vec[1, -1], abs=1e-3) == 0.0 assert pytest.approx(eig_vec[2, -1], abs=1e-3) == 1.0 @@ -103,11 +413,11 @@ def test_evaluator_setup_structure(clean_files): Setup.add_dummy_atoms_to_structure( host_structure, [ - np.array([0, 0, 0]), - np.array([0, 0, -3.0]), - np.array([0, 2.2, -5.2]), + numpy.array([0, 0, 0]), + numpy.array([0, 0, -3.0]), + numpy.array([0, 2.2, -5.2]), ], - np.zeros(3), + numpy.zeros(3), ) dummy_atoms = [] for atom in host_structure.topology.atoms(): @@ -121,9 +431,135 @@ def test_evaluator_setup_structure(clean_files): assert pytest.approx(host_structure[":DM3"].coordinates[0][2], abs=1e-3) == -5.2 +def test_evaluator_setup_restraints(clean_files, complex_file, restraints_schema): + complex_file_path = os.path.join(os.path.dirname(__file__), "complex.pdb") + complex_file.save(complex_file_path, overwrite=True) + + attach_lambdas = list(numpy.linspace(0, 1, 15)) + n_attach = 15 + n_pull = 46 + n_release = 15 + + static_restraints = Setup.build_static_restraints( + complex_file_path, + n_attach, + n_pull, + n_release, + restraints_schema["static"], + use_amber_indices=False, + ) + assert static_restraints[0].mask1 == ":DM1" + assert static_restraints[0].mask2 == ":CB6@O" + for k in static_restraints[0].phase["attach"]["force_constants"]: + assert k == restraints_schema["static"][0]["force_constant"] + + conformational_restraint = Setup.build_conformational_restraints( + complex_file_path, + attach_lambdas, + None, + None, + restraints_schema["conformational"], + use_amber_indices=False, + ) + assert conformational_restraint[0].mask1 == ":CB6@O" + assert conformational_restraint[0].mask2 == ":CB6@O2" + assert conformational_restraint[0].mask3 == ":CB6@O4" + assert conformational_restraint[0].mask4 == ":CB6@O6" + for i, (k, target) in enumerate( + zip( + conformational_restraint[0].phase["attach"]["force_constants"], + conformational_restraint[0].phase["attach"]["targets"], + ) + ): + assert ( + k + == attach_lambdas[i] + * restraints_schema["conformational"][0]["force_constant"] + ) + assert target == restraints_schema["conformational"][0]["target"] + + symmetry_restraints = Setup.build_symmetry_restraints( + complex_file_path, + n_attach, + restraints_schema["symmetry"], + use_amber_indices=False, + ) + assert symmetry_restraints[0].mask1 == ":DM2" + assert symmetry_restraints[0].mask2 == ":BUT@C" + assert symmetry_restraints[0].mask3 == ":BUT@C3" + assert symmetry_restraints[0].custom_restraint_values["r1"] is None + assert symmetry_restraints[0].custom_restraint_values["r2"] is None + assert ( + symmetry_restraints[0].custom_restraint_values["r3"] + == 0.0 * openff_unit.degrees + ) + assert ( + symmetry_restraints[0].custom_restraint_values["r4"] + == 0.0 * openff_unit.degrees + ) + assert ( + symmetry_restraints[0].custom_restraint_values["rk2"] + == restraints_schema["symmetry"][0]["force_constant"] + ) + assert ( + symmetry_restraints[0].custom_restraint_values["rk3"] + == restraints_schema["symmetry"][0]["force_constant"] + ) + + wall_restraints = Setup.build_wall_restraints( + complex_file_path, + n_attach, + restraints_schema["wall"], + use_amber_indices=False, + ) + assert wall_restraints[0].mask1 == ":CB6@O" + assert wall_restraints[0].mask2 == ":BUT@C" + assert ( + wall_restraints[0].custom_restraint_values["r1"] == 0.0 * openff_unit.angstrom + ) + assert ( + wall_restraints[0].custom_restraint_values["r2"] == 0.0 * openff_unit.angstrom + ) + assert wall_restraints[0].custom_restraint_values["r3"] is None + assert wall_restraints[0].custom_restraint_values["r4"] is None + assert ( + wall_restraints[0].custom_restraint_values["rk2"] + == restraints_schema["wall"][0]["force_constant"] + ) + assert ( + wall_restraints[0].custom_restraint_values["rk3"] + == restraints_schema["wall"][0]["force_constant"] + ) + + guest_restraints = Setup.build_guest_restraints( + complex_file_path, + attach_lambdas, + None, + restraints_schema["guest"], + use_amber_indices=False, + ) + assert guest_restraints[0].mask1 == ":DM1" + assert guest_restraints[0].mask2 == ":BUT@C" + for i, (k, target) in enumerate( + zip( + guest_restraints[0].phase["attach"]["force_constants"], + guest_restraints[0].phase["attach"]["targets"], + ) + ): + assert ( + k + == attach_lambdas[i] + * restraints_schema["guest"][0]["attach"]["force_constant"] + ) + assert target == restraints_schema["guest"][0]["attach"]["target"] + + def test_evaluator_analyze(clean_files): + # ---------------------------------------------------------------- # + # Configure system and restraints + # ---------------------------------------------------------------- # input_pdb = os.path.join(os.path.dirname(__file__), "../data/cb6-but/vac.pdb") - structure = pmd.load_file(input_pdb, structure=True) + structure = parmed.load_file(input_pdb, structure=True) guest_atom_indices = [] for atom in structure.topology.atoms(): @@ -141,13 +577,20 @@ def test_evaluator_analyze(clean_files): Setup.add_dummy_atoms_to_structure( host_guest_structure, [ - np.array([0, 0, 0]), - np.array([0, 0, -3.0]), - np.array([0, 2.2, -5.2]), + numpy.array([0, 0, 0]), + numpy.array([0, 0, -3.0]), + numpy.array([0, 2.2, -5.2]), ], - np.zeros(3), + numpy.zeros(3), ) + angle = 180 * openff_unit.degrees + k_angle = 100 * openff_unit.kcal / openff_unit.mole / openff_unit.radians**2 + r_initial = 6.0 * openff_unit.angstrom + r_final = 24.0 * openff_unit.angstrom + k_r = 5.0 * openff_unit.kcal / openff_unit.mole / openff_unit.angstrom**2 + attach_fractions = [0.00, 0.04, 0.181, 0.496, 1.000] + # Distance restraint rest1 = DAT_restraint() rest1.continuous_apr = True @@ -155,12 +598,12 @@ def test_evaluator_analyze(clean_files): rest1.topology = host_guest_structure rest1.mask1 = ":DM1" rest1.mask2 = ":BUT@C" - rest1.attach["target"] = 6.0 - rest1.attach["fraction_list"] = [0.00, 0.04, 0.181, 0.496, 1.000] - rest1.attach["fc_final"] = 5.0 + rest1.attach["target"] = r_initial + rest1.attach["fraction_list"] = attach_fractions + rest1.attach["fc_final"] = k_r rest1.pull["fc"] = rest1.attach["fc_final"] rest1.pull["target_initial"] = rest1.attach["target"] - rest1.pull["target_final"] = 24.0 + rest1.pull["target_final"] = r_final rest1.pull["num_windows"] = 19 rest1.initialize() @@ -172,9 +615,9 @@ def test_evaluator_analyze(clean_files): rest2.mask1 = ":DM2" rest2.mask2 = ":DM1" rest2.mask3 = ":BUT@C" - rest2.attach["target"] = 180.0 - rest2.attach["fraction_list"] = [0.00, 0.04, 0.181, 0.496, 1.000] - rest2.attach["fc_final"] = 100.0 + rest2.attach["target"] = angle + rest2.attach["fraction_list"] = attach_fractions + rest2.attach["fc_final"] = k_angle rest2.pull["fc"] = rest2.attach["fc_final"] rest2.pull["target_initial"] = rest2.attach["target"] rest2.pull["target_final"] = rest2.attach["target"] @@ -189,30 +632,94 @@ def test_evaluator_analyze(clean_files): rest3.mask1 = ":DM1" rest3.mask2 = ":BUT@C" rest3.mask3 = ":BUT@C3" - rest3.attach["target"] = 180.0 - rest3.attach["fraction_list"] = [0.00, 0.04, 0.181, 0.496, 1.000] - rest3.attach["fc_final"] = 100.0 + rest3.attach["target"] = angle + rest3.attach["fraction_list"] = attach_fractions + rest3.attach["fc_final"] = k_angle rest3.pull["fc"] = rest2.attach["fc_final"] rest3.pull["target_initial"] = rest2.attach["target"] rest3.pull["target_final"] = rest2.attach["target"] rest3.pull["num_windows"] = 19 rest3.initialize() - temperature = 298.15 + # Angle 2 + rest4 = DAT_restraint() + rest4.continuous_apr = True + rest4.amber_index = True + rest4.topology = input_pdb + rest4.mask1 = ":DM1" + rest4.mask2 = ":BUT@C" + rest4.mask3 = ":BUT@C2" + rest4.mask4 = ":BUT@C3" + rest4.attach["target"] = angle + rest4.attach["fraction_list"] = attach_fractions + rest4.attach["fc_final"] = k_angle + rest4.pull["fc"] = rest2.attach["fc_final"] + rest4.pull["target_initial"] = rest2.attach["target"] + rest4.pull["target_final"] = rest2.attach["target"] + rest4.pull["num_windows"] = 19 + rest4.initialize() + + temperature = 298.15 * openff_unit.kelvin guest_restraints = [rest1, rest2, rest3] - ref_state_work = Analyze.compute_ref_state_work(temperature, guest_restraints) + + # ---------------------------------------------------------------- # + # Test reference state work + # ---------------------------------------------------------------- # + ref_state_work = Analyze.compute_ref_state_work(temperature, [rest1, rest2, rest3]) assert ( pytest.approx( ref_state_work.to(openff_unit.kcal / openff_unit.mole).magnitude, abs=1e-3 ) == -7.14151 ) + ref_state_work = Analyze.compute_ref_state_work( + temperature, [rest1, rest2, rest3, rest4] + ) + assert ( + pytest.approx( + ref_state_work.to(openff_unit.kcal / openff_unit.mole).magnitude, abs=1e-3 + ) + == -9.41062 + ) + # ---------------------------------------------------------------- # + # Test symmetry correction + # ---------------------------------------------------------------- # fe_sym = Analyze.symmetry_correction(n_microstates=1, temperature=298.15) assert fe_sym == 0.0 fe_sym = Analyze.symmetry_correction(n_microstates=2, temperature=298.15) assert pytest.approx(fe_sym, abs=1e-3) == -0.410679 + # ---------------------------------------------------------------- # + # Test APR phase + # ---------------------------------------------------------------- # + guest_restraints[0].mask1 = ":CB6@O" + guest_restraints[1].mask1 = ":CB6@O2" + guest_restraints[1].mask2 = ":CB6@O" + guest_restraints[2].mask1 = ":CB6@O" + for restraint in guest_restraints: + restraint.initialize() + + apr_data_path = os.path.join(os.path.dirname(__file__), "../data/cb6-but-apr/") + tmp_path = os.path.join(os.path.dirname(__file__), "tmp") + window_list = ["a000", "a001", "a002", "a003"] + [f"p{i:03}" for i in range(19)] + for window in window_list: + shutil.copytree(f"{apr_data_path}/{window}", f"{tmp_path}/{window}") + shutil.copy(f"{apr_data_path}/vac.pdb", f"{tmp_path}/{window}/vac.pdb") + + results = Analyze.compute_phase_free_energy( + phase="attach", + restraints=guest_restraints, + windows_directory=tmp_path, + topology_name="vac.pdb", + trajectory_mask="*.nc", + analysis_method="ti-block", + ) + + # loose comparison due to short trajectory + assert pytest.approx(results["attach"]["ti-block"]["fe"].magnitude, abs=2) == 946 + assert pytest.approx(results["attach"]["ti-block"]["sem"].magnitude, abs=2) == 10 + def test_evaluator_gaff(clean_files): temporary_directory = os.path.join(os.path.dirname(__file__), "tmp") @@ -229,12 +736,15 @@ def test_evaluator_gaff(clean_files): directory_path=temporary_directory, ) - structure = pt.iterload( + structure = pytraj.iterload( os.path.join(temporary_directory, f"but.{gaff_version}.mol2") ) # fmt: off - butane_atom_type = ["c3", "hc", "hc", "hc", "c3", "hc", "c3", "hc", "hc", "hc", "c3", "hc", "hc", "hc"] + butane_atom_type = [ + "c3", "hc", "hc", "hc", "c3", "hc", "c3", + "hc", "hc", "hc", "c3", "hc", "hc", "hc", + ] # fmt: on residue_names = [] diff --git a/paprika/tests/test_restraints.py b/paprika/tests/test_restraints.py index 587a180d..cda75d12 100644 --- a/paprika/tests/test_restraints.py +++ b/paprika/tests/test_restraints.py @@ -1267,6 +1267,78 @@ def test_restraints_output_modules(clean_files): r.initialize() guest_restraints.append(r) + # Guest - Dihedral + r = DAT_restraint() + r.amber_index = True + r.continuous_apr = True + r.auto_apr = True + r.topology = os.path.join( + os.path.dirname(__file__), "../data/cb6-but/cb6-but-dum.pdb" + ) + r.mask1 = ":CB6@C13" + r.mask2 = ":CB6@C3" + r.mask3 = ":CB6@C" + r.mask4 = ":BUT@C1" + r.attach["target"] = 180.0 + r.attach["num_windows"] = 15 + r.attach["fc_initial"] = 0.0 + r.attach["fc_final"] = 100.0 + r.pull["fc"] = r.attach["fc_final"] + r.pull["num_windows"] = 46 + r.pull["target_initial"] = r.attach["target"] + r.pull["target_final"] = 180.0 + r.release["target"] = r.pull["target_final"] + r.release["num_windows"] = r.attach["num_windows"] + r.release["fc_initial"] = r.attach["fc_initial"] + r.release["fc_final"] = r.attach["fc_final"] + r.initialize() + guest_restraints.append(r) + + # Guest - Group + r = DAT_restraint() + r.amber_index = True + r.continuous_apr = True + r.auto_apr = True + r.topology = os.path.join( + os.path.dirname(__file__), "../data/cb6-but/cb6-but-dum.pdb" + ) + r.mask1 = ":CB6@C13,C12" + r.mask2 = ":CB6@C3,C2" + r.mask3 = ":CB6@C,C1" + r.mask4 = ":BUT@C1,C2" + r.attach["target"] = 180.0 + r.attach["num_windows"] = 15 + r.attach["fc_initial"] = 0.0 + r.attach["fc_final"] = 100.0 + r.pull["fc"] = r.attach["fc_final"] + r.pull["num_windows"] = 46 + r.pull["target_initial"] = r.attach["target"] + r.pull["target_final"] = 180.0 + r.release["target"] = r.pull["target_final"] + r.release["num_windows"] = r.attach["num_windows"] + r.release["fc_initial"] = r.attach["fc_initial"] + r.release["fc_final"] = r.attach["fc_final"] + r.initialize() + guest_restraints.append(r) + + # Guest - None + r_empty = DAT_restraint() + r_empty.amber_index = True + r_empty.continuous_apr = True + r_empty.auto_apr = True + r_empty.topology = os.path.join( + os.path.dirname(__file__), "../data/cb6-but/cb6-but-dum.pdb" + ) + r_empty.mask1 = ":CB6@C13,C12" + r_empty.mask2 = ":CB6@C3,C2" + r_empty.mask3 = ":CB6@C,C1" + r_empty.mask4 = ":BUT@C1,C2" + r_empty.attach["target"] = 180.0 + r_empty.attach["num_windows"] = 15 + r_empty.attach["fc_initial"] = 0.0 + r_empty.attach["fc_final"] = 100.0 + r_empty.initialize() + # Create OpenMM System system = structure.createSystem( nonbondedMethod=app.NoCutoff, @@ -1296,7 +1368,9 @@ def test_restraints_output_modules(clean_files): or isinstance(force, openmm.CustomTorsionForce) ] + # ---------------------------------------------------------------- # # Test dummy atom positional restraint + # ---------------------------------------------------------------- # for i, force in enumerate(positional_restraints): particle, parameters = force.getParticleParameters(0) assert pytest.approx(parameters[0], abs=1e-3) == k_pos.value_in_unit( @@ -1306,7 +1380,9 @@ def test_restraints_output_modules(clean_files): assert pytest.approx(parameters[2], abs=1e-3) == dummy_atoms[i]["y"] / 10 assert pytest.approx(parameters[3], abs=1e-3) == dummy_atoms[i]["z"] / 10 + # ---------------------------------------------------------------- # # Test Amber NMR-style restraints + # ---------------------------------------------------------------- # atom1, atom2, parameters = DAT_restraint_list[0].getBondParameters(0) assert pytest.approx(parameters[0]) == 2092.0 assert pytest.approx(parameters[1]) == 0.6 @@ -1319,7 +1395,9 @@ def test_restraints_output_modules(clean_files): assert pytest.approx(parameters[0]) == 418.4 assert pytest.approx(parameters[1]) == numpy.pi + # ---------------------------------------------------------------- # # Test Amber restraints + # ---------------------------------------------------------------- # r_string = amber_restraint_line(guest_restraints[0], window) assert r_string.split()[2].split(",")[0] == "123" assert r_string.split()[2].split(",")[1] == "119" @@ -1341,7 +1419,38 @@ def test_restraints_output_modules(clean_files): assert float(theta_string.split()[12].split(",")[0]) == 100.0 assert float(theta_string.split()[14].split(",")[0]) == 100.0 + dihedral_string = amber_restraint_line(guest_restraints[3], window) + assert dihedral_string.split()[2].split(",")[0] == "38" + assert dihedral_string.split()[2].split(",")[1] == "10" + assert dihedral_string.split()[2].split(",")[2] == "1" + assert dihedral_string.split()[2].split(",")[3] == "113" + + group_string = amber_restraint_line(guest_restraints[4], window) + assert group_string.split()[2].split(",")[0] == "-1" + assert group_string.split()[2].split(",")[1] == "-1" + assert group_string.split()[2].split(",")[2] == "-1" + assert group_string.split()[2].split(",")[3] == "-1" + assert float(group_string.split()[4].split(",")[0]) == 0.0 + assert float(group_string.split()[6].split(",")[0]) == 180.0 + assert float(group_string.split()[8].split(",")[0]) == 180.0 + assert float(group_string.split()[10].split(",")[0]) == 360.0 + assert float(group_string.split()[12].split(",")[0]) == 100.0 + assert float(group_string.split()[14].split(",")[0]) == 100.0 + assert group_string.split()[15] == "igr1=" + assert group_string.split()[16] == "37,38," + assert group_string.split()[17] == "igr2=" + assert group_string.split()[18] == "9,10," + assert group_string.split()[19] == "igr3=" + assert group_string.split()[20] == "1,2," + assert group_string.split()[21] == "igr4=" + assert group_string.split()[22] == "113,115," + + empty_string = amber_restraint_line(r_empty, window) + assert empty_string == "" + + # ---------------------------------------------------------------- # # Test Plumed output + # ---------------------------------------------------------------- # plumed = Plumed() plumed.path = "tmp" plumed.file_name = "plumed.dat" @@ -1383,7 +1492,9 @@ def test_restraints_output_modules(clean_files): assert float(restraint_line[2].split("=")[1]) == 3.1416 assert float(restraint_line[3].split("=")[1]) == 200.0 + # ---------------------------------------------------------------- # # Test Colvar output + # ---------------------------------------------------------------- # colvar = Colvars() colvar.path = "tmp" colvar.file_name = "colvars.dat" diff --git a/paprika/tests/test_utils.py b/paprika/tests/test_utils.py index 6e6b0b90..40e15ce8 100644 --- a/paprika/tests/test_utils.py +++ b/paprika/tests/test_utils.py @@ -28,6 +28,10 @@ def test_mkdirs(): for window in window_list: assert os.path.exists(os.path.join("tmp", "windows", window)) + make_window_dirs(window_list, stash_existing=True, path="tmp") + for window in window_list: + assert os.path.exists(os.path.join("tmp", "windows", window)) + def test_strip_prmtop(): """Test that we can remove items from structures.""" diff --git a/paprika/utils.py b/paprika/utils.py index 1aff72f6..9d9d23f1 100644 --- a/paprika/utils.py +++ b/paprika/utils.py @@ -1,12 +1,13 @@ import logging -import os as os +import os +import re import shutil from datetime import datetime from functools import lru_cache -import numpy as np -import parmed as pmd -import pytraj as pt +import numpy +import parmed +import pytraj from openff.units import unit as openff_unit from parmed.structure import Structure as ParmedStructureClass @@ -102,7 +103,7 @@ def return_parmed_structure(filename): # `parmed` can read both PDBs and # .inpcrd/.prmtop files with the same function call. try: - structure = pmd.load_file(filename) + structure = parmed.load_file(filename) logger.info("Loaded {}...".format(filename)) except IOError: logger.error("Unable to load file: {}".format(filename)) @@ -125,7 +126,7 @@ def index_from_mask(structure, mask, amber_index=False): Returns ------- - indices : int + indices : List[int] Atom index or indices corresponding to the mask. """ @@ -144,7 +145,8 @@ def index_from_mask(structure, mask, amber_index=False): ) # http://parmed.github.io/ParmEd/html/api/parmed/parmed.amber.mask.html?highlight=mask#module-parmed.amber.mask indices = [ - i + index_offset for i in pmd.amber.mask.AmberMask(structure, mask).Selected() + i + index_offset + for i in parmed.amber.mask.AmberMask(structure, mask).Selected() ] logger.debug("There are {} atoms in the mask {} ...".format(len(indices), mask)) return indices @@ -199,8 +201,8 @@ def strip_prmtop(prmtop, mask=":WAT,:Na+,:Cl-"): """ - structure = pt.load_topology(os.path.normpath(prmtop)) - stripped = pt.strip(mask, structure) + structure = pytraj.load_topology(os.path.normpath(prmtop)) + stripped = pytraj.strip(mask, structure) # stripped_name = os.path.join(os.path.splitext(prmtop)[0], '-stripped', os.path.splitext(prmtop)[1]) # stripped.save(filename=stripped_name) # logger.debug('Stripping {} from parameter file and writing {}...'.format(mask, stripped_name)) @@ -345,9 +347,21 @@ def check_unit(variable, base_unit): raise KeyError( "Please make my life easier by either specifying a list of all float or all openff.unit.Quantity." ) - elif isinstance(variable, np.ndarray): + elif isinstance(variable, numpy.ndarray): quantity = openff_unit.Quantity(variable, units=base_unit) else: raise KeyError("``variable`` should be a float or openff.unit.Quantity.") return quantity + + +def multiple_replace(dct, text): + """ + Create a regular expression to do multiple find and replace. + """ + + # Create a regular expression from the dictionary keys + regex = re.compile("(%s)" % "|".join(map(re.escape, dct.keys()))) + + # For each match, look-up corresponding value in dictionary + return regex.sub(lambda mo: dct[mo.string[mo.start() : mo.end()]], text)