diff --git a/.codecov.yml b/.codecov.yml index fc6bb83b..2fc9d7ca 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -13,5 +13,6 @@ comment: flags: null paths: null ignore: - "ensemble_md.exceptions.py" - "ensemble_md.utils.py" \ No newline at end of file + - "ensemble_md.exceptions.py" + - "ensemble_md.utils.py" + - "ensemble_md.ensemble_md.py" \ No newline at end of file diff --git a/docs/theory_implementation.rst b/docs/theory_implementation.rst index a55b9b32..b57745ef 100644 --- a/docs/theory_implementation.rst +++ b/docs/theory_implementation.rst @@ -134,10 +134,10 @@ Step 1: Set up parameters To run an ensemble of expanded ensemble in GROMACS using :code:`ensemble_md`, one at least needs to following four files: -- One GRO file of the system of interest -- One TOP file of the system of interest -- One MDP template for customizing different MDP files for different replicas. -- One YAML file that specify the EEXE-relevant parameters. +* One GRO file of the system of interest +* One TOP file of the system of interest +* One MDP template for customizing different MDP files for different replicas. +* One YAML file that specify the EEXE-relevant parameters. Notably, here we are assuming that all replicas start from the same configuration represented by the single GRO file, but the user should also be able to use the methods defined in @@ -151,23 +151,23 @@ Importantly, to instantiate the class :class:`.EnsembleEXE`, the input YAML file In this YAML file, the user needs to specify how the replicas should be set up or interact with each other during the simulation ensemble. Below we decribe the details of these parameters. -- Required parameters +* Required parameters - - :code:`parallel`: Whether the replicas of EEXE should be run in parallel or not. - - :code:`n_sim`: The number of replica simulations. - - :code:`n_iterations`: The number of iterations. - - :code:`s`: The shift in the alchemical ranges between adjacent replicas (e.g. :math:`s = 2` if :math:`λ_2 = (2, 3, 4)` and :math:`λ_3 = (4, 5, 6)`. - - :code:`mdp`: The MDP template that has the whole range of :math:`λ` values. + * :code:`parallel`: Whether the replicas of EEXE should be run in parallel or not. + * :code:`n_sim`: The number of replica simulations. + * :code:`n_iterations`: The number of iterations. + * :code:`s`: The shift in the alchemical ranges between adjacent replicas (e.g. :math:`s = 2` if :math:`λ_2 = (2, 3, 4)` and :math:`λ_3 = (4, 5, 6)`. + * :code:`mdp`: The MDP template that has the whole range of :math:`λ` values. -- Optional parameters +* Optional parameters - - :code:`nst_sim`: The number of simulation steps, i.e. exchange frequency. This option assumes replicas with homogeneous simulation lengths. If this option is not specified, the number of steps defined in the template MDP file will be used. - - :code:`mc_scheme`: The method for swapping simulations. Choices include :code:`same-state`/:code:`same_state`, :code:`metropolis`, and :code:`metropolis-eq`/:code:`metropolis_eq`. For more details, please refer to :ref:`doc_mc_schemes`. (Default: :code:`metropolis`) - - :code:`w_scheme`: The method for combining weights. Choices include :code:`None` (unspecified), :code:`exp-avg`/:code:`exp_avg`, and :code:`hist-exp-avg`/:code:`hist_exp_avg`. For more details, please refer to :ref:`doc_w_schemes`. (Default: :code:`hist-exp-avg`) - - :code:`N_cutoff`: The histogram cutoff. Only required if :code:`hist-exp-avg` is used. (Default: 0) - - :code:`n_ex`: The number of swaps to be proposed in one attempt. This works basically the same as :code:`-nex` flag in GROMACS. A recommended value is :math:`N^3`, where :math:`N` is the number of replicas. If `n_ex` is unspecified or specified as 0, neighboring swapping will be carried out. For more details, please refer to :ref:`doc_swap_basics`. (Default: 0) - - :code:`outfile`: The output file for logging how replicas interact with each other. - - :code:`verbose`: Whether a verbse log is wanted. + * :code:`nst_sim`: The number of simulation steps, i.e. exchange frequency. This option assumes replicas with homogeneous simulation lengths. If this option is not specified, the number of steps defined in the template MDP file will be used. + * :code:`mc_scheme`: The method for swapping simulations. Choices include :code:`same-state`/:code:`same_state`, :code:`metropolis`, and :code:`metropolis-eq`/:code:`metropolis_eq`. For more details, please refer to :ref:`doc_mc_schemes`. (Default: :code:`metropolis`) + * :code:`w_scheme`: The method for combining weights. Choices include :code:`None` (unspecified), :code:`exp-avg`/:code:`exp_avg`, and :code:`hist-exp-avg`/:code:`hist_exp_avg`. For more details, please refer to :ref:`doc_w_schemes`. (Default: :code:`hist-exp-avg`) + * :code:`N_cutoff`: The histogram cutoff. Only required if :code:`hist-exp-avg` is used. (Default: 0) + * :code:`n_ex`: The number of swaps to be proposed in one attempt. This works basically the same as :code:`-nex` flag in GROMACS. A recommended value is :math:`N^3`, where :math:`N` is the number of replicas. If `n_ex` is unspecified or specified as 0, neighboring swapping will be carried out. For more details, please refer to :ref:`doc_swap_basics`. (Default: 0) + * :code:`outfile`: The output file for logging how replicas interact with each other. + * :code:`verbose`: Whether a verbse log is wanted. Step 2: Run the 1st iteration ----------------------------- @@ -186,8 +186,8 @@ Step 3-1: Extract the final status of the previous iteration To calculate the acceptance ratio and modify the mdp files in later steps, we first need to extract the information of the final status of the previous iteration. Specifically, for all the replica simulations, we need to -- Find the last sampled state and the corresponding lambda values from the DHDL files -- Find the final Wang-Landau incrementors and weights from the LOG files. +* Find the last sampled state and the corresponding lambda values from the DHDL files +* Find the final Wang-Landau incrementors and weights from the LOG files. These two tasks are done by :obj:`.extract_final_dhdl_info` and :obj:`.extract_final_log_info`. @@ -232,8 +232,8 @@ when needed), the user should set up the input files for the next iteration. In status of the previous iteration. This means: -- For each replica, the input configuration for initializing a new iterations should be the output configuraiton of the previous iteration. For example, if the final configurations are represented by :code:`[1, 2, 0, 3]` (returned by :obj:`.get_swapped_configs`), then in the next iteration, replica 0 should be initialized by the output configuration of replica 1 in the previous iteration, while replica 3 can just inherit the output configuration from previous iteration of the same replica. Notably, instead of exchanging the MDP files, we recommend swapping out the coordinate files to exchange replicas. -- For each replica, the MDP file for the new iteration should be the same as the one used in the previous iteartion of the same replica except that parameters like :code:`tinit`, :code:`init-lambda-state`, :code:`init-wl-delta`, and :code:`init-lambda-weights` should be modified to the final values in the previous iteration. This can be done by :class:`.gmx_parser.MDP` and :obj:`.update_MDP`. +* For each replica, the input configuration for initializing a new iterations should be the output configuraiton of the previous iteration. For example, if the final configurations are represented by :code:`[1, 2, 0, 3]` (returned by :obj:`.get_swapped_configs`), then in the next iteration, replica 0 should be initialized by the output configuration of replica 1 in the previous iteration, while replica 3 can just inherit the output configuration from previous iteration of the same replica. Notably, instead of exchanging the MDP files, we recommend swapping out the coordinate files to exchange replicas. +* For each replica, the MDP file for the new iteration should be the same as the one used in the previous iteartion of the same replica except that parameters like :code:`tinit`, :code:`init-lambda-state`, :code:`init-wl-delta`, and :code:`init-lambda-weights` should be modified to the final values in the previous iteration. This can be done by :class:`.gmx_parser.MDP` and :obj:`.update_MDP`. Step 4: Run the new iteration ----------------------------- @@ -474,9 +474,9 @@ Each of these replicas sample 6 alchemical states. There alchemical ranges are d Replicas 1 and 2 have overlap at states 2 to 5 and replicas 2 and 3 have overlap at states 4 to 7. Notably, all three replicas sampled states 4 and 5. Therefore, what we are going to do is -- For states 2 and 3, combine weights across replicas 1 and 2. -- For states 4 and 5, combine weights across replicas 1, 2 and 3. -- For states 6 and 7, combine weights across replicas 2 and 3. +* For states 2 and 3, combine weights across replicas 1 and 2. +* For states 4 and 5, combine weights across replicas 1, 2 and 3. +* For states 6 and 7, combine weights across replicas 2 and 3. However, before we combine the weights, we should make sure the weights of all replicas have the same reference because now the references of the 3 replicas are states 0, 2, and 4, respectively. Therefore could be diff --git a/ensemble_md/ensemble_EXE.py b/ensemble_md/ensemble_EXE.py index a8238248..1de4d046 100644 --- a/ensemble_md/ensemble_EXE.py +++ b/ensemble_md/ensemble_EXE.py @@ -47,6 +47,12 @@ def __init__(self, yml_file): outfile : str The file name of the log file for documenting how different replicas interact during the process. + + Notes + ----- + For easier parsing of the mdp template file, please make sure that at least the following GROMACS mdp + parameters specified using dashes instead of underscores: :code:`ref-t`, :code:`vdw-lambdas`, + :code:`coul-lambdas`, :code:`restraint-lambdas`, and :code:`init-lambda-weights`. """ # Step 0: Set up constants k = 1.380649e-23 @@ -87,33 +93,33 @@ def __init__(self, yml_file): params_int = ['n_sim', 'n_iterations', 's', 'nst_sim', 'N_cutoff', 'n_ex'] # integer parameters for i in params_int: if type(getattr(self, i)) != int: - raise ParameterError(f"The parameter {i} should be an integer.") + raise ParameterError(f"The parameter '{i}' should be an integer.") params_pos = ['n_sim', 'n_iterations', 's', 'nst_sim'] # positive parameters for i in params_pos: if getattr(self, i) <= 0: - raise ParameterError(f"The parameter {i} should be positive.") + raise ParameterError(f"The parameter '{i}' should be positive.") params_non_neg = ['N_cutoff', 'n_ex'] # non-negative parameters for i in params_non_neg: if getattr(self, i) < 0: - raise ParameterError(f"The parameter {i} should be non-negative.") + raise ParameterError(f"The parameter '{i}' should be non-negative.") params_str = ['mdp', 'outfile'] for i in params_str: if type(getattr(self, i)) != str: - raise ParameterError(f"The parameter {i} should be a string.") + raise ParameterError(f"The parameter '{i}' should be a string.") params_bool = ['parallel', 'verbose'] for i in params_bool: if type(getattr(self, i)) != bool: - raise ParameterError(f"The parameter {i} should be a boolean variable.") + raise ParameterError(f"The parameter '{i}' should be a boolean variable.") # Step 4: Read in parameters from the MDP template self.template = gmx_parser.MDP(self.mdp) self.nsteps = self.template["nsteps"] # will be overwritten by self.nst_sim if nst_sim is specified. self.dt = self.template["dt"] # ps - self.temp = self.template["ref_t"] + self.temp = self.template["ref-t"] self.kT = k * NA * self.temp / 1000 # 1 kT in kJ/mol # Total # of states. n_tot = n_sub * n_sim - (n_overlap) * (n_sum - 1), where n_overlap = n_sub - s @@ -179,7 +185,7 @@ def map_lambda2state(self): ( self.template["coul-lambdas"][i], self.template["vdw-lambdas"][i], - self.template["restraint_lambdas"][i], + self.template["restraint-lambdas"][i], ) ] = i else: @@ -190,8 +196,9 @@ def map_lambda2state(self): ) ] = i else: - self.lambda_dict[(self.template["vdw-lambdas"][i])] = i - + self.lambda_dict[(self.template["vdw-lambdas"][i],)] = i + print(self.lambda_dict) + print(self.state_ranges) self.lambda_ranges = [[list(self.lambda_dict.keys())[j] for j in self.state_ranges[i]]for i in range(len(self.state_ranges))] # noqa: E501 def initialize_MDP(self, idx): @@ -380,12 +387,11 @@ def propose_swaps(self, states): print(f"Swappable pairs: {swappables}") - for i in range(n_ex): - try: - swap_list = random.choices(swappables, k=n_ex) - except IndexError: - # In the case that swappables is an empty list, i.e. no swappable pairs. - swap_list = None + try: + swap_list = random.choices(swappables, k=n_ex) + except IndexError: + # In the case that swappables is an empty list, i.e. no swappable pairs. + swap_list = None return swap_list diff --git a/ensemble_md/tests/data/EXE.log b/ensemble_md/tests/data/EXE.log new file mode 100644 index 00000000..50958172 --- /dev/null +++ b/ensemble_md/tests/data/EXE.log @@ -0,0 +1,796 @@ + :-) GROMACS - gmx mdrun, 2021.4-plumed-2.8.0 (-: + + GROMACS is written by: + Andrey Alekseenko Emile Apol Rossen Apostolov + Paul Bauer Herman J.C. Berendsen Par Bjelkmar + Christian Blau Viacheslav Bolnykh Kevin Boyd + Aldert van Buuren Rudi van Drunen Anton Feenstra + Gilles Gouaillardet Alan Gray Gerrit Groenhof + Anca Hamuraru Vincent Hindriksen M. Eric Irrgang + Aleksei Iupinov Christoph Junghans Joe Jordan + Dimitrios Karkoulis Peter Kasson Jiri Kraus + Carsten Kutzner Per Larsson Justin A. Lemkul + Viveca Lindahl Magnus Lundborg Erik Marklund + Pascal Merz Pieter Meulenhoff Teemu Murtola + Szilard Pall Sander Pronk Roland Schulz + Michael Shirts Alexey Shvetsov Alfons Sijbers + Peter Tieleman Jon Vincent Teemu Virolainen + Christian Wennberg Maarten Wolf Artem Zhmurov + and the project leaders: + Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel + +Copyright (c) 1991-2000, University of Groningen, The Netherlands. +Copyright (c) 2001-2019, The GROMACS development team at +Uppsala University, Stockholm University and +the Royal Institute of Technology, Sweden. +check out http://www.gromacs.org for more information. + +GROMACS is free software; you can redistribute it and/or modify it +under the terms of the GNU Lesser General Public License +as published by the Free Software Foundation; either version 2.1 +of the License, or (at your option) any later version. + +GROMACS: gmx mdrun, version 2021.4-plumed-2.8.0 +Executable: /usr/local/gromacs/bin/gmx +Data prefix: /usr/local/gromacs +Working dir: /Users/Wei-TseHsu/Documents/Life_in_CU_Bouler/Research_in_Shirts_Lab/EEXE_experiments/Preliminary_tests/Devop/gmxapi.commandline.cli1_i0 +Process ID: 64740 +Command line: + gmx mdrun -s /Users/Wei-TseHsu/Documents/Life_in_CU_Bouler/Research_in_Shirts_Lab/EEXE_experiments/Preliminary_tests/Devop/sim_0/iteration_0/sys_EE.tpr + +GROMACS version: 2021.4-plumed-2.8.0 +Precision: mixed +Memory model: 64 bit +MPI library: thread_mpi +OpenMP support: disabled +GPU support: disabled +SIMD instructions: AVX2_256 +FFT library: fftw-3.3.10-sse2-avx-avx2-avx2_128 +RDTSCP usage: enabled +TNG support: enabled +Hwloc support: disabled +Tracing support: disabled +C compiler: /usr/bin/gcc AppleClang 12.0.0.12000032 +C compiler flags: -mavx2 -mfma -Wno-unknown-pragmas -Wno-missing-field-initializers -fno-stack-check -fno-stack-check -O3 -DNDEBUG +C++ compiler: /usr/bin/g++ AppleClang 12.0.0.12000032 +C++ compiler flags: -mavx2 -mfma -Wno-unknown-pragmas -Wno-missing-field-initializers -fno-stack-check -fno-stack-check -Weverything -Wno-c++98-compat -Wno-c++98-compat-pedantic -Wno-source-uses-openmp -Wno-c++17-extensions -Wno-documentation-unknown-command -Wno-covered-switch-default -Wno-switch-enum -Wno-extra-semi-stmt -Wno-disabled-macro-expansion -Wno-cast-align -Wno-reserved-id-macro -Wno-global-constructors -Wno-exit-time-destructors -Wno-unused-macros -Wno-weak-vtables -Wno-conditional-uninitialized -Wno-format-nonliteral -Wno-shadow -Wno-cast-qual -Wno-documentation -Wno-used-but-marked-unused -Wno-padded -Wno-float-equal -Wno-old-style-cast -Wno-conversion -Wno-double-promotion -O3 -DNDEBUG + + +Running on 1 node with total 8 cores, 8 logical cores +Hardware detected: + CPU info: + Vendor: Intel + Brand: Intel(R) Core(TM) i5-8257U CPU @ 1.40GHz + Family: 6 Model: 142 Stepping: 10 + Features: aes apic avx avx2 clfsh cmov cx8 cx16 f16c fma htt intel lahf mmx msr nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdrnd rdtscp sse2 sse3 sse4.1 sse4.2 ssse3 tdt x2apic + Hardware topology: Only logical processor count + + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, E. +Lindahl +GROMACS: High performance molecular simulations through multi-level +parallelism from laptops to supercomputers +SoftwareX 1 (2015) pp. 19-25 +-------- -------- --- Thank You --- -------- -------- + + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +S. Páll, M. J. Abraham, C. Kutzner, B. Hess, E. Lindahl +Tackling Exascale Software Challenges in Molecular Dynamics Simulations with +GROMACS +In S. Markidis & E. Laure (Eds.), Solving Software Challenges for Exascale 8759 (2015) pp. 3-27 +-------- -------- --- Thank You --- -------- -------- + + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. +Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess, and E. Lindahl +GROMACS 4.5: a high-throughput and highly parallel open source molecular +simulation toolkit +Bioinformatics 29 (2013) pp. 845-54 +-------- -------- --- Thank You --- -------- -------- + + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +B. Hess and C. Kutzner and D. van der Spoel and E. Lindahl +GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable +molecular simulation +J. Chem. Theory Comput. 4 (2008) pp. 435-447 +-------- -------- --- Thank You --- -------- -------- + + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. +Berendsen +GROMACS: Fast, Flexible and Free +J. Comp. Chem. 26 (2005) pp. 1701-1719 +-------- -------- --- Thank You --- -------- -------- + + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +E. Lindahl and B. Hess and D. van der Spoel +GROMACS 3.0: A package for molecular simulation and trajectory analysis +J. Mol. Mod. 7 (2001) pp. 306-317 +-------- -------- --- Thank You --- -------- -------- + + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +H. J. C. Berendsen, D. van der Spoel and R. van Drunen +GROMACS: A message-passing parallel molecular dynamics implementation +Comp. Phys. Comm. 91 (1995) pp. 43-56 +-------- -------- --- Thank You --- -------- -------- + +Input Parameters: + integrator = md-vv + tinit = 0 + dt = 0.002 + nsteps = 500 + init-step = 0 + simulation-part = 1 + mts = false + comm-mode = Linear + nstcomm = 1 + bd-fric = 0 + ld-seed = -50397473 + emtol = 10 + emstep = 0.01 + niter = 20 + fcstep = 0 + nstcgsteep = 1000 + nbfgscorr = 10 + rtpi = 0.05 + nstxout = 0 + nstvout = 0 + nstfout = 0 + nstlog = 100 + nstcalcenergy = 1 + nstenergy = 1000 + nstxout-compressed = 1000 + compressed-x-precision = 1000 + cutoff-scheme = Verlet + nstlist = 10 + pbc = xyz + periodic-molecules = false + verlet-buffer-tolerance = 0.005 + rlist = 0.9 + coulombtype = PME + coulomb-modifier = Potential-shift + rcoulomb-switch = 0.89 + rcoulomb = 0.9 + epsilon-r = 1 + epsilon-rf = inf + vdw-type = Cut-off + vdw-modifier = Potential-switch + rvdw-switch = 0.85 + rvdw = 0.9 + DispCorr = AllEnerPres + table-extension = 1 + fourierspacing = 0.1 + fourier-nx = 36 + fourier-ny = 36 + fourier-nz = 36 + pme-order = 4 + ewald-rtol = 1e-05 + ewald-rtol-lj = 0.001 + lj-pme-comb-rule = Geometric + ewald-geometry = 0 + epsilon-surface = 0 + tcoupl = V-rescale + nsttcouple = 1 + nh-chain-length = 0 + print-nose-hoover-chain-variables = false + pcoupl = No + pcoupltype = Isotropic + nstpcouple = -1 + tau-p = 1 + compressibility (3x3): + compressibility[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + compressibility[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + compressibility[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + ref-p (3x3): + ref-p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + ref-p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + ref-p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + refcoord-scaling = No + posres-com (3): + posres-com[0]= 0.00000e+00 + posres-com[1]= 0.00000e+00 + posres-com[2]= 0.00000e+00 + posres-comB (3): + posres-comB[0]= 0.00000e+00 + posres-comB[1]= 0.00000e+00 + posres-comB[2]= 0.00000e+00 + QMMM = false +qm-opts: + ngQM = 0 + constraint-algorithm = Lincs + continuation = false + Shake-SOR = false + shake-tol = 0.0001 + lincs-order = 12 + lincs-iter = 2 + lincs-warnangle = 30 + nwall = 0 + wall-type = 9-3 + wall-r-linpot = -1 + wall-atomtype[0] = -1 + wall-atomtype[1] = -1 + wall-density[0] = 0 + wall-density[1] = 0 + wall-ewald-zfac = 3 + pull = false + awh = false + rotation = false + interactiveMD = false + disre = No + disre-weighting = Conservative + disre-mixed = false + dr-fc = 1000 + dr-tau = 0 + nstdisreout = 100 + orire-fc = 0 + orire-tau = 0 + nstorireout = 100 + free-energy = expanded + init-lambda = -1 + init-lambda-state = 0 + delta-lambda = 0 + nstdhdl = 10 + n-lambdas = 6 + separate-dvdl: + fep-lambdas = FALSE + mass-lambdas = FALSE + coul-lambdas = TRUE + vdw-lambdas = TRUE + bonded-lambdas = FALSE + restraint-lambdas = FALSE +temperature-lambdas = FALSE +all-lambdas: + fep-lambdas = 0 0 0 0 0 0 + mass-lambdas = 0 0 0 0 0 0 + coul-lambdas = 0 0.25 0.5 0.75 1 1 + vdw-lambdas = 0 0 0 0 0 0.25 + bonded-lambdas = 0 0 0 0 0 0 + restraint-lambdas = 0 0 0 0 0 0 +temperature-lambdas = 0 0 0 0 0 0 + calc-lambda-neighbors = -1 + dhdl-print-energy = total + sc-alpha = 0.5 + sc-power = 1 + sc-r-power = 6 + sc-sigma = 0.3 + sc-sigma-min = 0.3 + sc-coul = false + dh-hist-size = 0 + dh-hist-spacing = 0.1 + separate-dhdl-file = yes + dhdl-derivatives = yes + nstexpanded = 10 + lmc-stats = wang-landau + lmc-move = metropolized-gibbs + lmc-weights-equil = wl-delta + weight-equil-wl-delta = 0.001 + lmc-seed = 1000 + mc-temperature = 298 + lmc-repeats = 1 + lmc-gibbsdelta = -1 + lmc-forced-nstart = 0 + symmetrized-transition-matrix = true + nst-transition-matrix = 100 + mininum-var-min = 100 + weight-c-range = 0 + wl-scale = 0.8 + wl-ratio = 0.8 + init-wl-delta = 0.5 + wl-oneovert = false + init-lambda-weights (6): + init-lambda-weights[0]= 0.00000e+00 + init-lambda-weights[1]= 0.00000e+00 + init-lambda-weights[2]= 0.00000e+00 + init-lambda-weights[3]= 0.00000e+00 + init-lambda-weights[4]= 0.00000e+00 + init-lambda-weights[5]= 0.00000e+00 + init-weights = false + cos-acceleration = 0 + deform (3x3): + deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} + simulated-tempering = false + swapcoords = no + userint1 = 0 + userint2 = 0 + userint3 = 0 + userint4 = 0 + userreal1 = 0 + userreal2 = 0 + userreal3 = 0 + userreal4 = 0 + applied-forces: + electric-field: + x: + E0 = 0 + omega = 0 + t0 = 0 + sigma = 0 + y: + E0 = 0 + omega = 0 + t0 = 0 + sigma = 0 + z: + E0 = 0 + omega = 0 + t0 = 0 + sigma = 0 + density-guided-simulation: + active = false + group = protein + similarity-measure = inner-product + atom-spreading-weight = unity + force-constant = 1e+09 + gaussian-transform-spreading-width = 0.2 + gaussian-transform-spreading-range-in-multiples-of-width = 4 + reference-density-filename = reference.mrc + nst = 1 + normalize-densities = true + adaptive-force-scaling = false + adaptive-force-scaling-time-constant = 4 + shift-vector = + transformation-matrix = +grpopts: + nrdf: 5043 + ref-t: 298 + tau-t: 0.5 +annealing: No +annealing-npoints: 0 + acc: 0 0 0 + nfreeze: N N N + energygrp-flags[ 0]: 0 + +Changing nstlist from 10 to 40, rlist from 0.9 to 0.999 + + +Initializing Domain Decomposition on 8 ranks +Dynamic load balancing: auto +Using update groups, nr 843, average size 3.0 atoms, max. radius 0.078 nm +Minimum cell size due to atom displacement: 0.246 nm +Initial maximum distances in bonded interactions: + two-body bonded interactions: 0.317 nm, LJC-14 q, atoms 1 4 + multi-body bonded interactions: 0.317 nm, Proper Dih., atoms 1 4 +Minimum cell size due to bonded interactions: 0.348 nm +Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25 +Using 0 separate PME ranks, as there are too few total + ranks for efficient splitting +Optimizing the DD grid for 8 cells with a minimum initial size of 0.436 nm +The maximum allowed number of cells is: X 6 Y 6 Z 5 +Domain decomposition grid 4 x 2 x 1, separate PME ranks 0 +PME domain decomposition: 4 x 2 x 1 +Domain decomposition rank 0, coordinates 0 0 0 + +The initial number of communication pulses is: X 2 Y 1 +The initial domain decomposition cell size is: X 0.68 nm Y 1.36 nm + +The maximum allowed distance for atom groups involved in interactions is: + non-bonded interactions 1.155 nm + two-body bonded interactions (-rdd) 1.155 nm + multi-body bonded interactions (-rdd) 0.678 nm + +When dynamic load balancing gets turned on, these settings will change to: +The maximum number of communication pulses is: X 3 Y 1 +The minimum size for domain decomposition cells is 0.513 nm +The requested allowed shrink of DD cells (option -dds) is: 0.80 +The allowed shrink of domain decomposition cells is: X 0.76 Y 0.85 +The maximum allowed distance for atom groups involved in interactions is: + non-bonded interactions 1.155 nm + two-body bonded interactions (-rdd) 1.155 nm + multi-body bonded interactions (-rdd) 0.513 nm + +Using 8 MPI threads +System total charge, top. A: 0.000 top. B: 0.000 +Will do PME sum in reciprocal space for electrostatic interactions. + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen +A smooth particle mesh Ewald method +J. Chem. Phys. 103 (1995) pp. 8577-8592 +-------- -------- --- Thank You --- -------- -------- + +Using a Gaussian width (1/beta) of 0.288146 nm for Ewald +Potential shift: LJ r^-12: 0.000e+00 r^-6: 0.000e+00, Ewald -1.111e-05 +Initialized non-bonded Ewald tables, spacing: 8.85e-04 size: 2260 + +Generated table with 999 data points for 1-4 COUL. +Tabscale = 500 points/nm +Generated table with 999 data points for 1-4 LJ6. +Tabscale = 500 points/nm +Generated table with 999 data points for 1-4 LJ12. +Tabscale = 500 points/nm +Long Range LJ corr.: 2.7898e-04 2.7407e-07 + + +Using SIMD 4x8 nonbonded short-range kernels + +Using a dual 4x8 pair-list setup updated with dynamic pruning: + outer list: updated every 40 steps, buffer 0.099 nm, rlist 0.999 nm + inner list: updated every 11 steps, buffer 0.002 nm, rlist 0.902 nm +At tolerance 0.005 kJ/mol/ps per atom, equivalent classical 1x1 list would be: + outer list: updated every 40 steps, buffer 0.209 nm, rlist 1.109 nm + inner list: updated every 11 steps, buffer 0.042 nm, rlist 0.942 nm +There are 4 atoms and 2 charges for free energy perturbation +Removing pbc first time + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +S. Miyamoto and P. A. Kollman +SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithms for Rigid +Water Models +J. Comp. Chem. 13 (1992) pp. 952-962 +-------- -------- --- Thank You --- -------- -------- + + +Linking all bonded interactions to atoms + + +Intra-simulation communication will occur every 1 steps. +Initial vector of lambda components:[ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 ] + +++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ +G. Bussi, D. Donadio and M. Parrinello +Canonical sampling through velocity rescaling +J. Chem. Phys. 126 (2007) pp. 014101 +-------- -------- --- Thank You --- -------- -------- + +There are: 2521 Atoms +Atom distribution over 8 domains: av 315 stddev 19 min 294 max 342 + +Constraining the starting coordinates (step 0) +Center of mass motion removal mode is Linear +We have the following groups for center of mass motion removal: + 0: rest +RMS relative constraint deviation after constraining: 0.00e+00 +Initial temperature: 300.024 K + +Started mdrun on rank 0 Sun Aug 7 21:24:32 2022 + + Step Time + 0 0.00000 + + Energies (kJ/mol) + Bond Angle Proper Dih. LJ-14 Coulomb-14 + 2.66745e+00 5.47638e+00 4.47967e+00 2.97400e+00 -1.46179e+01 + LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential + 5.22675e+03 -2.14179e+02 -3.83612e+04 3.35401e+02 -3.30122e+04 + Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar) + 6.28997e+03 -2.67223e+04 -2.67223e+04 3.00024e+02 -1.37273e+02 + Pressure (bar) dVcoul/dl dVvdw/dl + 2.46153e+02 6.48052e+00 2.33270e+01 + + +DD step 39 load imb.: force 21.8% + +step 80 Turning on dynamic load balancing, because the performance loss due to load imbalance is 2.4 %. + Step Time + 100 0.20000 + + MC-lambda information + Wang-Landau incrementor is: 0.5 + N CoulL VdwL Count G(in kT) dG(in kT) + 1 0.000 0.000 6 0.00000 2.00000 << + 2 0.250 0.000 2 2.00000 0.50000 + 3 0.500 0.000 1 2.50000 0.00000 + 4 0.750 0.000 1 2.50000 0.50000 + 5 1.000 0.000 0 3.00000 0.00000 + 6 1.000 0.250 0 3.00000 0.00000 + + Transition Matrix + 1 2 3 4 5 6 + 0.81909299 0.27110553 0.14064617 0.09602036 0.00064359 0.00022395 1 + 0.27110553 0.10953421 0.15520282 0.12907360 0.03806938 0.00291879 2 + 0.14064617 0.15520282 0.00007268 0.04758952 0.00022461 0.00005783 3 + 0.09602036 0.12907360 0.04758952 0.00008816 0.00483189 0.00039835 4 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 5 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 6 + Empirical Transition Matrix + 1 2 3 4 5 6 + 0.66666669 0.25000000 0.28571430 0.00000000 0.00000000 0.00000000 1 + 0.25000000 0.00000000 0.00000000 0.66666669 0.00000000 0.00000000 2 + 0.28571430 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 3 + 0.00000000 0.66666669 0.00000000 0.00000000 0.00000000 0.00000000 4 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 5 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 6 + Energies (kJ/mol) + Bond Angle Proper Dih. LJ-14 Coulomb-14 + 5.66995e+00 6.32659e+00 5.31139e+00 2.00192e+00 -1.43375e+01 + LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential + 5.10216e+03 -2.14179e+02 -3.84731e+04 3.12680e+02 -3.32675e+04 + Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar) + 6.33264e+03 -2.69348e+04 -2.67338e+04 3.02059e+02 -1.37273e+02 + Pressure (bar) dVcoul/dl dVvdw/dl + -1.38783e+02 3.93696e+01 -3.56912e+01 + + +DD step 199 vol min/aver 0.940 load imb.: force 19.1% + Step Time + 200 0.40000 + + MC-lambda information + Wang-Landau incrementor is: 0.5 + N CoulL VdwL Count G(in kT) dG(in kT) + 1 0.000 0.000 10 0.00000 2.00000 + 2 0.250 0.000 6 2.00000 2.00000 << + 3 0.500 0.000 2 4.00000 0.50000 + 4 0.750 0.000 1 4.50000 0.00000 + 5 1.000 0.000 1 4.50000 0.50000 + 6 1.000 0.250 0 5.00000 0.00000 + + Transition Matrix + 1 2 3 4 5 6 + 0.55030239 0.36768454 0.15483068 0.07404882 0.00446463 0.00093804 1 + 0.36768454 0.25592616 0.27140886 0.10001598 0.06546144 0.00258266 2 + 0.15483068 0.27140886 0.04135029 0.06092318 0.11434975 0.00359765 3 + 0.07404882 0.10001598 0.06092318 0.00008816 0.16502382 0.00039835 4 + 0.00446463 0.06546144 0.11434975 0.16502382 0.03636860 0.01974302 5 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 6 + Empirical Transition Matrix + 1 2 3 4 5 6 + 0.50000000 0.37500000 0.25000000 0.00000000 0.00000000 0.00000000 1 + 0.37500000 0.33333334 0.00000000 0.28571430 0.14285715 0.00000000 2 + 0.25000000 0.00000000 0.00000000 0.00000000 0.33333334 0.00000000 3 + 0.00000000 0.28571430 0.00000000 0.00000000 0.00000000 0.00000000 4 + 0.00000000 0.14285715 0.33333334 0.00000000 0.00000000 0.00000000 5 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 6 + Energies (kJ/mol) + Bond Angle Proper Dih. LJ-14 Coulomb-14 + 1.06594e+00 6.01026e+00 1.69249e+00 1.79315e+00 -1.42662e+01 + LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential + 5.13488e+03 -2.14179e+02 -3.85906e+04 3.20424e+02 -3.33532e+04 + Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar) + 6.26268e+03 -2.70905e+04 -2.67188e+04 2.98722e+02 -1.37273e+02 + Pressure (bar) dVcoul/dl dVvdw/dl + -1.64534e+02 2.25177e+01 1.17826e+01 + + Step Time + 300 0.60000 + + MC-lambda information + Wang-Landau incrementor is: 0.5 + N CoulL VdwL Count G(in kT) dG(in kT) + 1 0.000 0.000 12 0.00000 1.50000 + 2 0.250 0.000 9 1.50000 2.50000 << + 3 0.500 0.000 4 4.00000 1.00000 + 4 0.750 0.000 2 5.00000 -0.50000 + 5 1.000 0.000 3 4.50000 1.50000 + 6 1.000 0.250 0 6.00000 0.00000 + + Transition Matrix + 1 2 3 4 5 6 + 0.50359410 0.30349964 0.16633579 0.07859923 0.01222900 0.00350482 1 + 0.30349964 0.17466077 0.27996165 0.15325463 0.08172923 0.02491526 2 + 0.16633579 0.27996165 0.11296166 0.16483080 0.20322983 0.05683871 3 + 0.07859923 0.15325463 0.16483080 0.15837476 0.23805423 0.06371155 4 + 0.01222900 0.08172923 0.20322983 0.23805423 0.01809935 0.23811643 5 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 6 + Empirical Transition Matrix + 1 2 3 4 5 6 + 0.50000000 0.38095239 0.18750000 0.00000000 0.00000000 0.00000000 1 + 0.38095239 0.22222222 0.15384616 0.27272728 0.08333334 0.00000000 2 + 0.18750000 0.15384616 0.00000000 0.00000000 0.42857143 0.00000000 3 + 0.00000000 0.27272728 0.00000000 0.00000000 0.40000001 0.00000000 4 + 0.00000000 0.08333334 0.42857143 0.40000001 0.00000000 0.00000000 5 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 6 + Energies (kJ/mol) + Bond Angle Proper Dih. LJ-14 Coulomb-14 + 2.24621e+00 2.63433e-02 8.11236e-01 8.82040e+00 -1.55741e+01 + LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential + 5.03718e+03 -2.14179e+02 -3.83012e+04 3.38728e+02 -3.31432e+04 + Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar) + 6.29869e+03 -2.68445e+04 -2.67215e+04 3.00440e+02 -1.37273e+02 + Pressure (bar) dVcoul/dl dVvdw/dl + -1.58800e+02 5.78631e-01 1.00147e+01 + + +DD step 399 vol min/aver 0.868 load imb.: force 8.4% + Step Time + 400 0.80000 + + MC-lambda information + Wang-Landau incrementor is: 0.5 + N CoulL VdwL Count G(in kT) dG(in kT) + 1 0.000 0.000 12 0.00000 1.50000 + 2 0.250 0.000 9 1.50000 1.50000 + 3 0.500 0.000 6 3.00000 0.00000 << + 4 0.750 0.000 6 3.00000 0.50000 + 5 1.000 0.000 5 3.50000 1.50000 + 6 1.000 0.250 2 5.00000 0.00000 + + Transition Matrix + 1 2 3 4 5 6 + 0.50359410 0.30349964 0.16638361 0.06404528 0.01149559 0.00359449 1 + 0.30349964 0.17466077 0.26737976 0.12295914 0.07353259 0.02331973 2 + 0.16638361 0.26737976 0.16847293 0.21531832 0.21830393 0.09132640 3 + 0.06404528 0.12295914 0.21531832 0.09865310 0.29120502 0.16403839 4 + 0.01149559 0.07353259 0.21830393 0.29120502 0.08283626 0.30287614 5 + 0.00359449 0.02331973 0.09132640 0.16403839 0.30287614 0.05424957 6 + Empirical Transition Matrix + 1 2 3 4 5 6 + 0.50000000 0.38095239 0.16666667 0.00000000 0.00000000 0.00000000 1 + 0.38095239 0.22222222 0.20000000 0.20000000 0.07142857 0.00000000 2 + 0.16666667 0.20000000 0.00000000 0.16666667 0.36363637 0.00000000 3 + 0.00000000 0.20000000 0.16666667 0.16666667 0.36363637 0.12500000 4 + 0.00000000 0.07142857 0.36363637 0.36363637 0.00000000 0.14285715 5 + 0.00000000 0.00000000 0.00000000 0.12500000 0.14285715 0.50000000 6 + Energies (kJ/mol) + Bond Angle Proper Dih. LJ-14 Coulomb-14 + 4.14196e+00 1.12954e-01 4.51524e+00 5.23232e+00 -1.50821e+01 + LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential + 4.85381e+03 -2.14179e+02 -3.81588e+04 3.39099e+02 -3.31811e+04 + Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar) + 6.34858e+03 -2.68326e+04 -2.67208e+04 3.02819e+02 -1.37273e+02 + Pressure (bar) dVcoul/dl dVvdw/dl + -6.41110e+02 1.47531e+01 6.64636e+00 + + +DD step 479 vol min/aver 0.908 load imb.: force 10.7% + Step Time + 500 1.00000 + +Writing checkpoint, step 500 at Sun Aug 7 21:24:35 2022 + + + MC-lambda information + Wang-Landau incrementor is: 0.5 + N CoulL VdwL Count G(in kT) dG(in kT) + 1 0.000 0.000 12 0.00000 0.50000 + 2 0.250 0.000 11 0.50000 1.50000 + 3 0.500 0.000 8 2.00000 1.00000 + 4 0.750 0.000 6 3.00000 -1.00000 + 5 1.000 0.000 8 2.00000 1.50000 + 6 1.000 0.250 5 3.50000 0.00000 << + + Transition Matrix + 1 2 3 4 5 6 + 0.50359410 0.27816516 0.15014392 0.06404528 0.01056806 0.00325544 1 + 0.27816516 0.14302190 0.22873016 0.12587893 0.07826510 0.07551135 2 + 0.15014392 0.22873016 0.12687601 0.22781774 0.19645289 0.14166102 3 + 0.06404528 0.12587893 0.22781774 0.09865310 0.27133784 0.18215007 4 + 0.01056806 0.07826510 0.19645289 0.27133784 0.06328789 0.36931726 5 + 0.00325544 0.07551135 0.14166102 0.18215007 0.36931726 0.29827550 6 + Empirical Transition Matrix + 1 2 3 4 5 6 + 0.50000000 0.34782609 0.15000001 0.00000000 0.00000000 0.00000000 1 + 0.34782609 0.18181819 0.15789473 0.17647059 0.10526316 0.12500000 2 + 0.15000001 0.15789473 0.00000000 0.14285715 0.43750000 0.07692308 3 + 0.00000000 0.17647059 0.14285715 0.16666667 0.28571430 0.18181819 4 + 0.00000000 0.10526316 0.43750000 0.28571430 0.00000000 0.23076923 5 + 0.00000000 0.12500000 0.07692308 0.18181819 0.23076923 0.20000000 6 + Energies (kJ/mol) + Bond Angle Proper Dih. LJ-14 Coulomb-14 + 5.44325e+00 2.28126e+00 8.23533e+00 2.04241e+00 -1.43507e+01 + LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential + 5.04307e+03 -2.13560e+02 -3.84729e+04 3.23321e+02 -3.33164e+04 + Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar) + 6.27854e+03 -2.70379e+04 -2.67193e+04 2.99478e+02 -1.36876e+02 + Pressure (bar) dVcoul/dl dVvdw/dl + -4.02436e+02 -5.16640e+00 2.10354e+01 + + +Energy conservation over simulation part #1 of length 1 ns, time 0 to 1 ns + Conserved energy drift: 1.17e-03 kJ/mol/ps per atom + + + <====== ############### ==> + <==== A V E R A G E S ====> + <== ############### ======> + + Statistics over 501 steps using 501 frames + + Energies (kJ/mol) + Bond Angle Proper Dih. LJ-14 Coulomb-14 + 2.69082e+00 1.82488e+00 4.28850e+00 4.59528e+00 -1.48978e+01 + LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential + 5.00771e+03 -2.14117e+02 -3.83758e+04 3.33078e+02 -3.32506e+04 + Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar) + 6.33754e+03 -2.69130e+04 -2.67216e+04 3.02292e+02 -1.37233e+02 + Pressure (bar) dVcoul/dl dVvdw/dl + -3.65296e+02 8.97698e+00 6.05381e+00 + + Total Virial (kJ/mol) + 2.43120e+03 -6.98762e+00 9.51030e+01 + -6.75444e+00 2.24529e+03 -4.06969e+01 + 9.50336e+01 -4.08237e+01 2.51504e+03 + + Pressure (bar) + -4.25137e+02 1.37702e+01 -1.55007e+02 + 1.34710e+01 -1.68439e+02 3.60660e+01 + -1.54918e+02 3.62287e+01 -5.02311e+02 + + + M E G A - F L O P S A C C O U N T I N G + + NB=Group-cutoff nonbonded kernels NxN=N-by-N cluster Verlet kernels + RF=Reaction-Field VdW=Van der Waals QSTab=quadratic-spline table + W3=SPC/TIP3p W4=TIP4p (single or pairs) + V&F=Potential and force V=Potential only F=Force only + + Computing: M-Number M-Flops % Flops +----------------------------------------------------------------------------- + NB Free energy kernel 469.832028 469.832 0.5 + Pair Search distance check 15.104240 135.938 0.2 + NxN Ewald Elec. + LJ [V&F] 357.853216 45447.358 50.8 + NxN Ewald Elec. [V&F] 353.795616 29718.832 33.2 + 1,4 nonbonded interactions 0.000501 0.045 0.0 + Calc Weights 3.789063 136.406 0.2 + Spread Q Bspline 161.666688 323.333 0.4 + Gather F Bspline 161.666688 970.000 1.1 + 3D-FFT 1450.142496 11601.140 13.0 + Solve PME 2.597184 166.220 0.2 + Reset In Box 0.032773 0.098 0.0 + CG-CoM 0.035294 0.106 0.0 + Bonds 0.001503 0.089 0.0 + Angles 0.001002 0.168 0.0 + Propers 0.000501 0.115 0.0 + Virial 1.443381 25.981 0.0 + Stop-CM 1.265542 12.655 0.0 + Calc-Ekin 1.263021 34.102 0.0 + Constraint-V 1.261017 11.349 0.0 + Constraint-Vir 1.261017 30.264 0.0 + Settle 0.842356 311.672 0.3 +----------------------------------------------------------------------------- + Total 89395.704 100.0 +----------------------------------------------------------------------------- + + + D O M A I N D E C O M P O S I T I O N S T A T I S T I C S + + av. #atoms communicated per step for force: 2 x 9223.4 + + +Dynamic load balancing report: + DLB was turned on during the run due to measured imbalance. + Average load imbalance: 12.6%. + The balanceable part of the MD step is 16%, load imbalance is computed from this. + Part of the total run time spent waiting due to load imbalance: 2.0%. + Steps where the load balancing was limited by -rdd, -rcon and/or -dds: X 0 % Y 0 % + + + R E A L C Y C L E A N D T I M E A C C O U N T I N G + +On 8 MPI ranks + + Computing: Num Num Call Wall time Giga-Cycles + Ranks Threads Count (s) total sum % +----------------------------------------------------------------------------- + Domain decomp. 8 1 13 0.011 0.120 0.3 + DD comm. load 8 1 12 0.000 0.000 0.0 + DD comm. bounds 8 1 11 0.000 0.001 0.0 + Neighbor search 8 1 13 0.016 0.176 0.4 + Comm. coord. 8 1 488 0.101 1.126 2.7 + Force 8 1 501 0.354 3.934 9.5 + Wait + Comm. F 8 1 501 0.079 0.876 2.1 + PME mesh 8 1 501 2.976 33.084 80.0 + NB X/F buffer ops. 8 1 1477 0.008 0.089 0.2 + Write traj. 8 1 2 0.001 0.016 0.0 + Update 8 1 2004 0.008 0.087 0.2 + Constraints 8 1 1004 0.005 0.058 0.1 + Comm. energies 8 1 1002 0.150 1.672 4.0 + Rest 0.011 0.121 0.3 +----------------------------------------------------------------------------- + Total 3.721 41.362 100.0 +----------------------------------------------------------------------------- + Breakdown of PME mesh computation +----------------------------------------------------------------------------- + PME redist. X/F 8 1 1503 1.057 11.752 28.4 + PME spread 8 1 1002 0.512 5.688 13.8 + PME gather 8 1 1002 0.287 3.188 7.7 + PME 3D-FFT 8 1 2004 0.153 1.703 4.1 + PME 3D-FFT Comm. 8 1 4008 0.890 9.889 23.9 + PME solve Elec 8 1 1002 0.077 0.853 2.1 +----------------------------------------------------------------------------- + + Core t (s) Wall t (s) (%) + Time: 29.713 3.721 798.5 + (ns/day) (hour/ns) +Performance: 23.267 1.032 +Finished mdrun on rank 0 Sun Aug 7 21:24:35 2022 + diff --git a/ensemble_md/tests/data/expanded.mdp b/ensemble_md/tests/data/expanded.mdp index 93a5e11b..ea65bc59 100755 --- a/ensemble_md/tests/data/expanded.mdp +++ b/ensemble_md/tests/data/expanded.mdp @@ -52,7 +52,7 @@ tcoupl = v-rescale nsttcouple = 1 tc_grps = System tau_t = 0.5 -ref_t = 298 +ref-t = 298 ; Pressure coupling is on for NPT pcoupl = no diff --git a/ensemble_md/tests/data/other_mdps/1.mdp b/ensemble_md/tests/data/other_mdps/1.mdp new file mode 100755 index 00000000..f92dadd1 --- /dev/null +++ b/ensemble_md/tests/data/other_mdps/1.mdp @@ -0,0 +1,108 @@ +; Run control +integrator = md-vv +tinit = 0 +dt = 0.002 +nsteps = 500 ; 10 ps +comm-mode = Linear +nstcomm = 1 +nstfout = 0 + +; Output control +nstlog = 100 +nstcalcenergy = 1 +nstenergy = 1000 +nstxout-compressed = 1000 + +; Neighborsearching and short-range nonbonded interactions +nstlist = 10 +ns_type = grid +pbc = xyz +rlist = 1.0 + +; Electrostatics +cutoff-scheme = verlet +coulombtype = PME +coulomb-modifier = Potential-shift-Verlet +rcoulomb-switch = 0.89 +rcoulomb = 0.9 + +; van der Waals +vdw-type = Cut-off +vdw-modifier = Potential-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.10 + +; EWALD/PME/PPPM parameters +fourier_nx = 0 +fourier_ny = 0 +fourier_nz = 0 +pme_order = 4 +ewald_rtol = 1e-05 +ewald_geometry = 3d +epsilon_surface = 0 + +; Temperature coupling +tcoupl = v-rescale +nsttcouple = 1 +tc_grps = System +tau_t = 0.5 +ref-t = 298 + +; Pressure coupling is on for NPT +pcoupl = no + +; refcoord_scaling should do nothing since there are no position restraints. + +gen_vel = yes +gen-temp = 298 +gen-seed = 6722267; need to randomize the seed each time. + +; options for bonds +constraints = h-bonds ; we only have C-H bonds here + +; Type of constraint algorithm +constraint-algorithm = lincs +continuation = no + +; Highest order in the expansion of the constraint coupling matrix +lincs-order = 12 +lincs-iter = 2 + +; Free energy calculation +free_energy = expanded +calc-lambda-neighbors = -1 +sc-alpha = 0.5 +couple-moltype = LIG +couple-lambda0 = vdw-q +couple-lambda1 = none +couple-intramol = no +init-lambda-state = 0 + +nstdhdl = 10 +dhdl-print-energy = total + +; Seed for Monte Carlo in lambda space +lmc-seed = 1000 +lmc-gibbsdelta = -1 +lmc-forced-nstart = 0 +symmetrized-transition-matrix = yes +nst-transition-matrix = 100000 +wl-scale = 0.8 +wl-ratio = 0.8 +init-wl-delta = 0.5 + +; expanded ensemble variables +nstexpanded = 10 +lmc-stats = wang-landau +lmc-move = metropolized-gibbs +lmc-weights-equil = wl-delta +weight-equil-wl-delta = 0.001 + +; lambda-states = 1 2 3 4 5 6 7 8 9 +vdw-lambdas = 0.1 0.20 0.30 0.40 0.50 0.60 0.70 0.8 1.00 diff --git a/ensemble_md/tests/data/other_mdps/2.mdp b/ensemble_md/tests/data/other_mdps/2.mdp new file mode 100755 index 00000000..945d65df --- /dev/null +++ b/ensemble_md/tests/data/other_mdps/2.mdp @@ -0,0 +1,111 @@ +; Run control +integrator = md-vv +tinit = 0 +dt = 0.002 +nsteps = 500 ; 10 ps +comm-mode = Linear +nstcomm = 1 +nstfout = 0 + +; Output control +nstlog = 100 +nstcalcenergy = 1 +nstenergy = 1000 +nstxout-compressed = 1000 + +; Neighborsearching and short-range nonbonded interactions +nstlist = 10 +ns_type = grid +pbc = xyz +rlist = 1.0 + +; Electrostatics +cutoff-scheme = verlet +coulombtype = PME +coulomb-modifier = Potential-shift-Verlet +rcoulomb-switch = 0.89 +rcoulomb = 0.9 + +; van der Waals +vdw-type = Cut-off +vdw-modifier = Potential-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.10 + +; EWALD/PME/PPPM parameters +fourier_nx = 0 +fourier_ny = 0 +fourier_nz = 0 +pme_order = 4 +ewald_rtol = 1e-05 +ewald_geometry = 3d +epsilon_surface = 0 + +; Temperature coupling +tcoupl = v-rescale +nsttcouple = 1 +tc_grps = System +tau_t = 0.5 +ref-t = 298 + +; Pressure coupling is on for NPT +pcoupl = no + +; refcoord_scaling should do nothing since there are no position restraints. + +gen_vel = yes +gen-temp = 298 +gen-seed = 6722267; need to randomize the seed each time. + +; options for bonds +constraints = h-bonds ; we only have C-H bonds here + +; Type of constraint algorithm +constraint-algorithm = lincs +continuation = no + +; Highest order in the expansion of the constraint coupling matrix +lincs-order = 12 +lincs-iter = 2 + +; Free energy calculation +free_energy = expanded +calc-lambda-neighbors = -1 +sc-alpha = 0.5 +couple-moltype = LIG +couple-lambda0 = vdw-q +couple-lambda1 = none +couple-intramol = no +init-lambda-state = 0 + +nstdhdl = 10 +dhdl-print-energy = total + +; Seed for Monte Carlo in lambda space +lmc-seed = 1000 +lmc-gibbsdelta = -1 +lmc-forced-nstart = 0 +symmetrized-transition-matrix = yes +nst-transition-matrix = 100000 +wl-scale = 0.8 +wl-ratio = 0.8 +init-wl-delta = 0.5 + +; expanded ensemble variables +nstexpanded = 10 +lmc-stats = wang-landau +lmc-move = metropolized-gibbs +lmc-weights-equil = wl-delta +weight-equil-wl-delta = 0.001 + +; lambda-states = 1 2 3 4 5 6 7 8 9 +coul-lambdas = 0.00 0.25 0.50 0.75 1.00 1.00 1.00 1.00 1.00 +vdw-lambdas = 0.00 0.00 0.00 0.00 0.00 0.25 0.50 0.75 1.00 +restraint-lambdas = 0.00 0.00 0.00 0.00 0.00 0.20 0.60 0.80 1.00 +init-lambda-weights = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 diff --git a/ensemble_md/tests/data/other_yamls/1.yaml b/ensemble_md/tests/data/other_yamls/1.yaml new file mode 100644 index 00000000..b23b9530 --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/1.yaml @@ -0,0 +1,6 @@ +# User-defined parameters +n_sim: 4 # Number of replica simulations +n_iterations: 10 # Number of iterations +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 'ensemble_md/tests/data/expanded.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/2.yaml b/ensemble_md/tests/data/other_yamls/2.yaml new file mode 100644 index 00000000..127f5c2b --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/2.yaml @@ -0,0 +1,9 @@ +# User-defined parameters +parallel: False # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4 # Number of replica simulations +n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mc_scheme: 'cool_scheme' +mdp: 'ensemble_md/tests/data/expanded.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/3.yaml b/ensemble_md/tests/data/other_yamls/3.yaml new file mode 100644 index 00000000..b7b56e53 --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/3.yaml @@ -0,0 +1,8 @@ +# User-defined parameters +parallel: False # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4.1 # Number of replica simulations +n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 'ensemble_md/tests/data/expanded.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/4.yaml b/ensemble_md/tests/data/other_yamls/4.yaml new file mode 100644 index 00000000..9b4984ad --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/4.yaml @@ -0,0 +1,8 @@ +# User-defined parameters +parallel: False # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4 # Number of replica simulations +n_iterations: 0 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 'ensemble_md/tests/data/expanded.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/5.yaml b/ensemble_md/tests/data/other_yamls/5.yaml new file mode 100644 index 00000000..412995ad --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/5.yaml @@ -0,0 +1,9 @@ +# User-defined parameters +parallel: False # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4 # Number of replica simulations +n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +n_ex: -1 +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 'ensemble_md/tests/data/expanded.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/6.yaml b/ensemble_md/tests/data/other_yamls/6.yaml new file mode 100644 index 00000000..fd3d8b09 --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/6.yaml @@ -0,0 +1,8 @@ +# User-defined parameters +parallel: False # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4 # Number of replica simulations +n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 3 # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/7.yaml b/ensemble_md/tests/data/other_yamls/7.yaml new file mode 100644 index 00000000..21953dfe --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/7.yaml @@ -0,0 +1,8 @@ +# User-defined parameters +parallel: 3 # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4 # Number of replica simulations +n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 'ensemble_md/tests/data/expanded.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/8.yaml b/ensemble_md/tests/data/other_yamls/8.yaml new file mode 100644 index 00000000..edfd5c3b --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/8.yaml @@ -0,0 +1,8 @@ +# User-defined parameters +parallel: False # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4 # Number of replica simulations +n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 'ensemble_md/tests/data/other_mdps/1.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/other_yamls/9.yaml b/ensemble_md/tests/data/other_yamls/9.yaml new file mode 100644 index 00000000..d0901290 --- /dev/null +++ b/ensemble_md/tests/data/other_yamls/9.yaml @@ -0,0 +1,8 @@ +# User-defined parameters +parallel: False # Whether the replicas of EXEE should be run in parallel or not +n_sim: 4 # Number of replica simulations +n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica +s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) +nst_sim: 500 # Number of simulation steps for each homogeneous replica +mdp: 'ensemble_md/tests/data/other_mdps/2.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/data/params.yaml b/ensemble_md/tests/data/params.yaml index 9ab63752..ff3ae7d3 100644 --- a/ensemble_md/tests/data/params.yaml +++ b/ensemble_md/tests/data/params.yaml @@ -2,6 +2,7 @@ parallel: False # Whether the replicas of EXEE should be run in parallel or not n_sim: 4 # Number of replica simulations n_iterations: 10 # Number of iterations +nst_sim: 500 # Number of simulation steps for each homogeneous replica s: 1 # Shift in λ ranges (e.g. s = 2 if λ_2 = (2, 3, 4) & λ_3 = (4, 5, 6)) nst_sim: 500 # Number of simulation steps for each homogeneous replica mdp: 'ensemble_md/tests/data/expanded.mdp' # The MDP template that has the whole range of λ values diff --git a/ensemble_md/tests/test_analyze_EEXE.py b/ensemble_md/tests/test_analyze_EEXE.py new file mode 100644 index 00000000..db3db9c3 --- /dev/null +++ b/ensemble_md/tests/test_analyze_EEXE.py @@ -0,0 +1,70 @@ +#################################################################### +# # +# 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 analyze_EEXE.py. +""" +import os +import pytest +import numpy as np +from ensemble_md import analyze_EEXE +from ensemble_md.exceptions import ParseError + +current_path = os.path.dirname(os.path.abspath(__file__)) +input_path = os.path.join(current_path, "data") + + +def test_parse_transmtx(): + A1, B1, C1 = analyze_EEXE.parse_transmtx(os.path.join(input_path, 'EXE.log')) + A1_expected = np.array([[0.5035941, 0.27816516, 0.15014392, 0.06404528, 0.01056806, 0.00325544], + [0.27816516, 0.1430219, 0.22873016, 0.12587893, 0.0782651, 0.07551135], # noqa: E128 + [0.15014392, 0.22873016, 0.12687601, 0.22781774, 0.19645289, 0.14166102], + [0.06404528, 0.12587893, 0.22781774, 0.0986531, 0.27133784, 0.18215007], + [0.01056806, 0.0782651, 0.19645289, 0.27133784, 0.06328789, 0.36931726], + [0.00325544, 0.07551135, 0.14166102, 0.18215007, 0.36931726, 0.2982755]]) + + B1_expected = np.array([[0.5, 0.34782609, 0.15000001, 0, 0, 0], + [0.34782609, 0.18181819, 0.15789473, 0.17647059, 0.10526316, 0.125], # noqa: E128 + [0.15000001, 0.15789473, 0, 0.14285715, 0.4375, 0.07692308], + [0, 0.17647059, 0.14285715, 0.16666667, 0.2857143, 0.18181819], + [0, 0.10526316, 0.4375, 0.2857143, 0, 0.23076923], + [0, 0.125, 0.07692308, 0.18181819, 0.23076923, 0.2]]) + + C1_expected = np.array([[3.5941000e-03, -6.9660930e-02, 1.4391000e-04, 6.4045280e-02, 1.0568060e-02, 3.2554400e-03], # noqa: E501 + [-6.9660930e-02, -3.8796290e-02, 7.0835430e-02, -5.0591660e-02, -2.6998060e-02, -4.9488650e-02], # noqa: E128 + [1.4391000e-04, 7.0835430e-02, 1.2687601e-01, 8.4960590e-02, -2.4104711e-01, 6.4737940e-02], + [6.4045280e-02, -5.0591660e-02, 8.4960590e-02, -6.8013570e-02, -1.4376460e-02, 3.3188000e-04], + [1.0568060e-02, -2.6998060e-02, -2.4104711e-01, -1.4376460e-02, 6.3287890e-02, 1.3854803e-01], + [3.2554400e-03, -4.9488650e-02, 6.4737940e-02, 3.3188000e-04, 1.3854803e-01, 9.8275500e-02]]) + + np.testing.assert_array_almost_equal(A1, A1_expected) + np.testing.assert_array_almost_equal(B1, B1_expected) + np.testing.assert_array_almost_equal(C1, C1_expected) + + log = os.path.join(input_path, 'EXE_0.log') + with pytest.raises(ParseError, match=f'No transition matrices found in {log}.'): + A2, B2, C2 = analyze_EEXE.parse_transmtx(log) + + +def test_plot_matrix(): + """ + We can only check if the figures are generated. Not really able to check how they look like. + """ + A1, B1, C1 = analyze_EEXE.parse_transmtx(os.path.join(input_path, 'EXE.log')) + analyze_EEXE.plot_matrix(A1, 'test_1.png') + analyze_EEXE.plot_matrix(B1, 'test_2.png') + analyze_EEXE.plot_matrix(C1, 'test_3.png') + + assert os.path.exists('test_1.png') is True + assert os.path.exists('test_2.png') is True + assert os.path.exists('test_3.png') is True + + os.remove('test_1.png') + os.remove('test_2.png') + os.remove('test_3.png') diff --git a/ensemble_md/tests/test_ensemble_EXE.py b/ensemble_md/tests/test_ensemble_EXE.py index eeb6d3c4..7dc94da1 100644 --- a/ensemble_md/tests/test_ensemble_EXE.py +++ b/ensemble_md/tests/test_ensemble_EXE.py @@ -8,14 +8,16 @@ # # #################################################################### """ -Unit tests for the module ensemble_EXE.py. The working directory shoul be the home directory of the package. +Unit tests for the module ensemble_EXE.py. """ import os +import pytest import random import numpy as np import ensemble_md import ensemble_md.gmx_parser as gmx_parser from ensemble_md.ensemble_EXE import EnsembleEXE +from ensemble_md.exceptions import ParameterError current_path = os.path.dirname(os.path.abspath(__file__)) input_path = os.path.join(current_path, "data") @@ -23,7 +25,7 @@ class Test_EnsembleEXE: - def test_init(self): + def test_init_1(self): k = 1.380649e-23 NA = 6.0221408e23 assert EEXE.mc_scheme == "metropolis" @@ -45,9 +47,6 @@ def test_init(self): {3, 4, 5, 6, 7, 8}, ] assert EEXE.nst_sim == 500 - - def test_map_lambda2state(self): - EEXE.map_lambda2state() assert EEXE.lambda_dict == { (0, 0): 0, (0.25, 0): 1, @@ -66,6 +65,75 @@ def test_map_lambda2state(self): [(0.75, 0.0), (1.0, 0.0), (1.0, 0.25), (1.0, 0.5), (1.0, 0.75), (1.0, 1.0)], ] + def test_init_2(self): + yaml_1 = os.path.join(input_path, "other_yamls/1.yaml") + with pytest.raises(ParameterError, match=f"Required parameter 'parallel' not specified in {yaml_1}"): + E1 = EnsembleEXE(yaml_1) # noqa: F841 + + yaml_2 = os.path.join(input_path, "other_yamls/2.yaml") + with pytest.raises(ParameterError, match="The specified MC scheme is not available. Options include 'same-state', 'metropolis', and 'metropolis-eq'."): # noqa: E501 + E2 = EnsembleEXE(yaml_2) # noqa: F841 + + yaml_3 = os.path.join(input_path, "other_yamls/3.yaml") + with pytest.raises(ParameterError, match="The parameter 'n_sim' should be an integer."): + E3 = EnsembleEXE(yaml_3) # noqa: F841 + + yaml_4 = os.path.join(input_path, "other_yamls/4.yaml") + with pytest.raises(ParameterError, match="The parameter 'n_iterations' should be positive."): + E4 = EnsembleEXE(yaml_4) # noqa: F841 + + yaml_5 = os.path.join(input_path, "other_yamls/5.yaml") + with pytest.raises(ParameterError, match="The parameter 'n_ex' should be non-negative."): + E5 = EnsembleEXE(yaml_5) # noqa: F841 + + yaml_6 = os.path.join(input_path, "other_yamls/6.yaml") + with pytest.raises(ParameterError, match="The parameter 'mdp' should be a string."): + E6 = EnsembleEXE(yaml_6) # noqa: F841 + + yaml_7 = os.path.join(input_path, "other_yamls/7.yaml") + with pytest.raises(ParameterError, match="The parameter 'parallel' should be a boolean variable."): + E7 = EnsembleEXE(yaml_7) # noqa: F841 + + yaml_8 = os.path.join(input_path, "other_yamls/8.yaml") + E8 = EnsembleEXE(yaml_8) + assert E8.lambda_dict == { + (0.1,): 0, + (0.2,): 1, + (0.3,): 2, + (0.4,): 3, + (0.5,): 4, + (0.6,): 5, + (0.7,): 6, + (0.8,): 7, + (1,): 8, + } # noqa: E501 + assert E8.lambda_ranges == [ + [(0.1,), (0.2,), (0.3,), (0.4,), (0.5,), (0.6,)], + [(0.2,), (0.3,), (0.4,), (0.5,), (0.6,), (0.7,)], + [(0.3,), (0.4,), (0.5,), (0.6,), (0.7,), (0.8,)], + [(0.4,), (0.5,), (0.6,), (0.7,), (0.8,), (1.0,)], + ] + + yaml_9 = os.path.join(input_path, "other_yamls/9.yaml") + E9 = EnsembleEXE(yaml_9) + assert E9.lambda_dict == { + (0, 0, 0): 0, + (0.25, 0, 0): 1, + (0.5, 0, 0): 2, + (0.75, 0, 0): 3, + (1, 0, 0): 4, + (1, 0.25, 0.2): 5, + (1, 0.5, 0.6): 6, + (1, 0.75, 0.8): 7, + (1, 1, 1): 8, + } # noqa: E501 + assert E9.lambda_ranges == [ + [(0.0, 0.0, 0.0), (0.25, 0.0, 0.0), (0.5, 0.0, 0.0), (0.75, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.25, 0.2)], + [(0.25, 0.0, 0.0), (0.5, 0.0, 0.0), (0.75, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.25, 0.2), (1.0, 0.5, 0.6)], + [(0.5, 0.0, 0.0), (0.75, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.25, 0.2), (1.0, 0.5, 0.6), (1.0, 0.75, 0.8)], + [(0.75, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.25, 0.2), (1.0, 0.5, 0.6), (1.0, 0.75, 0.8), (1.0, 1.0, 1.0)], + ] + def test_print_params(self, capfd): # capfd is a fixture in pytest for testing STDOUT EEXE.print_params() @@ -182,24 +250,28 @@ def test_extract_final_log_info(self): assert equil_bools == [False, False, False, False] def test_propose_swaps(self): - pass - """ - # This needs to be rewritten for multiple swapping. random.seed(0) - EEXE.n_sim = 8 - EEXE.n_pairs = 5 - EEXE.state_ranges = [ - set(range(i, i + 5)) for i in range(EEXE.n_sim) - ] # 12 states, 5 for each replica - states = [5, 2, 2, 8] + EEXE.n_sim = 4 + EEXE.state_ranges = [set(range(i, i + 5)) for i in range(EEXE.n_sim)] # 5 states per replica + states = [5, 2, 2, 7] # This would lead to the swappables: [(0, 1), (0, 2), (1, 2)] + + # Case 1: Neighboring swapping (n_ex = 0 --> swappables = [(0, 1), (1, 2)]) + EEXE.n_ex = 0 swap_list = EEXE.propose_swaps(states) - assert EEXE.n_pairs == 4 - assert swap_list == [ - (3, 4), - (5, 6), - (0, 1), - ] # The remaining pair of (2, 7) is not swappable - """ + assert swap_list == [(1, 2)] + + # Case 2: Multiple swaps (n_ex = 3) + EEXE.n_ex = 5 + swap_list = EEXE.propose_swaps(states) + assert swap_list == [(1, 2), (0, 2), (0, 1), (0, 2), (0, 2)] + + # Case 3: Empty swappable list + states = [10, 10, 10, 10] + swap_list = EEXE.propose_swaps(states) + assert swap_list is None + + def test_gest_swapped_configus(self): + pass def test_calc_prob_acc(self): pass diff --git a/ensemble_md/tests/test_gmx_parser.py b/ensemble_md/tests/test_gmx_parser.py index e69de29b..3bddc4a7 100644 --- a/ensemble_md/tests/test_gmx_parser.py +++ b/ensemble_md/tests/test_gmx_parser.py @@ -0,0 +1,12 @@ +#################################################################### +# # +# 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 gmx_parser.py. +"""