deploy: b335657

_sources/chapters/chapter3.rst.txt

@@ -183,8 +183,8 @@ neighboring atoms, the simulation becomes more efficient. Add the following
 
 .. code-block:: python
 
-    def update_neighbor_lists(self):
-        if (self.step % self.neighbor == 0):
+    def update_neighbor_lists(self, force_update=False):
+        if (self.step % self.neighbor == 0) | force_update:
             matrix = distances.contact_matrix(self.atoms_positions,
                 cutoff=self.cut_off, #+2,
                 returntype="numpy",
@@ -218,8 +218,8 @@ below).
 
 .. code-block:: python
 
-    def update_cross_coefficients(self):
-        if (self.step % self.neighbor == 0):
+    def update_cross_coefficients(self, force_update=False):
+        if (self.step % self.neighbor == 0) | force_update:
             # Precalculte LJ cross-coefficients
             sigma_ij_list = []
             epsilon_ij_list = []

_sources/chapters/chapter6.rst.txt

@@ -1,7 +1,7 @@
 .. _chapter6-label:
 
-Monte Carlo move
-================
+Monte Carlo dispace
+===================
 
 .. figure:: ../projects/project1/avatar-dm.webp
     :alt: The fluid made of argon atoms and simulated using monte carlo and python.
@@ -15,7 +15,8 @@ Monte Carlo move
     :align: right
     :class: only-light
 
-Here, a *Monte Carlo move* simulation is implemented. The principle of the
+Here, a Monte Carlo simulation is implemented where the only allowed move
+is a dispacement of the particles. The principle of such Monte Carlo
 simulation is the following:
 
 - 1) We start from a given intial configuration, and measure the potential

_sources/chapters/chapter8.rst.txt

@@ -29,39 +29,40 @@ Let us add the following method called *monte_carlo_exchange* to the *MonteCarlo
 .. code-block:: python
 
     def monte_carlo_exchange(self):
-        # 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)
-        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 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)
+            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.successful_insert += 1
+                self.monte_carlo_insert()
             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
+                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:
-                self.failed_delete += 1
+                # 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
 
@@ -237,7 +238,7 @@ Test the code
 One can use a similar test as previously, but with an imposed chemical
 potential *desired_mu*:
 
-.. label:: start_test_8a_class
+.. label:: start_test_9a_class
 
 .. code-block:: python
 
@@ -296,7 +297,7 @@ potential *desired_mu*:
         # 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 are written in the *Outputs/Epot.dat*.

_sources/chapters/chapter9.rst.txt

@@ -0,0 +1,127 @@
+.. _chapter9-label:
+
+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.
+
+.. 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
+            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
+                
+.. label:: end_MonteCarlo_class
+
+Let us initialise swap counter:
+
+.. label:: start_MonteCarlo_class
+
+.. code-block:: python
+
+    class MonteCarlo(Outputs):
+        def __init__(self,
+            (...)
+            self.failed_move = 0
+            self.successful_swap = 0
+            self.failed_swap = 0
+
+.. label:: end_MonteCarlo_class
+
+Complete the *__init__* method as follows:
+
+.. label:: start_MonteCarlo_class
+
+.. code-block:: python
+
+    class MonteCarlo(Outputs):
+        def __init__(self,
+                    (...)
+                    displace_mc = None,
+                    swap_type = [None, None],
+
+.. label:: end_MonteCarlo_class
+
+and
+
+.. label:: start_MonteCarlo_class
+
+.. code-block:: python
+
+    class MonteCarlo(Outputs):
+        def __init__(self,
+            (...)
+            self.displace_mc = displace_mc
+            self.swap_type = swap_type
+
+.. label:: end_MonteCarlo_class
+
+Finally, the *monte_carlo_exchange()* method must be included in the run:
+
+.. label:: start_MonteCarlo_class
+
+.. code-block:: python
+
+    def run(self):
+        (...)
+            self.monte_carlo_move()
+            self.monte_carlo_swap()
+
+.. label:: end_MonteCarlo_class

_sources/index.rst.txt

@@ -58,6 +58,7 @@ Learn Molecular Simulations with Python can be followed by anyone with a compute
     chapters/chapter6
     chapters/chapter7
     chapters/chapter8
+    chapters/chapter9
 
 .. toctree::
     :maxdepth: 2

chapters/chapter1.html

@@ -209,9 +209,10 @@
 <li class="toctree-l1"><a class="reference internal" href="chapter3.html">Initialize the simulation</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter4.html">Minimize the energy</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter5.html">Output the simulation state</a></li>
-<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo move</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo dispace</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter7.html">Pressure Measurement</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter8.html">Monte Carlo insert</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter9.html">Monte Carlo swap</a></li>
 </ul>
 <p class="caption" role="heading"><span class="caption-text">Illustration</span></p>
 <ul>
@@ -479,7 +480,7 @@ <h2>Create the Classes<a class="headerlink" href="#create-the-classes" title="Li
 <p>The <em>ignore warnings</em> commands are optional; they were added to avoid the
 annoying message “<em>RuntimeWarning: overflow encountered in exp</em>” that is sometimes
 triggered by the exponential function of <em>acceptation_probability</em> (see the
-<a class="reference internal" href="chapter6.html#chapter6-label"><span class="std std-ref">Monte Carlo move</span></a> chapter).</p>
+<a class="reference internal" href="chapter6.html#chapter6-label"><span class="std std-ref">Monte Carlo dispace</span></a> chapter).</p>
 <p>Finally, within the <em>MolecularDynamics.py</em> file, copy the following lines:</p>
 <div class="highlight-python notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
 <span class="kn">from</span> <span class="nn">Measurements</span> <span class="kn">import</span> <span class="n">Measurements</span>

chapters/chapter2.html

@@ -209,9 +209,10 @@
 <li class="toctree-l1"><a class="reference internal" href="chapter3.html">Initialize the simulation</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter4.html">Minimize the energy</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter5.html">Output the simulation state</a></li>
-<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo move</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo dispace</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter7.html">Pressure Measurement</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter8.html">Monte Carlo insert</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter9.html">Monte Carlo swap</a></li>
 </ul>
 <p class="caption" role="heading"><span class="caption-text">Illustration</span></p>
 <ul>

chapters/chapter3.html

@@ -209,9 +209,10 @@
 <li class="toctree-l1 current current-page"><a class="current reference internal" href="#">Initialize the simulation</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter4.html">Minimize the energy</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter5.html">Output the simulation state</a></li>
-<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo move</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo dispace</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter7.html">Pressure Measurement</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter8.html">Monte Carlo insert</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter9.html">Monte Carlo swap</a></li>
 </ul>
 <p class="caption" role="heading"><span class="caption-text">Illustration</span></p>
 <ul>
@@ -389,8 +390,8 @@ <h2>Build Neighbor Lists<a class="headerlink" href="#build-neighbor-lists" title
 to save computational time. By focusing only on interactions between
 neighboring atoms, the simulation becomes more efficient. Add the following
 <em>update_neighbor_lists()</em> method to the <em>Utilities</em> class:</p>
-<div class="highlight-python notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">update_neighbor_lists</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
-    <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">step</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">neighbor</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
+<div class="highlight-python notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">update_neighbor_lists</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">force_update</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
+    <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">step</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">neighbor</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="o">|</span> <span class="n">force_update</span><span class="p">:</span>
         <span class="n">matrix</span> <span class="o">=</span> <span class="n">distances</span><span class="o">.</span><span class="n">contact_matrix</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">atoms_positions</span><span class="p">,</span>
             <span class="n">cutoff</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">cut_off</span><span class="p">,</span> <span class="c1">#+2,</span>
             <span class="n">returntype</span><span class="o">=</span><span class="s2">&quot;numpy&quot;</span><span class="p">,</span>
@@ -416,8 +417,8 @@ <h2>Update Cross Coefficients<a class="headerlink" href="#update-cross-coefficie
 <p>As the neighbor lists are being built, let us also pre-calculate the cross
 coefficients. This will make the force calculation more efficient (see
 below).</p>
-<div class="highlight-python notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">update_cross_coefficients</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
-    <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">step</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">neighbor</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
+<div class="highlight-python notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">update_cross_coefficients</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">force_update</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
+    <span class="k">if</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">step</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">neighbor</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="o">|</span> <span class="n">force_update</span><span class="p">:</span>
         <span class="c1"># Precalculte LJ cross-coefficients</span>
         <span class="n">sigma_ij_list</span> <span class="o">=</span> <span class="p">[]</span>
         <span class="n">epsilon_ij_list</span> <span class="o">=</span> <span class="p">[]</span>

chapters/chapter4.html

@@ -209,9 +209,10 @@
 <li class="toctree-l1"><a class="reference internal" href="chapter3.html">Initialize the simulation</a></li>
 <li class="toctree-l1 current current-page"><a class="current reference internal" href="#">Minimize the energy</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter5.html">Output the simulation state</a></li>
-<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo move</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter6.html">Monte Carlo dispace</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter7.html">Pressure Measurement</a></li>
 <li class="toctree-l1"><a class="reference internal" href="chapter8.html">Monte Carlo insert</a></li>
+<li class="toctree-l1"><a class="reference internal" href="chapter9.html">Monte Carlo swap</a></li>
 </ul>
 <p class="caption" role="heading"><span class="caption-text">Illustration</span></p>
 <ul>