From 837da49a527da5b6b6c43108383ca840ae22127d Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 22 Sep 2024 14:03:25 +0200 Subject: [PATCH] small improvement of the MC swap --- docs/source/chapters/chapter8.rst | 315 ++++++++++++------------------ docs/source/chapters/chapter9.rst | 306 +++++++++++++++++++---------- docs/source/index.rst | 8 +- 3 files changed, 328 insertions(+), 301 deletions(-) diff --git a/docs/source/chapters/chapter8.rst b/docs/source/chapters/chapter8.rst index 5655afa..bd77231 100644 --- a/docs/source/chapters/chapter8.rst +++ b/docs/source/chapters/chapter8.rst @@ -1,141 +1,114 @@ -.. _chapter8-label: +.. _chapter9-label: -Monte Carlo insert -================== +Monte Carlo swap +================ -Here, a *Monte Carlo* simulation is implemented in the Grand Canonical ensemble, -where the number of atoms in the system fluctuates. The principle of the -simulation resembles the Monte Carlo move from :ref:`chapter6-label`: +A Monte Carlo swap move consists of randomly attempting to exchange the positions of two +particles. Just like in the case of Monte Carlo displace move, the swap is accepted +based on energy criteria. For systems where the dynamics are slow, such as in glasses, a Monte +Carlo swap move can considerably speed up equilibration. -- 1) We start from a given initial configuration, and measure the potential +The principle of such Monte Carlo swap simulation is the following: + +- 1) We start from a given intial configuration, and measure the potential energy, :math:`E_\text{pot}^\text{initial}`. -- 2) A random number is chosen, and depending on its value, either a new particle - will try to be inserted, or an existing particle will try to be deleted. -- 3) The energy of the system after the insertion/deletion, - :math:`E_\text{pot}^\text{trial}`, is measured. -- 4) We then decide to keep or reject the move by calculating - the difference in energy between the trial and the initial configurations: +- 2) Two particles are chosen randomly, and their respecitive positions are exchanced. +- 3) The energy of the system after the move, :math:`E_\text{pot}^\text{trial}`, is measured. +- 4) Just like with Monte Carlo move, we then have to decide to keep or reject + the move. This is done by calculating the difference in energy between the + trial and the initial configurations: :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`. + + - If :math:`\Delta E < 0`, then the move is automatically accepted. + - If :math:`\Delta E > 0`, then the move is accepted with a probability given + by the Boltzmann factor :math:`\exp{- \beta \Delta E}`, where + :math:`\beta = 1 / k_\text{B} T` and :math:`T` is the imposed temperature. + - 5) Steps 1-4 are repeated a large number of times, generating a broad range of possible configurations. +In general, both Monte Carlo swap moves and displacements are being performed +within the same simulation. + Implementation -------------- -Let us add the following method called *monte_carlo_exchange* to the *MonteCarlo* class: +Let us add a method named *monte_carlo_swap* to the *MonteCarlo* class: .. label:: start_MonteCarlo_class .. code-block:: python - def monte_carlo_exchange(self): - if self.desired_mu is not None: - # The first step is to make a copy of the previous state - # Since atoms numbers are evolving, its also important to store the - # neighbor, sigma, and epsilon lists - self.Epot = self.compute_potential() # TOFIX: not necessary every time + def monte_carlo_swap(self): + if self.swap_type[0] is not None: + self.update_neighbor_lists() + self.update_cross_coefficients() + if hasattr(self, 'Epot') is False: + self.Epot = self.compute_potential() + initial_Epot = self.Epot initial_positions = copy.deepcopy(self.atoms_positions) - initial_number_atoms = copy.deepcopy(self.number_atoms) - initial_neighbor_lists = copy.deepcopy(self.neighbor_lists) - initial_sigma_lists = copy.deepcopy(self.sigma_ij_list) - initial_epsilon_lists = copy.deepcopy(self.epsilon_ij_list) - # Apply a 50-50 probability to insert or delete - insert_or_delete = np.random.random() - if np.random.random() < insert_or_delete: - self.monte_carlo_insert() - else: - self.monte_carlo_delete() - if np.random.random() < self.acceptation_probability: # accepted move - # Update the success counters - if np.random.random() < insert_or_delete: - self.successful_insert += 1 - else: - self.successful_delete += 1 - else: - # Reject the new position, revert to inital position + # Pick an atom of type one randomly + atom_id_1 = np.random.randint(self.number_atoms[self.swap_type[0]]) + # Pick an atom of type two randomly + atom_id_2 = np.random.randint(self.number_atoms[self.swap_type[1]]) + + shift_1 = 0 + for N in self.number_atoms[:self.swap_type[0]]: + shift_1 += N + shift_2 = 0 + for N in self.number_atoms[:self.swap_type[1]]: + shift_2 += N + # attempt to swap the position of the atoms + position1 = copy.deepcopy(self.atoms_positions[shift_1+atom_id_1]) + position2 = copy.deepcopy(self.atoms_positions[shift_2+atom_id_2]) + self.atoms_positions[shift_2+atom_id_2] = position1 + self.atoms_positions[shift_1+atom_id_1] = position2 + # force the recalculation of neighbor list + initial_atoms_sigma = self.atoms_sigma + initial_atoms_epsilon = self.atoms_epsilon + initial_atoms_mass = self.atoms_mass + initial_atoms_type = self.atoms_type + initial_sigma_ij_list = self.sigma_ij_list + initial_epsilon_ij_list = self.epsilon_ij_list + initial_neighbor_lists = self.neighbor_lists + self.update_neighbor_lists(force_update=True) + self.identify_atom_properties() + self.update_cross_coefficients(force_update=True) + # Measure the potential energy of the new configuration + trial_Epot = self.compute_potential() + # Evaluate whether the new configuration should be kept or not + beta = 1/self.desired_temperature + delta_E = trial_Epot-initial_Epot + random_number = np.random.random() # random number between 0 and 1 + acceptation_probability = np.min([1, np.exp(-beta*delta_E)]) + if random_number <= acceptation_probability: # Accept new position + self.Epot = trial_Epot + self.successful_swap += 1 + else: # Reject new position + self.atoms_positions = initial_positions # Revert to initial positions + self.failed_swap += 1 + self.atoms_sigma = initial_atoms_sigma + self.atoms_epsilon = initial_atoms_epsilon + self.atoms_mass = initial_atoms_mass + self.atoms_type = initial_atoms_type + self.sigma_ij_list = initial_sigma_ij_list + self.epsilon_ij_list = initial_epsilon_ij_list self.neighbor_lists = initial_neighbor_lists - self.sigma_ij_list = initial_sigma_lists - self.epsilon_ij_list = initial_epsilon_lists - self.atoms_positions = initial_positions - self.number_atoms = initial_number_atoms - # Update the failed counters - if np.random.random() < insert_or_delete: - self.failed_insert += 1 - else: - self.failed_delete += 1 - + .. label:: end_MonteCarlo_class -First, the potential energy is calculated, and the current state of the -simulation, such as positions and atom numbers, is saved. Then, an insertion -trial or a deletion trial is made, each with a probability of 0.5. The -*monte_carlo_insert* and *monte_carlo_delete* methods, which are implemented -below, both calculate the *acceptance_probability*. A random number is selected, -and if that number is larger than the acceptance probability, the system reverts -to its previous position. If the random number is lower than the acceptance -probability, nothing happens, meaning that the last trial is accepted. - -Then, let us write the *monte_carlo_insert()* method: +Let us initialise swap counter: .. label:: start_MonteCarlo_class .. code-block:: python - def monte_carlo_insert(self): - self.number_atoms[self.inserted_type] += 1 - new_atom_position = np.random.random(3)*np.diff(self.box_boundaries).T \ - - np.diff(self.box_boundaries).T/2 - shift_id = 0 - for N in self.number_atoms[:self.inserted_type]: - shift_id += N - self.atoms_positions = np.vstack([self.atoms_positions[:shift_id], - new_atom_position, - self.atoms_positions[shift_id:]]) - self.update_neighbor_lists() - self.identify_atom_properties() - self.update_cross_coefficients() - trial_Epot = self.compute_potential() - Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type]) - beta = 1/self.desired_temperature - Nat = np.sum(self.number_atoms) # Number atoms, should it really be N? of N (type) ? - Vol = np.prod(self.box_size[:3]) # box volume - # dimension of 3 is enforced in the power of the Lambda - self.acceptation_probability = np.min([1, Vol/(Lambda**3*Nat) \ - *np.exp(beta*(self.desired_mu-trial_Epot+self.Epot))]) - -.. label:: end_MonteCarlo_class - -After trying to insert a new particle, neighbor lists and cross coefficients -must be re-evaluated. Then, the acceptance probability is calculated. - -Let us add the very similar *monte_carlo_delete()* method: - -.. label:: start_MonteCarlo_class - -.. code-block:: python - - def monte_carlo_delete(self): - # Pick one atom to delete randomly - atom_id = np.random.randint(self.number_atoms[self.inserted_type]) - self.number_atoms[self.inserted_type] -= 1 - if self.number_atoms[self.inserted_type] > 0: - shift_id = 0 - for N in self.number_atoms[:self.inserted_type]: - shift_id += N - self.atoms_positions = np.delete(self.atoms_positions, shift_id+atom_id, axis=0) - self.update_neighbor_lists() - self.identify_atom_properties() - self.update_cross_coefficients() - trial_Epot = self.compute_potential() - Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type]) - beta = 1/self.desired_temperature - Nat = np.sum(self.number_atoms) # Number atoms, should it really be N? of N (type) ? - Vol = np.prod(self.box_size[:3]) # box volume - # dimension of 3 is enforced in the power of the Lambda - self.acceptation_probability = np.min([1, (Lambda**3 *(Nat-1)/Vol) \ - *np.exp(-beta*(self.desired_mu+trial_Epot-self.Epot))]) - else: - print("Error: no more atoms to delete") + class MonteCarlo(Measurements): + def __init__(self, + (...) + self.failed_move = 0 + self.successful_swap = 0 + self.failed_swap = 0 .. label:: end_MonteCarlo_class @@ -149,8 +122,7 @@ Complete the *__init__* method as follows: def __init__(self, (...) displace_mc = None, - desired_mu = None, - inserted_type = 0, + swap_type = [None, None], .. label:: end_MonteCarlo_class @@ -164,28 +136,7 @@ and def __init__(self, (...) self.displace_mc = displace_mc - self.desired_mu = desired_mu - self.inserted_type = inserted_type - -.. label:: end_MonteCarlo_class - -Let us also normalize the "desired_mu": - -.. label:: start_MonteCarlo_class - -.. code-block:: python - - class MonteCarlo(Measurements): - def __init__(self, - (...) - self.nondimensionalize_units(["desired_temperature", "displace_mc"]) - self.nondimensionalize_units(["desired_mu"]) - self.successful_move = 0 - self.failed_move = 0 - self.successful_insert = 0 - self.failed_insert = 0 - self.successful_delete = 0 - self.failed_delete = 0 + self.swap_type = swap_type .. label:: end_MonteCarlo_class @@ -198,47 +149,16 @@ Finally, the *monte_carlo_exchange()* method must be included in the run: def run(self): (...) self.monte_carlo_move() - self.monte_carlo_exchange() - self.wrap_in_box() - (...) + self.monte_carlo_swap() .. label:: end_MonteCarlo_class -We need to calculate :math:`\Lambda`: - -.. label:: start_MonteCarlo_class - -.. code-block:: python - - def calculate_Lambda(self, mass): - """Estimate the de Broglie wavelength.""" - T = self.desired_temperature # N - return 1/np.sqrt(2*np.pi*mass*T) - -.. label:: end_MonteCarlo_class - -To output the density, let us add the following method to the *Measurements* class: - -.. label:: start_Measurements_class - -.. code-block:: python - - def calculate_density(self): - """Calculate the mass density.""" - # TOFIX: not used yet - volume = np.prod(self.box_size[:3]) # Unitless - total_mass = np.sum(self.atoms_mass) # Unitless - return total_mass/volume # Unitless - -.. label:: end_Measurements_class - Test the code ------------- -One can use a similar test as previously, but with an imposed chemical -potential *desired_mu*: +Let's test the Monte Carlo swap. -.. label:: start_test_8a_class +.. label:: start_test_9a_class .. code-block:: python @@ -248,56 +168,65 @@ potential *desired_mu*: import os # Define atom number of each group - nmb_1= 50 + nmb_1 = 50 + nmb_2 = 50 # New group for testing swaps # Define LJ parameters (sigma) - sig_1 = 3*ureg.angstrom + sig_1 = 3 * ureg.angstrom + sig_2 = 4 * ureg.angstrom # Different sigma for group 2 # Define LJ parameters (epsilon) - eps_1 = 0.1*ureg.kcal/ureg.mol + eps_1 = 0.1 * ureg.kcal / ureg.mol + eps_2 = 0.15 * ureg.kcal / ureg.mol # Different epsilon for group 2 # Define atom mass - mss_1 = 10*ureg.gram/ureg.mol + mss_1 = 10 * ureg.gram / ureg.mol + mss_2 = 12 * ureg.gram / ureg.mol # Different mass for group 2 # Define box size - L = 20*ureg.angstrom + L = 20 * ureg.angstrom # Define a cut off - rc = 2.5*sig_1 + rc = 2.5 * sig_1 # Pick the desired temperature - T = 300*ureg.kelvin - # choose the desired_mu - desired_mu = -3*ureg.kcal/ureg.mol + T = 300 * ureg.kelvin # Initialize the prepare object mc = MonteCarlo( - ureg = ureg, + ureg=ureg, maximum_steps=100, thermo_period=10, dumping_period=10, - number_atoms=[nmb_1], - epsilon=[eps_1], # kcal/mol - sigma=[sig_1], # A - atom_mass=[mss_1], # g/mol - box_dimensions=[L, L, L], # A + number_atoms=[nmb_1, nmb_2], # Include two groups of atoms for swap + epsilon=[eps_1, eps_2], # kcal/mol + sigma=[sig_1, sig_2], # A + atom_mass=[mss_1, mss_2], # g/mol + box_dimensions=[L, L, L], # A cut_off=rc, thermo_outputs="Epot-press", - desired_temperature=T, # K + desired_temperature=T, # K neighbor=1, - desired_mu = desired_mu, + swap_type=[0, 1] # Enable Monte Carlo swap between groups 1 and 2 ) + + # Run the Monte Carlo simulation mc.run() # Test function using pytest def test_output_files(): assert os.path.exists("Outputs/dump.mc.lammpstrj"), \ - "Test failed: dump file was not created" + "Test failed: dump file was not created" assert os.path.exists("Outputs/simulation.log"), \ - "Test failed: log file was not created" + "Test failed: log file was not created" print("Test passed") + # Test the swap counters + def test_swap_counters(): + assert mc.successful_swap + mc.failed_swap > 0, \ + "Test failed: No swaps were attempted" + print("Swap test passed") + # If the script is run directly, execute the tests if __name__ == "__main__": import pytest # Run pytest programmatically pytest.main(["-s", __file__]) -.. label:: end_test_8a_class +.. label:: end_test_9a_class + -The evolution of the potential energy as a function of the -number of steps is written in the *Outputs/Epot.dat*. diff --git a/docs/source/chapters/chapter9.rst b/docs/source/chapters/chapter9.rst index 89c01f6..5655afa 100644 --- a/docs/source/chapters/chapter9.rst +++ b/docs/source/chapters/chapter9.rst @@ -1,87 +1,141 @@ -.. _chapter9-label: +.. _chapter8-label: -Monte Carlo swap -================ +Monte Carlo insert +================== -A Monte Carlo swap move consists of randomly attempting to exchange the positions of two -particles. Just like in the case of Monte Carlo displace move, the swap is accepted -based on energy criteria. For systems where the dynamics are slow, such as in glasses, a Monte -Carlo swap move can considerably speed up equilibration. +Here, a *Monte Carlo* simulation is implemented in the Grand Canonical ensemble, +where the number of atoms in the system fluctuates. The principle of the +simulation resembles the Monte Carlo move from :ref:`chapter6-label`: + +- 1) We start from a given initial configuration, and measure the potential + energy, :math:`E_\text{pot}^\text{initial}`. +- 2) A random number is chosen, and depending on its value, either a new particle + will try to be inserted, or an existing particle will try to be deleted. +- 3) The energy of the system after the insertion/deletion, + :math:`E_\text{pot}^\text{trial}`, is measured. +- 4) We then decide to keep or reject the move by calculating + the difference in energy between the trial and the initial configurations: + :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`. +- 5) Steps 1-4 are repeated a large number of times, generating a broad range of + possible configurations. + +Implementation +-------------- + +Let us add the following method called *monte_carlo_exchange* to the *MonteCarlo* class: .. label:: start_MonteCarlo_class .. code-block:: python - def monte_carlo_swap(self): - if self.swap_type[0] is not None: - self.update_neighbor_lists() - self.update_cross_coefficients() - if hasattr(self, 'Epot') is False: - self.Epot = self.compute_potential() - initial_Epot = self.Epot + def monte_carlo_exchange(self): + if self.desired_mu is not None: + # The first step is to make a copy of the previous state + # Since atoms numbers are evolving, its also important to store the + # neighbor, sigma, and epsilon lists + self.Epot = self.compute_potential() # TOFIX: not necessary every time initial_positions = copy.deepcopy(self.atoms_positions) - # Pick an atom of type one randomly - atom_id_1 = np.random.randint(self.number_atoms[self.swap_type[0]]) - # Pick an atom of type two randomly - atom_id_2 = np.random.randint(self.number_atoms[self.swap_type[1]]) - - shift_1 = 0 - for N in self.number_atoms[:self.swap_type[0]]: - shift_1 += N - shift_2 = 0 - for N in self.number_atoms[:self.swap_type[1]]: - shift_2 += N - # attempt to swap the position of the atoms - position1 = copy.deepcopy(self.atoms_positions[shift_1+atom_id_1]) - position2 = copy.deepcopy(self.atoms_positions[shift_2+atom_id_2]) - self.atoms_positions[shift_2+atom_id_2] = position1 - self.atoms_positions[shift_1+atom_id_1] = position2 - # force the recalculation of neighbor list - initial_atoms_sigma = self.atoms_sigma - initial_atoms_epsilon = self.atoms_epsilon - initial_atoms_mass = self.atoms_mass - initial_atoms_type = self.atoms_type - initial_sigma_ij_list = self.sigma_ij_list - initial_epsilon_ij_list = self.epsilon_ij_list - initial_neighbor_lists = self.neighbor_lists - self.update_neighbor_lists(force_update=True) - self.identify_atom_properties() - self.update_cross_coefficients(force_update=True) - # Measure the potential energy of the new configuration - trial_Epot = self.compute_potential() - # Evaluate whether the new configuration should be kept or not - beta = 1/self.desired_temperature - delta_E = trial_Epot-initial_Epot - random_number = np.random.random() # random number between 0 and 1 - acceptation_probability = np.min([1, np.exp(-beta*delta_E)]) - if random_number <= acceptation_probability: # Accept new position - self.Epot = trial_Epot - self.successful_swap += 1 - else: # Reject new position - self.atoms_positions = initial_positions # Revert to initial positions - self.failed_swap += 1 - self.atoms_sigma = initial_atoms_sigma - self.atoms_epsilon = initial_atoms_epsilon - self.atoms_mass = initial_atoms_mass - self.atoms_type = initial_atoms_type - self.sigma_ij_list = initial_sigma_ij_list - self.epsilon_ij_list = initial_epsilon_ij_list + initial_number_atoms = copy.deepcopy(self.number_atoms) + initial_neighbor_lists = copy.deepcopy(self.neighbor_lists) + initial_sigma_lists = copy.deepcopy(self.sigma_ij_list) + initial_epsilon_lists = copy.deepcopy(self.epsilon_ij_list) + # Apply a 50-50 probability to insert or delete + insert_or_delete = np.random.random() + if np.random.random() < insert_or_delete: + self.monte_carlo_insert() + else: + self.monte_carlo_delete() + if np.random.random() < self.acceptation_probability: # accepted move + # Update the success counters + if np.random.random() < insert_or_delete: + self.successful_insert += 1 + else: + self.successful_delete += 1 + else: + # Reject the new position, revert to inital position self.neighbor_lists = initial_neighbor_lists - + self.sigma_ij_list = initial_sigma_lists + self.epsilon_ij_list = initial_epsilon_lists + self.atoms_positions = initial_positions + self.number_atoms = initial_number_atoms + # Update the failed counters + if np.random.random() < insert_or_delete: + self.failed_insert += 1 + else: + self.failed_delete += 1 + .. label:: end_MonteCarlo_class -Let us initialise swap counter: +First, the potential energy is calculated, and the current state of the +simulation, such as positions and atom numbers, is saved. Then, an insertion +trial or a deletion trial is made, each with a probability of 0.5. The +*monte_carlo_insert* and *monte_carlo_delete* methods, which are implemented +below, both calculate the *acceptance_probability*. A random number is selected, +and if that number is larger than the acceptance probability, the system reverts +to its previous position. If the random number is lower than the acceptance +probability, nothing happens, meaning that the last trial is accepted. + +Then, let us write the *monte_carlo_insert()* method: .. label:: start_MonteCarlo_class .. code-block:: python - class MonteCarlo(Measurements): - def __init__(self, - (...) - self.failed_move = 0 - self.successful_swap = 0 - self.failed_swap = 0 + def monte_carlo_insert(self): + self.number_atoms[self.inserted_type] += 1 + new_atom_position = np.random.random(3)*np.diff(self.box_boundaries).T \ + - np.diff(self.box_boundaries).T/2 + shift_id = 0 + for N in self.number_atoms[:self.inserted_type]: + shift_id += N + self.atoms_positions = np.vstack([self.atoms_positions[:shift_id], + new_atom_position, + self.atoms_positions[shift_id:]]) + self.update_neighbor_lists() + self.identify_atom_properties() + self.update_cross_coefficients() + trial_Epot = self.compute_potential() + Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type]) + beta = 1/self.desired_temperature + Nat = np.sum(self.number_atoms) # Number atoms, should it really be N? of N (type) ? + Vol = np.prod(self.box_size[:3]) # box volume + # dimension of 3 is enforced in the power of the Lambda + self.acceptation_probability = np.min([1, Vol/(Lambda**3*Nat) \ + *np.exp(beta*(self.desired_mu-trial_Epot+self.Epot))]) + +.. label:: end_MonteCarlo_class + +After trying to insert a new particle, neighbor lists and cross coefficients +must be re-evaluated. Then, the acceptance probability is calculated. + +Let us add the very similar *monte_carlo_delete()* method: + +.. label:: start_MonteCarlo_class + +.. code-block:: python + + def monte_carlo_delete(self): + # Pick one atom to delete randomly + atom_id = np.random.randint(self.number_atoms[self.inserted_type]) + self.number_atoms[self.inserted_type] -= 1 + if self.number_atoms[self.inserted_type] > 0: + shift_id = 0 + for N in self.number_atoms[:self.inserted_type]: + shift_id += N + self.atoms_positions = np.delete(self.atoms_positions, shift_id+atom_id, axis=0) + self.update_neighbor_lists() + self.identify_atom_properties() + self.update_cross_coefficients() + trial_Epot = self.compute_potential() + Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type]) + beta = 1/self.desired_temperature + Nat = np.sum(self.number_atoms) # Number atoms, should it really be N? of N (type) ? + Vol = np.prod(self.box_size[:3]) # box volume + # dimension of 3 is enforced in the power of the Lambda + self.acceptation_probability = np.min([1, (Lambda**3 *(Nat-1)/Vol) \ + *np.exp(-beta*(self.desired_mu+trial_Epot-self.Epot))]) + else: + print("Error: no more atoms to delete") .. label:: end_MonteCarlo_class @@ -95,7 +149,8 @@ Complete the *__init__* method as follows: def __init__(self, (...) displace_mc = None, - swap_type = [None, None], + desired_mu = None, + inserted_type = 0, .. label:: end_MonteCarlo_class @@ -109,7 +164,28 @@ and def __init__(self, (...) self.displace_mc = displace_mc - self.swap_type = swap_type + self.desired_mu = desired_mu + self.inserted_type = inserted_type + +.. label:: end_MonteCarlo_class + +Let us also normalize the "desired_mu": + +.. label:: start_MonteCarlo_class + +.. code-block:: python + + class MonteCarlo(Measurements): + def __init__(self, + (...) + self.nondimensionalize_units(["desired_temperature", "displace_mc"]) + self.nondimensionalize_units(["desired_mu"]) + self.successful_move = 0 + self.failed_move = 0 + self.successful_insert = 0 + self.failed_insert = 0 + self.successful_delete = 0 + self.failed_delete = 0 .. label:: end_MonteCarlo_class @@ -122,16 +198,47 @@ Finally, the *monte_carlo_exchange()* method must be included in the run: def run(self): (...) self.monte_carlo_move() - self.monte_carlo_swap() + self.monte_carlo_exchange() + self.wrap_in_box() + (...) .. label:: end_MonteCarlo_class +We need to calculate :math:`\Lambda`: + +.. label:: start_MonteCarlo_class + +.. code-block:: python + + def calculate_Lambda(self, mass): + """Estimate the de Broglie wavelength.""" + T = self.desired_temperature # N + return 1/np.sqrt(2*np.pi*mass*T) + +.. label:: end_MonteCarlo_class + +To output the density, let us add the following method to the *Measurements* class: + +.. label:: start_Measurements_class + +.. code-block:: python + + def calculate_density(self): + """Calculate the mass density.""" + # TOFIX: not used yet + volume = np.prod(self.box_size[:3]) # Unitless + total_mass = np.sum(self.atoms_mass) # Unitless + return total_mass/volume # Unitless + +.. label:: end_Measurements_class + Test the code ------------- -Let's test the Monte Carlo swap. +One can use a similar test as previously, but with an imposed chemical +potential *desired_mu*: -.. label:: start_test_9a_class +.. label:: start_test_8a_class .. code-block:: python @@ -141,65 +248,56 @@ Let's test the Monte Carlo swap. import os # Define atom number of each group - nmb_1 = 50 - nmb_2 = 50 # New group for testing swaps + nmb_1= 50 # Define LJ parameters (sigma) - sig_1 = 3 * ureg.angstrom - sig_2 = 4 * ureg.angstrom # Different sigma for group 2 + sig_1 = 3*ureg.angstrom # Define LJ parameters (epsilon) - eps_1 = 0.1 * ureg.kcal / ureg.mol - eps_2 = 0.15 * ureg.kcal / ureg.mol # Different epsilon for group 2 + eps_1 = 0.1*ureg.kcal/ureg.mol # Define atom mass - mss_1 = 10 * ureg.gram / ureg.mol - mss_2 = 12 * ureg.gram / ureg.mol # Different mass for group 2 + mss_1 = 10*ureg.gram/ureg.mol # Define box size - L = 20 * ureg.angstrom + L = 20*ureg.angstrom # Define a cut off - rc = 2.5 * sig_1 + rc = 2.5*sig_1 # Pick the desired temperature - T = 300 * ureg.kelvin + T = 300*ureg.kelvin + # choose the desired_mu + desired_mu = -3*ureg.kcal/ureg.mol # Initialize the prepare object mc = MonteCarlo( - ureg=ureg, + ureg = ureg, maximum_steps=100, thermo_period=10, dumping_period=10, - number_atoms=[nmb_1, nmb_2], # Include two groups of atoms for swap - epsilon=[eps_1, eps_2], # kcal/mol - sigma=[sig_1, sig_2], # A - atom_mass=[mss_1, mss_2], # g/mol - box_dimensions=[L, L, L], # A + number_atoms=[nmb_1], + epsilon=[eps_1], # kcal/mol + sigma=[sig_1], # A + atom_mass=[mss_1], # g/mol + box_dimensions=[L, L, L], # A cut_off=rc, thermo_outputs="Epot-press", - desired_temperature=T, # K + desired_temperature=T, # K neighbor=1, - swap_type=[0, 1] # Enable Monte Carlo swap between groups 1 and 2 + desired_mu = desired_mu, ) - - # Run the Monte Carlo simulation mc.run() # Test function using pytest def test_output_files(): assert os.path.exists("Outputs/dump.mc.lammpstrj"), \ - "Test failed: dump file was not created" + "Test failed: dump file was not created" assert os.path.exists("Outputs/simulation.log"), \ - "Test failed: log file was not created" + "Test failed: log file was not created" print("Test passed") - # Test the swap counters - def test_swap_counters(): - assert mc.successful_swap + mc.failed_swap > 0, \ - "Test failed: No swaps were attempted" - print("Swap test passed") - # If the script is run directly, execute the tests if __name__ == "__main__": import pytest # Run pytest programmatically pytest.main(["-s", __file__]) -.. label:: end_test_9a_class - +.. label:: end_test_8a_class +The evolution of the potential energy as a function of the +number of steps is written in the *Outputs/Epot.dat*. diff --git a/docs/source/index.rst b/docs/source/index.rst index 55aec27..36009c3 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -20,10 +20,10 @@ Goals :align: right :class: only-light -By following *Step-by-Step Molecular Simulations with Python* (MoleSimPy), you will write a simple Python -code containing the most basic functionalities of molecular dynamics and Monte Carlo -simulations. The main goal is to help users understand the basics of molecular simulation -algorithms. +By following *Step-by-Step Molecular Simulations with Python* (MoleSimPy), you +will write a simple Python code containing the most basic functionalities of +molecular dynamics and Monte Carlo simulations. The main goal is to help users +understand the basics of molecular simulation algorithms. In a second part of the project, the scripts will be used for small scientific projects. In particular, some of the earliest results obtained by molecular