From e70770cf53d5bf8d289b44b51074b7150c4748cb Mon Sep 17 00:00:00 2001 From: Wei-Tse Hsu Date: Wed, 11 Sep 2024 18:16:11 +0800 Subject: [PATCH 1/3] Added simple unit tests for get_dimensions and compute_angle in coordinate_swap --- ensemble_md/tests/test_coordinate_swap.py | 55 +++++++++++++++++++++++ ensemble_md/utils/coordinate_swap.py | 4 +- 2 files changed, 57 insertions(+), 2 deletions(-) create mode 100644 ensemble_md/tests/test_coordinate_swap.py diff --git a/ensemble_md/tests/test_coordinate_swap.py b/ensemble_md/tests/test_coordinate_swap.py new file mode 100644 index 0000000..2bd7540 --- /dev/null +++ b/ensemble_md/tests/test_coordinate_swap.py @@ -0,0 +1,55 @@ +#################################################################### +# # +# ensemble_md, # +# a python package for running GROMACS simulation ensembles # +# # +# Written by Wei-Tse Hsu # +# Copyright (c) 2022 University of Colorado Boulder # +# # +#################################################################### +""" +Unit tests for the module test_coordainte_swap.py. +""" +import os +import numpy as np +from ensemble_md.utils import coordinate_swap + +current_path = os.path.dirname(os.path.abspath(__file__)) +input_path = os.path.join(current_path, "data") + + +def test_get_dimensions(): + gro = os.path.join(input_path, 'sys.gro') + f = open(gro, 'r') + lines = f.readlines() + f.close() + vec = coordinate_swap.get_dimensions(lines) + assert vec == [3.32017, 3.32017, 2.34772, 0.00000, 0.00000, 0.00000, 0.00000, 1.66009, 1.66009] + + # Write a flat file with cubic box dimensions + f = open('test.gro', 'w') + f.write('test\n') + f.write(' 1.00000 2.00000 3.00000\n') + f.close() + + f = open('test.gro', 'r') + lines = f.readlines() + f.close() + vec = coordinate_swap.get_dimensions(lines) + assert vec == [1.0, 2.0, 3.0] + + os.remove('test.gro') + + +def test_compute_angle(): + coords_1 = [ + np.array([0.0, 0.0, 0.0]), + np.array([1.0, 0.0, 0.0]), + np.array([0.0, 1.0, 0.0]) + ] + coords_2 = coords_1[-1::-1] + coords_3 = [coords_1[1], coords_1[0], coords_1[2]] + + assert np.isclose(coordinate_swap.compute_angle(coords_1), np.sqrt(2) / 2) + assert np.isclose(coordinate_swap.compute_angle(coords_2), np.sqrt(2) / 2) + assert np.isclose(coordinate_swap.compute_angle(coords_3), 0) diff --git a/ensemble_md/utils/coordinate_swap.py b/ensemble_md/utils/coordinate_swap.py index e70f841..a6b0a11 100644 --- a/ensemble_md/utils/coordinate_swap.py +++ b/ensemble_md/utils/coordinate_swap.py @@ -977,8 +977,8 @@ def compute_angle(coords): Parameters ---------- - coords : numpy.ndarray - Four points which define two vectors to compute the angle between. + coords : list + A list of numpy arrays containing the XYZ coordinates of 3 points, for which the angle 1-2-3 is to be computed. Returns ------- From c5cbd33fc21c4766c7194d46e27f688ca142bf65 Mon Sep 17 00:00:00 2001 From: Wei-Tse Hsu Date: Wed, 11 Sep 2024 18:36:15 +0800 Subject: [PATCH 2/3] Tweaked the docstrings of the new functions in replica_exchange_EE.py --- ensemble_md/replica_exchange_EE.py | 36 +++++++++--------------------- 1 file changed, 11 insertions(+), 25 deletions(-) diff --git a/ensemble_md/replica_exchange_EE.py b/ensemble_md/replica_exchange_EE.py index fb63993..91ee776 100644 --- a/ensemble_md/replica_exchange_EE.py +++ b/ensemble_md/replica_exchange_EE.py @@ -19,6 +19,8 @@ import warnings import importlib import subprocess +import mdtraj as md +import pandas as pd import numpy as np from mpi4py import MPI from itertools import combinations @@ -1511,24 +1513,17 @@ def run_REXEE(self, n, swap_pattern=None): # rank 3 that has not been generated, which will lead to an I/O error. comm.barrier() - def default_coords_fn(self, molA_file_name: str, molB_file_name: str): + def default_coords_fn(self, molA_file_name, molB_file_name): """ - Swap coordinates between two GRO files + Swaps coordinates between two GRO files. Parameters ---------- molA_file_name : str - GRO file name for the moleucle to be swapped + GRO file name for the moleucle to be swapped. molB_file_name : str - GRO file name for the other moleucle to be swapped - Return - ------ - None + GRO file name for the other moleucle to be swapped. """ - # Step 1: Load necessary files - import mdtraj as md - import pandas as pd - # Determine name for transformed residue molA_dir = molA_file_name.rsplit('/', 1)[0] + '/' molB_dir = molB_file_name.rsplit('/', 1)[0] + '/' @@ -1541,7 +1536,7 @@ def default_coords_fn(self, molA_file_name: str, molB_file_name: str): connection_map = pd.read_csv('residue_connect.csv') swap_map = pd.read_csv('residue_swap_map.csv') - # Step 2: Read the GRO input coordinate files and open temporary Output files + # Step 1: Read the GRO input coordinate files and open temporary Output files molA_file = open(molA_file_name, 'r').readlines() # open input file molB_new_file_name = 'B_hybrid_swap.gro' molB_new = open(molB_new_file_name, 'w') @@ -1549,19 +1544,19 @@ def default_coords_fn(self, molA_file_name: str, molB_file_name: str): molA_new_file_name = 'A_hybrid_swap.gro' molA_new = open(molA_new_file_name, 'w') - # Step 3: Determine atoms for alignment and swapping + # Step 2: Determine atoms for alignment and swapping nameA = coordinate_swap.identify_res(molA.topology, swap_map['Swap A'].to_list() + swap_map['Swap B'].to_list()) # noqa: E501 nameB = coordinate_swap.identify_res(molB.topology, swap_map['Swap A'].to_list() + swap_map['Swap B'].to_list()) # noqa: E501 df_atom_swap = coordinate_swap.find_common(molA_file, molB_file, nameA, nameB) - # Step 4: Fix break if present for solvated systems only + # Step 3: Fix break if present for solvated systems only if len(molA.topology.select('water')) != 0: A_dimensions = coordinate_swap.get_dimensions(molA_file) B_dimensions = coordinate_swap.get_dimensions(molB_file) molA = coordinate_swap.fix_break(molA, nameA, A_dimensions, connection_map[connection_map['Resname'] == nameA]) # noqa: E501 molB = coordinate_swap.fix_break(molB, nameB, B_dimensions, connection_map[connection_map['Resname'] == nameB]) # noqa: E501 - # Step 5: Determine coordinates of atoms which need to be reconstructed as we swap coordinates between molecules # noqa: E501 + # Step 4: Determine coordinates of atoms which need to be reconstructed as we swap coordinates between molecules # noqa: E501 miss_B = df_atom_swap[(df_atom_swap['Swap'] == 'B2A') & (df_atom_swap['Direction'] == 'miss')]['Name'].to_list() # noqa: E501 miss_A = df_atom_swap[(df_atom_swap['Swap'] == 'A2B') & (df_atom_swap['Direction'] == 'miss')]['Name'].to_list() # noqa: E501 if len(miss_B) != 0: @@ -1588,18 +1583,9 @@ def default_coords_fn(self, molA_file_name: str, molB_file_name: str): def process_top(self): """ - Process the input topologies in order to determine the atoms for alignment in the default GRO swapping + Processes the input topologies in order to determine the atoms for alignment in the default GRO swapping function. Output as csv files to prevent needing to re-run this step. - - Parameters - ---------- - None - Return - ------ - None """ - import pandas as pd - if not os.path.exists('residue_connect.csv'): df_top = pd.DataFrame() for f, file_name in enumerate(self.top): From b6d5fe146fe97a10763bccbbdbe21736d2640db0 Mon Sep 17 00:00:00 2001 From: Wei-Tse Hsu Date: Thu, 12 Sep 2024 00:54:50 +0800 Subject: [PATCH 3/3] Fixed the unit test for compute_angle --- ensemble_md/tests/test_coordinate_swap.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ensemble_md/tests/test_coordinate_swap.py b/ensemble_md/tests/test_coordinate_swap.py index 2bd7540..aa6d677 100644 --- a/ensemble_md/tests/test_coordinate_swap.py +++ b/ensemble_md/tests/test_coordinate_swap.py @@ -50,6 +50,6 @@ def test_compute_angle(): coords_2 = coords_1[-1::-1] coords_3 = [coords_1[1], coords_1[0], coords_1[2]] - assert np.isclose(coordinate_swap.compute_angle(coords_1), np.sqrt(2) / 2) - assert np.isclose(coordinate_swap.compute_angle(coords_2), np.sqrt(2) / 2) - assert np.isclose(coordinate_swap.compute_angle(coords_3), 0) + assert np.isclose(coordinate_swap.compute_angle(coords_1), np.pi / 4) + assert np.isclose(coordinate_swap.compute_angle(coords_2), np.pi / 4) + assert np.isclose(coordinate_swap.compute_angle(coords_3), np.pi / 2)