Added test relevant to the MDP parameter 'pull'.

ensemble_md/replica_exchange_EE.py

@@ -372,9 +372,7 @@ def set_params(self, analysis):
                 if self.template[f'pull_coord{i+1}_geometry'] == 'distance':
                     if self.template[f'pull_coord{i+1}_start'] == 'yes':
                         self.set_ref_dist.append(True)  # starting from the second iteration, set pull_coord*_init.
-                        if 'pull_nstxout' not in self.template:
-                            self.warnings.append('A non-zero value should be specified for pull_nstxout if pull_coord*_start is set to yes.')  # noqa: E501
-                        if self.template['pull_nstxout'] == 0:
+                        if 'pull_nstxout' not in self.template or self.template['pull_nstxout'] == 0:
                             self.warnings.append('A non-zero value should be specified for pull_nstxout if pull_coord*_start is set to yes.')  # noqa: E501
                     else:
                         self.set_ref_dist.append(False)  # Here we assume that the user know what reference distance to use.  # noqa: E501

ensemble_md/tests/data/expanded_pull.mdp

@@ -0,0 +1,122 @@
+; Run control
+integrator = md-vv
+tinit = 0
+dt = 0.002
+nsteps = 100000000
+nstcomm = 10
+
+; Output control
+nstlog = 1000
+nstcalcenergy = 10
+nstenergy = 1000
+nstxout_compressed = 1000
+
+; Neighborsearching and short-range nonbonded interactions
+nstlist = 10
+ns_type = grid
+pbc = xyz
+rlist = 1.3
+
+; Electrostatics
+coulombtype = PME-switch
+rcoulomb_switch = 0.88
+rcoulomb = 0.9
+
+; van der Waals
+vdw_type = switch
+rvdw_switch = 0.85
+rvdw = 0.9
+
+; Apply long range dispersion corrections for Energy and Pressure
+DispCorr = AllEnerPres
+
+; Spacing for the PME/PPPM FFT grid
+fourierspacing = 0.12
+
+; EWALD/PME/PPPM parameters
+pme_order = 4
+ewald_rtol = 1e-05
+ewald_geometry = 3d
+epsilon_surface = 0
+optimize_fft = yes
+
+; Temperature coupling
+tcoupl = nose-hoover
+nsttcouple = 10
+tc_grps = System
+tau_t = 1.0
+ref_t = 300
+
+; Pressure coupling is on for NPT
+pcoupl = no
+
+gen_vel = yes
+gen_temp = 300
+gen_seed = -1
+
+; options for bonds
+constraints = h-bonds
+
+; Type of constraint algorithm
+constraint_algorithm = shake
+shake_tol = 1e-05
+
+; Free energy calculation
+free_energy = expanded
+calc_lambda_neighbors = -1
+sc_alpha = 0.5
+sc_power = 1
+sc_sigma = 0.5
+couple_moltype = MOL
+couple_lambda0 = vdw-q
+couple_lambda1 = none
+couple_intramol = no
+init_lambda_state = 0
+nstdhdl = 100
+dhdl_print_energy = total
+
+; Seed for Monte Carlo in lambda space
+symmetrized_transition_matrix = no
+nst_transition_matrix = 100000
+; wl-scale                 = 0.8
+; wl-ratio                 = 0.7
+; init-wl-delta            = 10
+
+; expanded ensemble variables
+nstexpanded = 100
+lmc_stats = no
+lmc_move = metropolized-gibbs
+; lmc-weights-equil       = wl-delta
+; weight-equil-wl-delta   = 0.001
+; wl-oneovert             = yes
+
+; lambda-states          = 1      2      3      4      5      6      7      8      9      10     11     12     13     14     15     16     17     18     19     20     21     22     23     24     25     26     27     28     29     30     31     32     33     34     35     36     37     38     39     40
+
+coul_lambdas = 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+vdw_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.45 0.5 0.55 0.6 0.63 0.66 0.69 0.72 0.75 0.78 0.81 0.84 0.88 0.92 1.0
+
+; PULL CODE
+pull = yes
+pull_ngroups = 2
+pull_ncoords = 1
+pull_group1_name = HOS
+pull_group2_name = MOL
+pull_pbc_ref_prev_step_com = yes
+
+pull_coord1_groups = 1 2
+pull_coord1_type = umbrella
+pull_coord1_geometry = distance
+pull_coord1_dim = Y Y Y
+pull_coord1_origin = 0.0 0.0 0.0
+pull_coord1_vec = 0.0 0.0 0.0
+pull_coord1_start = yes
+pull_coord1_init = 0
+pull_coord1_rate = 0
+pull_coord1_k = 0
+pull_coord1_kB = 1000
+pull_nstfout = 400000
+pull_nstxout = 1000
+pull_print_ref_value = yes
+
+restraint_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+init_lambda_weights = 0.0 57.88597 112.71883 163.84425 210.48097 253.80261 294.79849 333.90408 370.82669 406.02515 438.53116 468.53751 496.24649 521.58417 544.57404 565.26697 583.7337 599.60651 613.43958 624.70471 633.95947 638.29785 642.44977 646.33551 649.91626 651.54779 652.93359 654.13263 654.94073 655.13086 655.07239 654.66443 653.68683 652.32123 650.72308 649.2381 647.94586 646.599 645.52063 643.99133

ensemble_md/tests/test_replica_exchange_EE.py

@@ -186,6 +186,31 @@ def test_set_params_error(self, params_dict):
         params_dict['mdp_args'] = {'nstdhdl': [20, 10, 10, 10], 'nstexpanded': [10, 50, 10, 10]}
         with pytest.raises(ParameterError, match='In REXEE, the parameter "nstdhdl" must be a factor of the parameter "nstexpanded", or the calculation of acceptance ratios may be wrong.'):  # noqa: E501
             get_REXEE_instance(params_dict)
+        params_dict['mdp_args'] = None  # set back to a normal value
+
+        # 13. Test the parameter 'pull'
+        params_dict['nst_sim'] = 2000
+        mdp = gmx_parser.MDP(os.path.join(input_path, "expanded_pull.mdp"))
+        params_dict['mdp'] = 'ensemble_md/tests/data/expanded_test.mdp'
+
+        mdp['pull_coord1_geometry'] = 'direction'
+        mdp.write(params_dict['mdp'])
+        get_REXEE_instance(params_dict)
+
+        mdp['pull_coord1_geometry'] = 'distance'  # set back to the original value
+        mdp['pull_coord1_start'] = 'no'
+        mdp.write(params_dict['mdp'])
+        get_REXEE_instance(params_dict)
+
+        mdp['pull_coord1_start'] = 'yes'
+        mdp.write(params_dict['mdp'])
+        mdp['pull_nstxout'] = 0
+        mdp.write(params_dict['mdp'])
+        REXEE_pull = get_REXEE_instance(params_dict)
+        warning = 'A non-zero value should be specified for pull_nstxout if pull_coord*_start is set to yes.'
+        assert warning in REXEE_pull.warnings
+
+        os.remove(os.path.join(input_path, "expanded_test.mdp"))
 
     def test_set_params_warnings(self, params_dict):
         # 1. Non-recognizable parameter in the YAML file