Bug fix

ensemble_md/cli/run_REXEE.py

@@ -311,14 +311,18 @@ def main():
                             gro_2 = f'{REXEE.working_dir}/sim_{swap_list[j][1]}/iteration_{i-1}/confout.gro'
                             print(f'  - {gro_1}\n  - {gro_2}')
 
-                            # Now we rename gro_1 and gro_2 to back them up
-                            gro_1_backup = gro_1.split('.gro')[0] + '_backup.gro'
-                            gro_2_backup = gro_2.split('.gro')[0] + '_backup.gro'
-                            os.rename(gro_1, gro_1_backup)
-                            os.rename(gro_2, gro_2_backup)
-
-                            # Here we input gro_1_backup and gro_2_backup and modify_coords_fn will save the modified gro files as gro_1 and gro_2  # noqa: E501
-                            REXEE.modify_coords_fn(gro_1_backup, gro_2_backup)  # the order should not matter
+                            # Check that swap was not performed before checkpoint was created
+                            if os.path.exists(gro_1.split('.gro')[0] + '_backup.gro') and os.path.exists(gro_2.split('.gro')[0] + '_backup.gro'):  # noqa: E501
+                                print('\nSwap already performed')
+                            else:
+                                # Now we rename gro_1 and gro_2 to back them up
+                                gro_1_backup = gro_1.split('.gro')[0] + '_backup.gro'
+                                gro_2_backup = gro_2.split('.gro')[0] + '_backup.gro'
+                                os.rename(gro_1, gro_1_backup)
+                                os.rename(gro_2, gro_2_backup)
+
+                                # Here we input gro_1_backup and gro_2_backup and modify_coords_fn will save the modified gro files as gro_1 and gro_2  # noqa: E501
+                                REXEE.modify_coords_fn(gro_1_backup, gro_2_backup)  # the order should not matter
                 except Exception:
                     print('\n--------------------------------------------------------------------------\n')
                     print(f'\nAn error occurred on rank 0:\n{traceback.format_exc()}')

ensemble_md/replica_exchange_EE.py

@@ -1608,7 +1608,7 @@ def process_top(self):
                 input_file = gmx_parser.read_top(file_name, self.resname_list[f])
 
                 # Determine the atom names corresponding to the atom numbers
-                start_line, atom_name, state = coordinate_swap.get_names(input_file)
+                start_line, atom_name, atom_num, state = coordinate_swap.get_names(input_file, self.resname_list[f])
 
                 # Determine the connectivity of all atoms
                 connect_1, connect_2, state_1, state_2 = [], [], [], []  # Atom 1 and atom 2 which are connected and which state they are dummy atoms  # noqa: E501
@@ -1620,10 +1620,14 @@ def process_top(self):
                         break
                     while '' in line_sep:
                         line_sep.remove('')
-                    connect_1.append(atom_name[int(line_sep[0])-1])
-                    connect_2.append(atom_name[int(line_sep[1])-1])
-                    state_1.append(state[int(line_sep[0])-1])
-                    state_2.append(state[int(line_sep[1])-1])
+                    if int(line_sep[0]) in atom_num and int(line_sep[1]) in atom_num:
+                        num_1 = np.where(atom_num == int(line_sep[0]))[0][0]
+                        num_2 = np.where(atom_num == int(line_sep[1]))[0][0]
+                        connect_1.append(atom_name[num_1])
+                        connect_2.append(atom_name[num_2])
+                        state_1.append(state[num_1])
+                        state_2.append(state[num_2])
+
                 df = pd.DataFrame({'Resname': self.resname_list[f], 'Connect 1': connect_1, 'Connect 2': connect_2, 'State 1': state_1, 'State 2': state_2})  # noqa: E501
                 df_top = pd.concat([df_top, df])
             df_top.to_csv('residue_connect.csv')
@@ -1638,9 +1642,9 @@ def process_top(self):
                 lam = {X: int(swap[0][1]), Y: int(swap[1][1])}
                 for A, B in zip([X, Y], [Y, X]):
                     input_A = gmx_parser.read_top(self.top[A], self.resname_list[A])
-                    start_line, A_name, state = coordinate_swap.get_names(input_A)
+                    start_line, A_name, A_num, state = coordinate_swap.get_names(input_A, self.resname_list[A])
                     input_B = gmx_parser.read_top(self.top[B], self.resname_list[B])
-                    start_line, B_name, state = coordinate_swap.get_names(input_B)
+                    start_line, B_name, B_num, state = coordinate_swap.get_names(input_B, self.resname_list[B])
 
                     A_only = [x for x in A_name if x not in B_name]
                     B_only = [x for x in B_name if x not in A_name]

ensemble_md/tests/data/coord_swap/output_D.gro

@@ -21,9 +21,9 @@ Protein in water
     1J2K    H11   19   4.7855296   5.1634459   2.6281769   0.000   0.000   0.000
     1J2K    H12   20   4.8469658   5.2264252   2.4736989   0.000   0.000   0.000
     1J2K    H13   21   4.9016318   5.2935543   2.6352797   0.000   0.000   0.000
-    1J2K   HV14   21   5.0237503   5.3428926   2.6540799   0.000   0.000   0.000
-    1J2K   HV15   22   4.8513546   5.4049191   2.6414304   0.000   0.000   0.000
-    1J2K   HV16   23   4.9280410   5.2066817   2.8707309   0.000   0.000   0.000
-    1J2K   HV17   24   4.9172225   5.3824968   2.8867180   0.000   0.000   0.000
-    1J2K   HV18   25   4.7816763   5.2949572   2.8441288   0.000   0.000   0.000
+    1J2K   HV14   22   5.0237503   5.3428926   2.6540799   0.000   0.000   0.000
+    1J2K   HV15   23   4.8513546   5.4049191   2.6414304   0.000   0.000   0.000
+    1J2K   HV16   24   4.9280410   5.2066817   2.8707309   0.000   0.000   0.000
+    1J2K   HV17   25   4.9172225   5.3824968   2.8867180   0.000   0.000   0.000
+    1J2K   HV18   26   4.7816763   5.2949572   2.8441288   0.000   0.000   0.000
    6.81420   6.81420   4.81837   0.00000   0.00000   0.00000   0.00000   3.40710   3.40710

ensemble_md/tests/test_coordinate_swap.py

@@ -346,6 +346,7 @@ def test_swap_name():
 
 def test_get_names():
     top_files = ['A-B.itp', 'B-C.itp', 'C-D.itp', 'D-E.itp', 'E-F.itp']
+    resnames = ['A2B', 'B2C', 'C2D', 'D2E', 'E2F']
 
     start_lines = [26, 29, 33, 32, 36]
     names = [['S1', 'C2', 'N3', 'C4', 'C5', 'C6', 'H1', 'H2', 'H3', 'H4', 'H17', 'DC7', 'HV5', 'HV6', 'HV7'], ['S1', 'C2', 'N3', 'C4', 'C5', 'C6', 'C7', 'H1', 'H2', 'H3', 'H4', 'H5', 'H6', 'H7', 'DC8', 'HV8', 'HV9', 'HV10'], ['S1', 'C2', 'N3', 'C4', 'C5', 'C6', 'C7', 'C8', 'H1', 'H2', 'H3', 'H4', 'H6', 'H7', 'H8', 'H9', 'H10', 'DC9', 'HV5', 'HV11', 'HV12', 'HV13'], ['S1', 'C2', 'N3', 'C4', 'C5', 'C6', 'C7', 'C9', 'H1', 'H2', 'H3', 'H5', 'H6', 'H7', 'H11', 'H12', 'H13', 'DC8', 'HV8', 'HV9', 'HV10'], ['S1', 'C2', 'N3', 'C4', 'C5', 'C6', 'C7', 'C8', 'C9', 'H1', 'H2', 'H3', 'H6', 'H7', 'H8', 'H9', 'H10', 'H11', 'H12', 'H13', 'DC10', 'HV4', 'HV14', 'HV15', 'HV16']]  # noqa: E501
@@ -357,7 +358,7 @@ def test_get_names():
                      [-1, -1, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 0, 0, 0, 0, 0]]
     for i, top_file in enumerate(top_files):
         top = open(f'{input_path}/coord_swap/{top_file}', 'r').readlines()
-        test_start_line, test_names, test_lambda_states = coordinate_swap.get_names(top)
+        test_start_line, test_names, test_nums, test_lambda_states = coordinate_swap.get_names(top, resnames[i])
         assert test_start_line == start_lines[i]
         assert test_names == names[i]
         assert test_lambda_states == lambda_states[i]

ensemble_md/tests/test_replica_exchange_EE.py

@@ -460,9 +460,9 @@ def test_print_params(self, capfd, params_dict):
         L += "Additional runtime arguments: None\n"
         L += "External modules for coordinate manipulation: None\n"
         L += "MDP parameters differing across replicas: None\n"
-        L += "Alchemical ranges of each replica in REXEE:\n  - Replica 0: States [0, 1, 2, 3, 4, 5]\n"
-        L += "  - Replica 1: States [1, 2, 3, 4, 5, 6]\n  - Replica 2: States [2, 3, 4, 5, 6, 7]\n"
-        L += "  - Replica 3: States [3, 4, 5, 6, 7, 8]\n"
+        L += f"Alchemical ranges of each replica in REXEE:\n  - Replica 0: States {[0, 1, 2, 3, 4, 5]}\n"
+        L += f"  - Replica 1: States {[1, 2, 3, 4, 5, 6]}\n  - Replica 2: States {[2, 3, 4, 5, 6, 7]}\n"
+        L += f"  - Replica 3: States {[3, 4, 5, 6, 7, 8]}\n"
         assert out == L
 
         REXEE.reformatted_mdp = True  # Just to test the case where REXEE.reformatted_mdp is True