Monte Carlo swap¶
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.
The principle of such Monte Carlo swap simulation is the following:
1) We start from a given intial configuration, and measure the potential energy, \(E_\text{pot}^\text{initial}\).
Two particles are chosen randomly, and their respecitive positions are exchanced.
The energy of the system after the move, \(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: \(\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}\).
If \(\Delta E < 0\), then the move is automatically accepted.
If \(\Delta E > 0\), then the move is accepted with a probability given by the Boltzmann factor \(\exp{- \beta \Delta E}\), where \(\beta = 1 / k_\text{B} T\) and \(T\) is the imposed temperature.
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 a method named monte_carlo_swap to the MonteCarlo class:
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)
# 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
Let us initialise swap counter:
class MonteCarlo(Measurements):
def __init__(self,
(...)
self.failed_move = 0
self.successful_swap = 0
self.failed_swap = 0
Complete the __init__ method as follows:
class MonteCarlo(Measurements):
def __init__(self,
(...)
displace_mc = None,
swap_type = [None, None],
and
class MonteCarlo(Measurements):
def __init__(self,
(...)
self.displace_mc = displace_mc
self.swap_type = swap_type
Finally, the monte_carlo_exchange() method must be included in the run:
def run(self):
(...)
self.monte_carlo_move()
self.monte_carlo_swap()
Test the code¶
Let’s test the Monte Carlo swap.
from MonteCarlo import MonteCarlo
from pint import UnitRegistry
ureg = UnitRegistry()
import os
# Define atom number of each group
nmb_1 = 50
nmb_2 = 50 # New group for testing swaps
# Define LJ parameters (sigma)
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_2 = 0.15 * ureg.kcal / ureg.mol # Different epsilon for group 2
# Define atom mass
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
# Define a cut off
rc = 2.5 * sig_1
# Pick the desired temperature
T = 300 * ureg.kelvin
# Initialize the prepare object
mc = MonteCarlo(
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
cut_off=rc,
thermo_outputs="Epot-press",
desired_temperature=T, # K
neighbor=1,
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"
assert os.path.exists("Outputs/simulation.log"), \
"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__])