diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst
index 92bbed6..632d616 100644
--- a/docs/source/chapters/chapter1.rst
+++ b/docs/source/chapters/chapter1.rst
@@ -1,8 +1,13 @@
.. _chapter1-label:
-Start coding
+Start Coding
============
+Let's start with the Python script. During this first chapter, some Python
+files will be created and filled in a minimal fashion. At the end of this
+chapter, a small test will be set up to ensure that the files were correctly
+created.
+
Presentation
------------
@@ -34,41 +39,47 @@ containing either Python functions or classes:
* - File Name
- Content
- * - *Prepare.py*
+ * - Prepare.py
- *Prepare* class: Methods for preparing the non-dimensionalization of the
units
- * - *Utilities.py*
- - *Utilities* class: General-purpose methods, inherited by all other classes
- * - *InitializeSimulation.py*
+ * - Utilities.py
+ - *Utilities* class: General-purpose methods, inherited by all other
+ classes
+ * - InitializeSimulation.py
- *InitializeSimulation* class: Methods necessary to set up the system and
prepare the simulation, inherited by all the classes below
- * - *MinimizeEnergy.py*
+ * - MinimizeEnergy.py
- *MinimizeEnergy* class: Methods for performing energy minimization
- * - *MonteCarlo.py*
+ * - MonteCarlo.py
- *MonteCarlo* class: Methods for performing Monte Carlo simulations in
different ensembles (e.g., Grand Canonical, Canonical)
- * - *MolecularDynamics.py*
+ * - MolecularDynamics.py
- *MolecularDynamics* class: Methods for performing molecular dynamics in
different ensembles (NVE, NPT, NVT)
- * - *measurements.py*
- - Functions for performing specific measurements on the system
+ * - Measurements.py
+ - *Measurements* class: Methods for for performing specific measurements on the system
* - *potentials.py*
- Functions for calculating the potentials and forces between atoms
- * - *tools.py*
+ * - logger.py
- Functions for outputting data into text files
+ * - dumper.py
+ - Functions for outputting data into trajectory files for visualization
+ * - reader.py
+ - Functions for importing data from text files
+Some of these files are created in this chapter; others will be created later
+on. All of these files must be created within the same folder.
-Potential for inter-atomic interaction
+Potential for Inter-Atomic Interaction
--------------------------------------
-In molecular simulations, potential functions are used to mimic the interaction
-between atoms. Although more complicated options exist, potentials are usually
+In molecular simulations, potential functions are used to model the interaction
+between atoms. Although more complex options exist, potentials are usually
defined as functions of the distance :math:`r` between two atoms.
-Within a dedicated folder, create the first file named *potentials.py*. This
-file will contain a function called *potentials*. Two types of potential can
-be returned by this function: the Lennard-Jones potential (LJ), and the
-hard-sphere potential.
+Create a file named *potentials.py*. This file will contain a function called
+*potentials*. For now, the only potential that can be returned by this function
+is the Lennard-Jones (LJ) potential, but this may change in the future.
Copy the following lines into *potentials.py*:
@@ -76,52 +87,41 @@ Copy the following lines into *potentials.py*:
.. code-block:: python
- import numpy as np
-
- def potentials(potential_type, epsilon, sigma, r, derivative=False):
- if potential_type == "Lennard-Jones":
- if derivative:
- return 48 * epsilon * ((sigma / r) ** 12 - 0.5 * (sigma / r) ** 6) / r
- else:
- return 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
- elif potential_type == "Hard-Sphere":
- if derivative:
- # Derivative is not defined for Hard-Sphere potential.
- # --> return 0
- return np.zeros(len(r))
- else:
- return np.where(r > sigma, 0, 1000)
+ def potentials(epsilon, sigma, r, derivative=False):
+ if derivative:
+ return 48 * epsilon * ((sigma / r) ** 12 - 0.5 * (sigma / r) ** 6) / r
else:
- raise ValueError(f"Unknown potential type: {potential_type}")
+ return 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
.. label:: end_potentials_class
-The hard-sphere potential either returns a value of 0 when the distance between
-the two particles is larger than the parameter, :math:`r > \sigma`, or 1000 when
-:math:`r < \sigma`. The value of *1000* was chosen to be large enough to ensure
-that any Monte Carlo move that would lead to the two particles to overlap will
-be rejected.
-
-In the case of the LJ potential, depending on the value of the optional
-argument *derivative*, which can be either *False* or *True*, the *LJ_potential*
-function will return the force:
+Depending on the value of the optional argument *derivative*, which can be
+either *False* or *True*, this function returns the derivative of the potential,
+i.e., the force, :math:`F_\text{LJ} = - \mathrm{d} U_\text{LJ} / \mathrm{d} r`:
.. math::
- F_\text{LJ} = 48 \dfrac{\epsilon}{r} \left[ \left( \frac{\sigma}{r} \right)^{12} - \frac{1}{2} \left( \frac{\sigma}{r} \right)^6 \right],
+ F_\text{LJ} = 48 \dfrac{\epsilon}{r} \left[ \left( \frac{\sigma}{r} \right)^{12}
+ - \frac{1}{2} \left( \frac{\sigma}{r} \right)^6 \right], ~ \text{for} ~ r < r_\text{c},
or the potential energy:
.. math::
- U_\text{LJ} = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right].
+ U_\text{LJ} = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12}
+ - \left( \frac{\sigma}{r} \right)^6 \right], ~ \text{for} ~ r < r_\text{c}.
+
+Here, :math:`\sigma` is the distance at which the potential :math:`U_\text{LJ}`
+is zero, :math:`\epsilon` is the depth of the potential well, and
+:math:`r_\text{c}` is a cutoff distance. For :math:`r > r_\text{c}`,
+:math:`U_\text{LJ} = 0` and :math:`F_\text{LJ} = 0`.
Create the Classes
------------------
-Let's create the files with the minimal information about the classes and
-their inheritance. The classes will be developed progressively in the
-following chapters.
+Let's create the files with some minimal details about the classes and their
+inheritance. The classes will be developed progressively in the following
+chapters.
The first class is the *Prepare* class, which will be used for the
nondimensionalization of the parameters. In the same folder as *potentials.py*,
@@ -157,8 +157,8 @@ copy the following lines:
.. label:: end_Utilities_class
-The line *from potentials import LJ_potential* is used to import the
-*LJ_potential* function.
+The line *from potentials import potentials* is used to import the
+previously created *potentials* function.
Within the *InitializeSimulation.py* file, copy the following lines:
@@ -166,11 +166,13 @@ Within the *InitializeSimulation.py* file, copy the following lines:
.. code-block:: python
+ import os
import numpy as np
from Prepare import Prepare
+ from Utilities import Utilities
- class InitializeSimulation(Prepare):
+ class InitializeSimulation(Prepare, Utilities):
def __init__(self,
*args,
**kwargs,
@@ -180,7 +182,14 @@ Within the *InitializeSimulation.py* file, copy the following lines:
.. label:: end_InitializeSimulation_class
The *InitializeSimulation* class inherits from the previously created
-*Prepare* class. Additionally, we anticipate that *NumPy* will be required.
+*Prepare* and *Utilities* classes. Additionally, we anticipate that |NumPy|
+will be required :cite:`harris2020array`. We also anticipate that the *os*
+module, which provides a way to interact with the operating system, will
+be required :cite:`Rossum2009Python3`.
+
+.. |NumPy| raw:: html
+
+ NumPy
Within the *Measurements.py* file, copy the following lines:
@@ -188,11 +197,11 @@ Within the *Measurements.py* file, copy the following lines:
.. code-block:: python
+ import numpy as np
from InitializeSimulation import InitializeSimulation
- from Utilities import Utilities
- class Measurements(InitializeSimulation, Utilities):
+ class Measurements(InitializeSimulation):
def __init__(self,
*args,
**kwargs):
@@ -200,19 +209,23 @@ Within the *Measurements.py* file, copy the following lines:
.. label:: end_Measurements_class
-The *Measurements* class inherits both the *InitializeSimulation* and
-*Utilities* classes.
+The *Measurements* class inherits from *InitializeSimulation* (and thus
+also inherits from the *Prepare* and *Utilities* classes).
+
+Finally, let us create the three remaining classes: *MinimizeEnergy*,
+*MonteCarlo*, and *MolecularDynamics*. Each of these classes inherits
+from the *Measurements* class (and thus also from the *Prepare*, *Utilities*,
+and *InitializeSimulation* classes).
-Finally, let us create the three remaining classes, named *MinimizeEnergy*,
-*MonteCarlo*, and *MolecularDynamics*. Each of these three classes inherits
-from the *Measurements* class, and thus from the classes inherited by
-*Measurements*. Within the *MinimizeEnergy.py* file, copy the following lines:
+Within the *MinimizeEnergy.py* file, copy the following lines:
.. label:: start_MinimizeEnergy_class
.. code-block:: python
from Measurements import Measurements
+ import numpy as np
+ import copy
import os
@@ -224,8 +237,8 @@ from the *Measurements* class, and thus from the classes inherited by
.. label:: end_MinimizeEnergy_class
-We anticipate that the *os* module, which provides a way to interact with the
-operating system, will be required :cite:`Rossum2009Python3`.
+The *copy* library, which provides functions to create shallow or deep copies of
+objects, is imported, along with *NumPy* and *os*.
Within the *MonteCarlo.py* file, copy the following lines:
@@ -233,10 +246,8 @@ Within the *MonteCarlo.py* file, copy the following lines:
.. code-block:: python
- from scipy import constants as cst
import numpy as np
import copy
- import os
from Measurements import Measurements
import warnings
@@ -251,12 +262,10 @@ Within the *MonteCarlo.py* file, copy the following lines:
.. label:: end_MonteCarlo_class
-Several libraries were imported, namely *Constants* from *SciPy*, *NumPy*, *copy*
-and *os*.
-
-The *warnings* was placed to avoid the anoying message "*RuntimeWarning: overflow
-encountered in exp*" that is sometimes triggered by the exponential of the
-*acceptation_probability* (see :ref:`chapter6-label`).
+The *ignore warnings* commands are optional; they were added to avoid the
+annoying message "*RuntimeWarning: overflow encountered in exp*" that is sometimes
+triggered by the exponential function of *acceptation_probability* (see the
+:ref:`chapter6-label` chapter).
Finally, within the *MolecularDynamics.py* file, copy the following lines:
@@ -294,12 +303,14 @@ and copy the following lines into it:
# Make sure that MonteCarlo correctly inherits from Utilities
def test_montecarlo_inherits_from_utilities():
- assert issubclass(MonteCarlo, Utilities), "MonteCarlo should inherit from Utilities"
+ assert issubclass(MonteCarlo, Utilities), \
+ "MonteCarlo should inherit from Utilities"
print("MonteCarlo correctly inherits from Utilities")
# Make sure that Utilities does not inherit from MonteCarlo
def test_utilities_does_not_inherit_from_montecarlo():
- assert not issubclass(Utilities, MonteCarlo), "Utilities should not inherit from MonteCarlo"
+ assert not issubclass(Utilities, MonteCarlo), \
+ "Utilities should not inherit from MonteCarlo"
print("Utilities does not inherit from MonteCarlo, as expected")
# In the script is launched with Python, call Pytest
@@ -317,15 +328,15 @@ any *AssertionError*:
Utilities does not inherit from MonteCarlo, as expected
MonteCarlo correctly inherits from Utilities
-Alternatively, this test can also be launched using Pytest by typing in a terminal:
+Alternatively, this test can also be launched using *Pytest* by typing in a terminal:
.. code-block:: bash
pytest .
-We can also test that calling the *__init__*
-method of the *MonteCarlo* class does not return any error. In new Python file
-called *test_1b.py*, copy the following lines:
+We can also test that calling the *__init__* method of the *MonteCarlo* class
+does not return any error. In new Python file called *test_1b.py*, copy the
+following lines:
.. label:: start_test_1b_class
diff --git a/docs/source/chapters/chapter2.rst b/docs/source/chapters/chapter2.rst
index 5752c44..3b1ccdd 100644
--- a/docs/source/chapters/chapter2.rst
+++ b/docs/source/chapters/chapter2.rst
@@ -29,7 +29,7 @@ The *real* unit system follows the conventions outlined in the |lammps-unit-syst
- Forces are in (kcal/mol)/Ångström,
- Temperature is in Kelvin,
- Pressure is in atmospheres,
-- Density is in g/cm\ :sup:`3` (in 3D).
+- Density is in g/cm\ :sup:`3`.
.. |lammps-unit-systems| raw:: html
@@ -56,12 +56,8 @@ where :math:`k_\text{B}` is the Boltzmann constant.
Start coding
------------
-Let's fill in the previously created class named Prepare. To facilitate unit
-conversion, we will import |NumPy| and the constants module from |SciPy|.
-
-.. |NumPy| raw:: html
-
- NumPy
+Let's fill in the previously created class named *Prepare*. To facilitate unit
+conversion, we will import NumPy and the constants module from |SciPy|.
.. |SciPy| raw:: html
@@ -78,7 +74,7 @@ In the file named *Prepare.py*, add the following lines:
.. label:: end_Prepare_class
-Four parameters are provided to the *Prepare* class:
+Four atom parameters are provided to the *Prepare* class:
- the atom masses :math:`m`,
- the LJ parameters :math:`\sigma` and :math:`\epsilon`,
@@ -95,29 +91,31 @@ Modify the *Prepare* class as follows:
class Prepare:
def __init__(self,
- number_atoms=[10], # List - no unit
- epsilon=[0.1], # List - Kcal/mol
- sigma=[3], # List - Angstrom
- atom_mass=[10], # List - g/mol
- potential_type="Lennard-Jones",
+ ureg, # Pint unit registry
+ number_atoms, # List - no unit
+ epsilon, # List - Kcal/mol
+ sigma, # List - Angstrom
+ atom_mass, # List - g/mol
*args,
**kwargs):
+ self.ureg = ureg
self.number_atoms = number_atoms
self.epsilon = epsilon
self.sigma = sigma
self.atom_mass = atom_mass
- self.potential_type = potential_type
super().__init__(*args, **kwargs)
.. label:: end_Prepare_class
-Here, the four lists *number_atoms* :math:`N`, *epsilon* :math:`\epsilon`,
-*sigma* :math:`\sigma`, and *atom_mass* :math:`m` are given default values of
-:math:`10`, :math:`0.1~\text{[Kcal/mol]}`, :math:`3~\text{[Å]}`, and
-:math:`10~\text{[g/mol]}`, respectively.
+Here, *number_atoms* :math:`N`, *epsilon* :math:`\epsilon`,
+*sigma* :math:`\sigma`, and *atom_mass* :math:`m` must be provided as lists
+where the elements have no units, kcal/mol, angstrom, and g/mol units,
+respectively. The units will be enforced with the |Pint| unit registry, *ureg*,
+which must also be provided as a parameter (more details later).
+
+.. |Pint| raw:: html
-The type of potential is also specified, with Lennard-Jones being closen as
-the default option.
+ Pint
All the parameters are assigned to *self*, allowing other methods to access
them. The *args* and *kwargs* parameters are used to accept an arbitrary number
@@ -126,9 +124,9 @@ of positional and keyword arguments, respectively.
Calculate LJ units prefactors
-----------------------------
-Within the *Prepare* class, let us create a method called *calculate_LJunits_prefactors*
-that will be used to calculate the prefactors necessary to convert units from the *real*
-unit system to the *LJ* unit system:
+Within the *Prepare* class, create a method called *calculate_LJunits_prefactors*
+that will be used to calculate the prefactors necessary to convert units from the
+*real* unit system to the *LJ* unit system:
.. label:: start_Prepare_class
@@ -136,37 +134,47 @@ unit system to the *LJ* unit system:
def calculate_LJunits_prefactors(self):
"""Calculate the Lennard-Jones units prefactors."""
- # Define the reference distance, energy, and mass
- self.reference_distance = self.sigma[0] # Angstrom
- self.reference_energy = self.epsilon[0] # kcal/mol
- self.reference_mass = self.atom_mass[0] # g/mol
-
- # Calculate the prefactor for the time
- mass_kg = self.atom_mass[0]/cst.kilo/cst.Avogadro # kg
- epsilon_J = self.epsilon[0]*cst.calorie*cst.kilo/cst.Avogadro # J
- sigma_m = self.sigma[0]*cst.angstrom # m
- time_s = np.sqrt(mass_kg*sigma_m**2/epsilon_J) # s
- self.reference_time = time_s / cst.femto # fs
-
- # Calculate the prefactor for the temperature
+ # First define Boltzmann and Avogadro constants
kB = cst.Boltzmann*cst.Avogadro/cst.calorie/cst.kilo # kcal/mol/K
- self.reference_temperature = self.epsilon[0]/kB # K
-
- # Calculate the prefactor for the pressure
- pressure_pa = epsilon_J/sigma_m**3 # Pa
- self.reference_pressure = pressure_pa/cst.atm # atm
+ kB *= self.ureg.kcal/self.ureg.mol/self.ureg.kelvin
+ Na = cst.Avogadro/self.ureg.mol
+ # Define the reference distance, energy, and mass
+ self.ref_length = self.sigma[0] # Angstrom
+ self.ref_energy = self.epsilon[0] # kcal/mol
+ self.ref_mass = self.atom_mass[0] # g/mol
+ # Optional: assert that units were correctly provided by users
+ assert self.ref_length.units == self.ureg.angstrom, \
+ f"Error: Provided sigma has wrong units, should be angstrom"
+ assert self.ref_energy.units == self.ureg.kcal/self.ureg.mol, \
+ f"Error: Provided epsilon has wrong units, should be kcal/mol"
+ assert self.ref_mass.units == self.ureg.g/self.ureg.mol, \
+ f"Error: Provided mass has wrong units, should be g/mol"
+ # Calculate the prefactor for the time (in femtosecond)
+ self.ref_time = np.sqrt(self.ref_mass \
+ *self.ref_length**2/self.ref_energy).to(self.ureg.femtosecond)
+ # Calculate the prefactor for the temperature (in Kelvin)
+ self.ref_temperature = self.ref_energy/kB # Kelvin
+ # Calculate the prefactor for the pressure (in Atmosphere)
+ self.ref_pressure = (self.ref_energy \
+ /self.ref_length**3/Na).to(self.ureg.atmosphere)
+ # Group all the reference quantities into a list for practicality
+ self.ref_quantities = [self.ref_length, self.ref_energy,
+ self.ref_mass, self.ref_time, self.ref_pressure, self.ref_temperature]
+ self.ref_units = [ref.units for ref in self.ref_quantities]
.. label:: end_Prepare_class
-This method defines the reference distance as the first element in the
-*sigma* list, i.e., :math:`\sigma_{11}`. Therefore, atoms of type one will
-always be used for the normalization. Similarly, the first element
-in the *epsilon* list (:math:`\epsilon_{11}`) is used as the reference energy,
-and the first element in the *atom_mass* list (:math:`m_1`) is used as the
-reference mass. Then, the reference_time in femtoseconds is calculated
-as :math:`\sqrt{m_1 \sigma_{11}^2 / \epsilon_{11}}`, the reference temperature
-in Kelvin as :math:`\epsilon_{11} / k_\text{B}`, and the reference_pressure
-in atmospheres is calculated as :math:`\epsilon_{11}/\sigma_{11}^3`.
+This method defines the reference length as the first element in the
+*sigma* list, i.e., :math:`\sigma_{11}`. Therefore, atoms of type 1 will
+always be used for normalization. Similarly, the first element in the
+*epsilon* list (:math:`\epsilon_{11}`) is used as the reference energy, and
+the first element in the *atom_mass* list (:math:`m_1`) is used as the
+reference mass.
+
+The reference time in femtoseconds is then calculated as
+:math:`\sqrt{m_1 \sigma_{11}^2 / \epsilon_{11}}`, the reference
+temperature in Kelvin as :math:`\epsilon_{11} / k_\text{B}`, and the
+reference pressure in atmospheres as :math:`\epsilon_{11} / \sigma_{11}^3`.
Finally, let us ensure that the *calculate_LJunits_prefactors* method is
called systematically by adding the following line to the *__init__()* method:
@@ -192,33 +200,48 @@ Let us take advantage of the calculated reference values and normalize the
three inputs of the *Prepare* class that have physical dimensions, i.e.,
*epsilon*, *sigma*, and *atom_mass*.
-Create a new method called *nondimensionalize_units_0* within the *Prepare*
+Create a new method called *nondimensionalize_units* within the *Prepare*
class:
.. label:: start_Prepare_class
.. code-block:: python
- def nondimensionalize_units_0(self):
- # Normalize LJ properties
- epsilon, sigma, atom_mass = [], [], []
- for e0, s0, m0 in zip(self.epsilon, self.sigma, self.atom_mass):
- epsilon.append(e0/self.reference_energy)
- sigma.append(s0/self.reference_distance)
- atom_mass.append(m0/self.reference_mass)
- self.epsilon = epsilon
- self.sigma = sigma
- self.atom_mass = atom_mass
+ def nondimensionalize_units(self, quantities_to_normalise):
+ for name in quantities_to_normalise:
+ quantity = getattr(self, name) # Get the attribute by name
+ if isinstance(quantity, list):
+ for i, element in enumerate(quantity):
+ assert element.units in self.ref_units, \
+ f"Error: Units not part of the reference units"
+ ref_value = self.ref_quantities[self.ref_units.index(element.units)]
+ quantity[i] = element/ref_value
+ assert quantity[i].units == self.ureg.dimensionless, \
+ f"Error: Quantities are not properly nondimensionalized"
+ quantity[i] = quantity[i].magnitude # get rid of ureg
+ setattr(self, name, quantity)
+ else:
+ if quantity is not None:
+ assert np.shape(quantity) == (), \
+ f"Error: The quantity is a list or an array"
+ assert quantity.units in self.ref_units, \
+ f"Error: Units not part of the reference units"
+ ref_value = self.ref_quantities[self.ref_units.index(quantity.units)]
+ quantity = quantity/ref_value
+ assert quantity.units == self.ureg.dimensionless, \
+ f"Error: Quantities are not properly nondimensionalized"
+ quantity = quantity.magnitude # get rid of ureg
+ setattr(self, name, quantity)
.. label:: end_Prepare_class
-The index *0* is used to differentiate this method from other methods that
-will be used to nondimensionalize units in future classes. We anticipate that
-*epsilon*, *sigma*, and *atom_mass* may contain more than one element, so
-each element is normalized with the corresponding reference value. The
-*zip()* function allows us to loop over all three lists at once.
+When a *quantities_to_normalise* list containing parameter names is provided
+to the *nondimensionalize_units* method, a loop is performed over all the
+quantities. The value of each quantity is extracted using *getattr*. The units
+of the quantity of interest are then detected and normalized by the appropriate
+reference quantities defined by *calculate_LJunits_prefactors*.
-Let us also call the *nondimensionalize_units_0* from the *__init__()* method
+Let us also call the *nondimensionalize_units* from the *__init__()* method
of the *Prepare* class:
.. label:: start_Prepare_class
@@ -228,26 +251,29 @@ of the *Prepare* class:
def __init__(self,
(...)
self.calculate_LJunits_prefactors()
- self.nondimensionalize_units_0()
+ self.nondimensionalize_units(["epsilon", "sigma", "atom_mass"])
.. label:: end_Prepare_class
-Identify atom properties
+Here, the *epsilon*, *sigma*, and *atom_mass* parameters will be
+nondimensionalized.
+
+Identify Atom Properties
------------------------
Anticipating the future use of multiple atom types, where each type will be
associated with its own :math:`\sigma`, :math:`\epsilon`, and :math:`m`, let
us create arrays containing the properties of each atom in the simulation. For
-instance, in the case of a simulation with two atoms of type 1 and three atoms
-of type 2, the corresponding *atoms_sigma* array will be:
+instance, in a simulation with two atoms of type 1 and three atoms of type 2,
+the corresponding *atoms_sigma* array will be:
.. math::
\text{atoms_sigma} = [\sigma_{11}, \sigma_{11}, \sigma_{22}, \sigma_{22}, \sigma_{22}]
where :math:`\sigma_{11}` and :math:`\sigma_{22}` are the sigma values for
-atoms of type 1 and 2, respectively. The *atoms_sigma* array will allow for
-future calculations of force.
+atoms of type 1 and 2, respectively. The *atoms_sigma* array will facilitate
+future calculations of forces.
Create a new method called *identify_atom_properties*, and place it
within the *Prepare* class:
@@ -257,8 +283,7 @@ within the *Prepare* class:
.. code-block:: python
def identify_atom_properties(self):
- """Identify the atom properties for each atom."""
- self.total_number_atoms = np.sum(self.number_atoms)
+ """Identify the properties for each atom."""
atoms_sigma = []
atoms_epsilon = []
atoms_mass = []
@@ -288,7 +313,7 @@ Let us call the *identify_atom_properties* from the *__init__()* method:
def __init__(self,
(...)
- self.nondimensionalize_units_0()
+ self.nondimensionalize_units(["epsilon", "sigma", "atom_mass"])
self.identify_atom_properties()
.. label:: end_Prepare_class
@@ -296,9 +321,10 @@ Let us call the *identify_atom_properties* from the *__init__()* method:
Test the code
-------------
-Let us test the *Prepare* class to make sure that it does what is expected.
-Just like in the previous example, let us initiates a system with 2 atoms of
-type 1, and 3 atoms of type 2:
+Let's test the *Prepare* class to make sure that it does what is expected.
+Here, a system containing 2 atoms of type 1, and 3 atoms of type 2 is
+prepared. LJs parameters and masses for each groups are also defined, and given
+physical units thanks to the *UnitRegistry* of Pint.
.. label:: start_test_2a_class
@@ -306,20 +332,33 @@ type 1, and 3 atoms of type 2:
import numpy as np
from Prepare import Prepare
-
- # Initialize the Prepare object
+ from pint import UnitRegistry
+ ureg = UnitRegistry()
+
+ # Define atom number of each group
+ nmb_1, nmb_2= [2, 3]
+ # Define LJ parameters (sigma)
+ sig_1, sig_2 = [3, 4]*ureg.angstrom
+ # Define LJ parameters (epsilon)
+ eps_1, eps_2 = [0.2, 0.4]*ureg.kcal/ureg.mol
+ # Define atom mass
+ mss_1, mss_2 = [10, 20]*ureg.gram/ureg.mol
+
+ # Initialize the prepare object
prep = Prepare(
- number_atoms=[2, 3],
- epsilon=[0.2, 0.4], # kcal/mol
- sigma=[3, 4], # A
- atom_mass=[10, 20], # g/mol
+ ureg = ureg,
+ number_atoms=[nmb_1, nmb_2],
+ epsilon=[eps_1, eps_2], # kcal/mol
+ sigma=[sig_1, sig_2], # A
+ atom_mass=[mss_1, mss_2], # g/mol
)
# Test function using pytest
def test_atoms_epsilon():
expected = np.array([1., 1., 2., 2., 2.])
result = prep.atoms_epsilon
- assert np.array_equal(result, expected), f"Test failed: {result} != {expected}"
+ assert np.array_equal(result, expected), \
+ f"Test failed: {result} != {expected}"
print("Test passed")
# In the script is launched with Python, call Pytest
@@ -329,5 +368,6 @@ type 1, and 3 atoms of type 2:
.. label:: end_test_2a_class
-This test assert that the generated *atoms_epsilon* array is consistent with
-its expected value (see the previous paragraphs).
+This test asserts that the generated *atoms_epsilon* array is consistent
+with its expected value (see the previous paragraphs). Note that the expected
+values are in LJ units, i.e., they have been nondimensionalized by :math:`\epsilon_1`.
diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst
index 4731e76..72a6bf6 100644
--- a/docs/source/chapters/chapter3.rst
+++ b/docs/source/chapters/chapter3.rst
@@ -15,105 +15,54 @@ Initialize the simulation
:align: right
:class: only-light
-Here, the *InitializeSimulation* class is created. This class is used to
-prepare the simulation box and populate the box randomly with atoms.
-Let us improve the previously created *InitializeSimulation* class:
+Here, the *InitializeSimulation* class is completed. This class is used to
+prepare the simulation box and populate it randomly with atoms. Improve the
+previously created *InitializeSimulation* class as follows:
.. label:: start_InitializeSimulation_class
.. code-block:: python
- class InitializeSimulation(Prepare):
+ class InitializeSimulation(Prepare, Utilities):
def __init__(self,
- box_dimensions=[10, 10, 10], # List - Angstroms
- seed=None, # Int
+ box_dimensions, # List - Angstroms
+ cut_off, # Angstroms
initial_positions=None, # Array - Angstroms
- thermo_period=None,
- dumping_period=None,
+ neighbor=1, # Integer
*args,
**kwargs,
):
super().__init__(*args, **kwargs)
self.box_dimensions = box_dimensions
- # If a box dimension was entered as 0, dimensions will be 2
- self.dimensions = len(list(filter(lambda x: x > 0, box_dimensions)))
- self.seed = seed
+ self.cut_off = cut_off
+ self.neighbor = neighbor
+ self.step = 0 # initialize simulation step
self.initial_positions = initial_positions
- self.thermo_period = thermo_period
- self.dumping_period = dumping_period
.. label:: end_InitializeSimulation_class
Several parameters are provided to the *InitializeSimulation* class:
-- The first parameter is the box dimensions, which is a list with a length
- corresponding to the dimension of the system. Each element of the list
- corresponds to a dimension of the box in Ångström in the x, y, and z
- directions, respectively.
-- Optionally, a seed can be provided as an integer. If the seed is given
- by the user, the random number generator will always produce the same
- displacements.
+- The box dimensions, which is a list with 3 elements. Each element of
+ the list corresponds to a dimension of the box in Ångström in the :math:`x`,
+ :math:`y`, and :math:`z`` directions, respectively.
+- The cutoff :math:`r_\text{c}` for the potential interaction.
- Optionally, initial positions for the atoms can be provided as an array
of length corresponding to the number of atoms. If *initial_positions*
- is left equal to *None*, positions will be randomly assigned to the
- atoms (see below).
-- Optionally, a thermo period and a dumping_period can be provided to
- control the outputs from the simulation (it will be implemented
- in :ref:`chapter5-label`).
+ is set to *None*, positions will be randomly assigned to the atoms (see
+ below).
+- The *neighbor* parameter determines the interval between recalculations
+ of the neighbor lists. By default, a value of 1 is used, meaning that the
+ neighbor list will be reconstructed every step. In certain cases, this
+ number may be increased to reduce computation time.
-Finally, the *dimensions* of the system are calculated as the length of
-*box_dimensions*.
-
-Set a random seed
------------------
-
-Providing a *seed* allows for reproducible simulations. When a seed is
-specified, the simulation will produce the exact same results each time it
-is run, including identical atom positions and velocities. This can be
-particularly useful for debugging purposes.
-
-Add the following *if* condition to the *__init__()* method:
-
-.. label:: start_InitializeSimulation_class
-
-.. code-block:: python
-
- def __init__(self,
- (...)
- self.initial_positions = initial_positions
- if self.seed is not None:
- np.random.seed(self.seed)
-
-.. label:: end_InitializeSimulation_class
-
-If a *seed* is provided, it is used to initialize the *random.seed()* function
-from *NumPy*. If *seed* is set to its default value of *None*, the simulation
-will proceed with randomization.
+Note that the simulation step is initialized to 0.
Nondimensionalize units
-----------------------
Just like we did in :ref:`chapter2-label`, let us nondimensionalize the provided
-parameters, here the *box_dimensions* as well as the *initial_positions*:
-
-.. label:: start_InitializeSimulation_class
-
-.. code-block:: python
-
- def nondimensionalize_units_1(self):
- """Use LJ prefactors to convert units into non-dimensional."""
- # Normalize box dimensions
- box_dimensions = []
- for L in self.box_dimensions:
- box_dimensions.append(L/self.reference_distance)
- self.box_dimensions = box_dimensions # errase the previously defined box_dimensions
- # Normalize the box dimensions
- if self.initial_positions is not None:
- self.initial_positions = self.initial_positions/self.reference_distance
-
-.. label:: end_InitializeSimulation_class
-
-Let us call the *nondimensionalize_units_1* method from the *__init__* class:
+parameters, here the *box_dimensions*, the *cut_off*, as well as the *initial_positions*:
.. label:: start_InitializeSimulation_class
@@ -121,31 +70,27 @@ Let us call the *nondimensionalize_units_1* method from the *__init__* class:
def __init__(self,
(...)
- if self.seed is not None:
- np.random.seed(self.seed)
- self.nondimensionalize_units_1()
+ self.initial_positions = initial_positions
+ self.nondimensionalize_units(["box_dimensions", "cut_off",
+ "initial_positions"])
.. label:: end_InitializeSimulation_class
Define the box
--------------
-Let us define the simulation box using the values from *box_dimensions*. Add the following
-method to the *InitializeSimulation* class:
+Let us define the simulation box using the values from *box_dimensions*. Add the
+following method to the *InitializeSimulation* class:
.. label:: start_InitializeSimulation_class
.. code-block:: python
def define_box(self):
- """Define the simulation box.
- For 2D simulations, the third dimensions only contains 0.
- """
+ """Define the simulation box. Only 3D boxes are supported."""
box_boundaries = np.zeros((3, 2))
- dim = 0
- for L in self.box_dimensions:
+ for dim, L in enumerate(self.box_dimensions):
box_boundaries[dim] = -L/2, L/2
- dim += 1
self.box_boundaries = box_boundaries
box_size = np.diff(self.box_boundaries).reshape(3)
box_geometry = np.array([90, 90, 90])
@@ -156,12 +101,16 @@ method to the *InitializeSimulation* class:
The *box_boundaries* are calculated from the *box_dimensions*. They
represent the lowest and highest coordinates in all directions. By symmetry,
the box is centered at 0 along all axes. A *box_size* is also defined,
-following the MDAnalysis conventions: Lx, Ly, Lz, 90, 90, 90, where the
-last three numbers are angles in degrees. Values different from *90* for
-the angles would define a triclinic (non-orthogonal) box, which is not
-currently supported by the existing code.
+following the MDAnalysis |mdanalysis_box|: Lx, Ly, Lz, 90, 90, 90, where
+the last three numbers are angles in degrees :cite:`michaud2011mdanalysis`.
+Values different from *90* for the angles would define a triclinic (non-orthogonal)
+box, which is not currently supported by the existing code.
+
+.. |mdanalysis_box| raw:: html
-Let us call the *define_box* method from the *__init__* class:
+ conventions
+
+Let us call *define_box()* from the *__init__()* method:
.. label:: start_InitializeSimulation_class
@@ -169,18 +118,19 @@ Let us call the *define_box* method from the *__init__* class:
def __init__(self,
(...)
- self.nondimensionalize_units_1()
+ self.nondimensionalize_units(["box_dimensions", "cut_off",
+ "initial_positions"])
self.define_box()
.. label:: end_InitializeSimulation_class
-Populate the box
+Populate the Box
----------------
Here, the atoms are placed within the simulation box. If initial
positions were not provided (i.e., *initial_positions = None*), atoms
-are placed randomly within the box. If *initial_positions* was provided
-as an array, the provided positions are used instead. Note that, in this
+are placed randomly within the box. If *initial_positions* is provided
+as an array, the provided positions are used instead. Note that in this
case, the array must be of size 'number of atoms' times 'number of dimensions'.
.. label:: start_InitializeSimulation_class
@@ -188,11 +138,12 @@ case, the array must be of size 'number of atoms' times 'number of dimensions'.
.. code-block:: python
def populate_box(self):
+ Nat = np.sum(self.number_atoms) # total number of atoms
if self.initial_positions is None:
- atoms_positions = np.zeros((self.total_number_atoms, 3))
+ atoms_positions = np.zeros((Nat, 3))
for dim in np.arange(3):
diff_box = np.diff(self.box_boundaries[dim])
- random_pos = np.random.random(self.total_number_atoms)
+ random_pos = np.random.random(Nat)
atoms_positions[:, dim] = random_pos*diff_box-diff_box/2
self.atoms_positions = atoms_positions
else:
@@ -200,14 +151,14 @@ case, the array must be of size 'number of atoms' times 'number of dimensions'.
.. label:: end_InitializeSimulation_class
-In case *initial_positions is None*, and array is first created. Then, random
-positions that are constrained within the box boundaries are defined using the
-random function of NumPy. Note that, here, the newly added atoms are added
-randomly within the box, without taking care of avoiding overlaps with
-existing atoms. Overlaps will be dealt with using energy minimization,
-see :ref:`chapter4-label`.
+In case *initial_positions* is None, an array is first created. Then,
+random positions constrained within the box boundaries are defined using
+the random function from NumPy. Note that newly added atoms are placed
+randomly within the box, without considering overlaps with existing
+atoms. Overlaps will be addressed using energy minimization (see
+the :ref:`chapter4-label` chapter).
-Let us call the *populate_box* method from the *__init__* class:
+Let us call *populate_box* from the *__init__* method:
.. label:: start_InitializeSimulation_class
@@ -219,13 +170,114 @@ Let us call the *populate_box* method from the *__init__* class:
self.populate_box()
.. label:: end_InitializeSimulation_class
+
+Build Neighbor Lists
+--------------------
+
+In molecular simulations, it is common practice to identify neighboring atoms
+to save computational time. By focusing only on interactions between
+neighboring atoms, the simulation becomes more efficient. Add the following
+*update_neighbor_lists()* method to the *Utilities* class:
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+ def update_neighbor_lists(self):
+ if (self.step % self.neighbor == 0):
+ matrix = distances.contact_matrix(self.atoms_positions,
+ cutoff=self.cut_off, #+2,
+ returntype="numpy",
+ box=self.box_size)
+ neighbor_lists = []
+ for cpt, array in enumerate(matrix[:-1]):
+ list = np.where(array)[0].tolist()
+ list = [ele for ele in list if ele > cpt]
+ neighbor_lists.append(list)
+ self.neighbor_lists = neighbor_lists
+
+.. label:: end_Utilities_class
+
+The *update_neighbor_lists()* method generates neighbor lists for each
+atom, ensuring that only relevant interactions are considered in the
+calculations. These lists will be recalculated at intervals specified by
+the *neighbor* input parameter.
+
+For efficiency, the *contact_matrix* function from MDAnalysis is
+used :cite:`michaud2011mdanalysis`. The *contact_matrix* function returns
+information about atoms located at a distance less than the cutoff from one another.
+
+Update Cross Coefficients
+-------------------------
+
+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).
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+ def update_cross_coefficients(self):
+ if (self.step % self.neighbor == 0):
+ # Precalculte LJ cross-coefficients
+ sigma_ij_list = []
+ epsilon_ij_list = []
+ for Ni in np.arange(np.sum(self.number_atoms)-1): # tofix error for GCMC
+ # Read information about atom i
+ sigma_i = self.atoms_sigma[Ni]
+ epsilon_i = self.atoms_epsilon[Ni]
+ neighbor_of_i = self.neighbor_lists[Ni]
+ # Read information about neighbors j
+ sigma_j = self.atoms_sigma[neighbor_of_i]
+ epsilon_j = self.atoms_epsilon[neighbor_of_i]
+ # Calculare cross parameters
+ sigma_ij_list.append((sigma_i+sigma_j)/2)
+ epsilon_ij_list.append((epsilon_i+epsilon_j)/2)
+ self.sigma_ij_list = sigma_ij_list
+ self.epsilon_ij_list = epsilon_ij_list
+
+.. label:: end_Utilities_class
+
+Here, the values of the cross coefficients between atom of type 1 and 2,
+:math:`\sigma_{12}` and :math:`\epsilon_{12}`, are assumed to follow the arithmetic mean:
+
+.. math::
+
+ \sigma_{12} = (\sigma_{11}+\sigma_{22})/2 \\
+ \epsilon_{12} = (\epsilon_{11}+\epsilon_{22})/2
+
+Finally, import the following libraries in the *Utilities.py* file:
+
+.. label:: start_Utilities_class
+
+.. code-block:: python
+
+ import numpy as np
+ from MDAnalysis.analysis import distances
+
+.. label:: end_Utilities_class
+
+Let us call *update_neighbor_lists* and *update_cross_coefficients*
+from the *__init__* method:
+
+.. label:: start_InitializeSimulation_class
+
+.. code-block:: python
+
+ def __init__(self,
+ (...)
+ self.populate_box()
+ self.update_neighbor_lists()
+ self.update_cross_coefficients()
+
+.. label:: end_InitializeSimulation_class
-Test the code
+Test the Code
-------------
-Let us test the *InitializeSimulation* class to make sure that it does what
-is expected, i.e. that it does create atoms that are located within the box
-boundaries along all 3 dimensions of space:
+Let us test the *InitializeSimulation* class to ensure that it behaves as
+expected, i.e., that it creates atoms located within the box boundaries:
.. label:: start_test_3a_class
@@ -233,14 +285,31 @@ boundaries along all 3 dimensions of space:
import numpy as np
from InitializeSimulation import InitializeSimulation
-
- # Initialize the InitializeSimulation object
+ from pint import UnitRegistry
+ ureg = UnitRegistry()
+
+ # Define atom number of each group
+ nmb_1, nmb_2= [2, 3]
+ # Define LJ parameters (sigma)
+ sig_1, sig_2 = [3, 4]*ureg.angstrom
+ # Define LJ parameters (epsilon)
+ eps_1, eps_2 = [0.2, 0.4]*ureg.kcal/ureg.mol
+ # Define atom mass
+ mss_1, mss_2 = [10, 20]*ureg.gram/ureg.mol
+ # Define box size
+ L = 20*ureg.angstrom
+ # Define a cut off
+ rc = 2.5*sig_1
+
+ # Initialize the prepare object
init = InitializeSimulation(
- number_atoms=[2, 3],
- epsilon=[0.2, 0.4], # kcal/mol
- sigma=[3, 4], # A
- atom_mass=[10, 20], # g/mol
- box_dimensions=[20, 20, 20], # A
+ ureg = ureg,
+ number_atoms=[nmb_1, nmb_2],
+ 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,
)
# Test function using pytest
@@ -249,7 +318,8 @@ boundaries along all 3 dimensions of space:
atoms_positions = init.atoms_positions
for atom_position in atoms_positions:
for x, boundary in zip(atom_position, box_boundaries):
- assert (x >= boundary[0]) and (x <= boundary[1]), f"Test failed: Atoms outside of the box at position {atom_position}"
+ assert (x >= boundary[0]) and (x <= boundary[1]), \
+ f"Test failed: Atoms outside of the box at position {atom_position}"
print("Test passed")
# If the script is run directly, execute the tests
@@ -259,3 +329,6 @@ boundaries along all 3 dimensions of space:
pytest.main(["-s", __file__])
.. label:: end_test_3a_class
+
+The value of the cutoff, chosen as :math:`2.5 \sigma_{11}`, is relatively
+common for Lennard-Jones fluids and will often be the default choice for us.
diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst
index e460b16..e471092 100644
--- a/docs/source/chapters/chapter4.rst
+++ b/docs/source/chapters/chapter4.rst
@@ -32,8 +32,10 @@ The steepest descent method is used for energy minimization, following these ste
- 3) Move the atoms in the opposite direction of the maximum
force to minimize the potential energy by a displacement step.
The size of the step is adjusted iteratively based on the reduction in energy.
-- 4) Compute the new potential energy after the displacement, :math:`E_\text{pot}^\text{trial}`.
-- 5) Evaluate the change in energy: :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`.
+- 4) Compute the new potential energy after the displacement,
+ :math:`E_\text{pot}^\text{trial}`.
+- 5) Evaluate the change in energy:
+ :math:`\Delta E = E_\text{pot}^\text{trial} - E_\text{pot}^\text{initial}`.
- If :math:`\Delta E < 0`, the new configuration is accepted as it results in
lower energy, and the step size is increased.
@@ -47,19 +49,7 @@ minimized energy state.
Prepare the minimization
------------------------
-Let us start by importing NumPy and the copy libraries. Add the following
-to the beginning of the *MinimizeEnergy.py* file:
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
- import numpy as np
- import copy
-
-.. label:: end_MinimizeEnergy_class
-
-Then, let us fill the *__init__()* method:
+Let us fill the *__init__()* method:
.. label:: start_MinimizeEnergy_class
@@ -68,67 +58,15 @@ Then, let us fill the *__init__()* method:
class MinimizeEnergy(Measurements):
def __init__(self,
maximum_steps,
- cut_off=9,
- neighbor=1,
- displacement=0.01,
- thermo_outputs="MaxF",
- data_folder=None,
*args,
**kwargs):
- self.neighbor = neighbor
- self.cut_off = cut_off
- self.displacement = displacement
self.maximum_steps = maximum_steps
- self.thermo_outputs = thermo_outputs
- self.data_folder = data_folder
- if self.data_folder is not None:
- if os.path.exists(self.data_folder) is False:
- os.mkdir(self.data_folder)
super().__init__(*args, **kwargs)
-
-.. label:: end_MinimizeEnergy_class
-
-An important parameter is *maximum_steps*, which sets the maximum number
-of steps for the energy minimization process. A *cut_off* value with a
-default of 9 Ångströms is also defined. The *neighbor* parameter determines
-the interval between recalculations of the neighbor lists, and the *displacement*
-parameter, with a default value of 0.01 Ångström, sets the initial atom
-displacement value.
-
-The *thermo_outputs* and *data_folder* parameters are used for printing data
-to files. These two parameters will be useful in the next chapter, :ref:`chapter5-label`.
-
-Nondimensionalize units
------------------------
-
-As was done previously, some parameters from the *MinimizeEnergy* class
-must be non-dimensionalized: *cut_off* and *displacement*. Add the following
-method to the *MinimizeEnergy* class:
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
- def nondimensionalize_units_2(self):
- """Use LJ prefactors to convert units into non-dimensional."""
- self.cut_off = self.cut_off/self.reference_distance
- self.displacement = self.displacement/self.reference_distance
-
+
.. label:: end_MinimizeEnergy_class
-Let us call the *nondimensionalize_units_2()* method from the *__init__()*
-method:
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
- def __init__(self,
- (...)
- super().__init__(*args, **kwargs)
- self.nondimensionalize_units_2()
-
-.. label:: end_MinimizeEnergy_class
+The parameter *maximum_steps* sets the maximum number
+of steps for the energy minimization process.
Energy minimizer
----------------
@@ -141,14 +79,20 @@ following *run()* method to the *MinimizeEnergy* class:
.. code-block:: python
def run(self):
+ self.displacement = 0.01 # pick a random initial displacement (dimentionless)
# *step* loops for 0 to *maximum_steps*+1
for self.step in range(0, self.maximum_steps+1):
# First, meevaluate the initial energy and max force
self.update_neighbor_lists() # Rebuild neighbor list, if necessary
self.update_cross_coefficients() # Recalculate the cross coefficients, if necessary
# Compute Epot/MaxF/force
- init_Epot = self.compute_potential()
- forces, init_MaxF = self.compute_force()
+ if hasattr(self, 'Epot') is False: # If self.Epot does not exists yet, calculate it
+ self.Epot = self.compute_potential()
+ if hasattr(self, 'MaxF') is False: # If self.MaxF does not exists yet, calculate it
+ forces = self.compute_force()
+ self.MaxF = np.max(np.abs(forces))
+ init_Epot = self.Epot
+ init_MaxF = self.MaxF
# Save the current atom positions
init_positions = copy.deepcopy(self.atoms_positions)
# Move the atoms in the opposite direction of the maximum force
@@ -159,10 +103,9 @@ following *run()* method to the *MinimizeEnergy* class:
# Keep the more favorable energy
if trial_Epot < init_Epot: # accept new position
self.Epot = trial_Epot
- # calculate the new max force and save it
- forces, init_MaxF = self.compute_force()
+ forces = self.compute_force() # calculate the new max force and save it
self.MaxF = np.max(np.abs(forces))
- self.wrap_in_box() # Wrap atoms in the box, if necessary
+ self.wrap_in_box() # Wrap atoms in the box
self.displacement *= 1.2 # Multiply the displacement by a factor 1.2
else: # reject new position
self.Epot = init_Epot # Revert to old energy
@@ -171,100 +114,19 @@ following *run()* method to the *MinimizeEnergy* class:
.. label:: end_MinimizeEnergy_class
-The displacement, which has an initial value of 0.01, is adjusted through energy
-minimization. When the trial is successful, its value is multiplied by 1.2. When
-the trial is rejected, its value is multiplied by 0.2.
+The displacement, which has an initial value of 0.01, is adjusted through
+energy minimization. When the trial is successful, its value is multiplied
+by 1.2. When the trial is rejected, its value is multiplied by 0.2. Thus,
+as the minimization progresses, the displacements that the particles
+perform increase.
-Build neighbor lists
---------------------
+Compute the Potential
+---------------------
-In molecular simulations, it is common practice to identify neighboring atoms
-to save computational time. By focusing only on interactions between
-neighboring atoms, the simulation becomes more efficient. Add the following
-*update_neighbor_lists()* method to the *Utilities* class:
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
- def update_neighbor_lists(self):
- if (self.step % self.neighbor == 0):
- matrix = distances.contact_matrix(self.atoms_positions,
- cutoff=self.cut_off, #+2,
- returntype="numpy",
- box=self.box_size)
- neighbor_lists = []
- for cpt, array in enumerate(matrix[:-1]):
- list = np.where(array)[0].tolist()
- list = [ele for ele in list if ele > cpt]
- neighbor_lists.append(list)
- self.neighbor_lists = neighbor_lists
-
-.. label:: end_Utilities_class
-
-The *update_neighbor_lists()* method generates neighbor lists for each
-atom, ensuring that only relevant interactions are considered in the
-calculations. These lists will be recalculated at intervals specified by
-the *neighbor* input parameter.
-
-Update cross coefficients
--------------------------
-
-At the same time as the neighbor lists are getting build up, let us also
-pre-calculate the cross coefficients. This will make the force calculation
-more practical (see below).
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
- def update_cross_coefficients(self):
- if (self.step % self.neighbor == 0):
- # Precalculte LJ cross-coefficients
- sigma_ij_list = []
- epsilon_ij_list = []
- for Ni in np.arange(self.total_number_atoms-1): # tofix error for GCMC
- # Read information about atom i
- sigma_i = self.atoms_sigma[Ni]
- epsilon_i = self.atoms_epsilon[Ni]
- neighbor_of_i = self.neighbor_lists[Ni]
- # Read information about neighbors j
- sigma_j = self.atoms_sigma[neighbor_of_i]
- epsilon_j = self.atoms_epsilon[neighbor_of_i]
- # Calculare cross parameters
- sigma_ij_list.append((sigma_i+sigma_j)/2)
- epsilon_ij_list.append((epsilon_i+epsilon_j)/2)
- self.sigma_ij_list = sigma_ij_list
- self.epsilon_ij_list = epsilon_ij_list
-
-.. label:: end_Utilities_class
-
-Here, the values of the cross coefficients between atom of type 1 and 2,
-:math:`\sigma_{12}` and :math:`\epsilon_{12}`, are assumed to follow the arithmetic mean:
-
-.. math::
-
- \sigma_{12} = (\sigma_{11}+\sigma_{22})/2 \\
- \epsilon_{12} = (\epsilon_{11}+\epsilon_{22})/2
-
-Finally, import the following library in the *Utilities.py* file:
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
- import numpy as np
- from MDAnalysis.analysis import distances
-
-.. label:: end_Utilities_class
-
-Compute_potential
------------------
-
-Computing the potential energy of the system is central to the energy minimizer,
-as the value of the potential is used to decide if the trial is accepted or
-rejected. Add the following method called *compute_potential()* to the *Utilities*
-class:
+Computing the potential energy of the system is central to the energy
+minimizer, as the value of the potential is used to decide if the trial is
+accepted or rejected. Add the following method called *compute_potential()*
+to the *Utilities* class:
.. label:: start_Utilities_class
@@ -273,24 +135,23 @@ class:
def compute_potential(self):
"""Compute the potential energy by summing up all pair contributions."""
energy_potential = 0
- for Ni in np.arange(self.total_number_atoms-1):
+ for Ni in np.arange(np.sum(self.number_atoms)-1):
# Read neighbor list
neighbor_of_i = self.neighbor_lists[Ni]
# Measure distance
rij = self.compute_distance(self.atoms_positions[Ni],
self.atoms_positions[neighbor_of_i],
- self.box_size[:3])
- # Measure potential using information about cross coefficients
+ self.box_size)
+ # Measure potential using pre-calculated cross coefficients
sigma_ij = self.sigma_ij_list[Ni]
epsilon_ij = self.epsilon_ij_list[Ni]
- energy_potential += np.sum(potentials(self.potential_type,
- epsilon_ij, sigma_ij, rij))
+ energy_potential += np.sum(potentials(epsilon_ij, sigma_ij, rij))
return energy_potential
.. label:: end_Utilities_class
-Measuring the distance is an important step of computing the potential. Let us
-do it using a dedicated method. Add the following method to the *Utilities*
+Measuring the distance is a central step in computing the potential. Let us
+do this using a dedicated method. Add the following method to the *Utilities*
class as well:
.. label:: start_Utilities_class
@@ -300,11 +161,10 @@ class as well:
def compute_distance(self,position_i, positions_j, box_size, only_norm = True):
"""
Measure the distances between two particles.
- The nan_to_num is crutial in 2D to avoid nan value along third dimension.
- # TOFIX: Move as function instead of a method?
+ # TOFIX: Move as a function instead of a method?
"""
rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j
- + box_size[:3]/2.0, box_size) - box_size[:3]/2.0)
+ + box_size[:3]/2.0, box_size[:3]) - box_size[:3]/2.0)
if only_norm:
return np.linalg.norm(rij_xyz, axis=1)
else:
@@ -314,7 +174,9 @@ class as well:
Finally, the energy minimization requires the computation of the minimum
force in the system. Although not very different from the potential measurement,
-let us create a new method that is dedicated solely to measuring forces:
+let us create a new method that is dedicated solely to measuring forces.
+The force can be returned as a vector (default), or as a matrix (which will be
+useful for pressure calculation, see the :ref:`chapter7-label` chapter):
.. label:: start_Utilities_class
@@ -322,22 +184,21 @@ let us create a new method that is dedicated solely to measuring forces:
def compute_force(self, return_vector = True):
if return_vector: # return a N-size vector
- force_vector = np.zeros((self.total_number_atoms,3))
+ force_vector = np.zeros((np.sum(self.number_atoms),3))
else: # return a N x N matrix
- force_matrix = np.zeros((self.total_number_atoms,
- self.total_number_atoms,3))
- for Ni in np.arange(self.total_number_atoms-1):
+ force_matrix = np.zeros((np.sum(self.number_atoms),
+ np.sum(self.number_atoms),3))
+ for Ni in np.arange(np.sum(self.number_atoms)-1):
# Read neighbor list
neighbor_of_i = self.neighbor_lists[Ni]
# Measure distance
rij, rij_xyz = self.compute_distance(self.atoms_positions[Ni],
self.atoms_positions[neighbor_of_i],
- self.box_size[:3], only_norm = False)
+ self.box_size, only_norm = False)
# Measure force using information about cross coefficients
sigma_ij = self.sigma_ij_list[Ni]
epsilon_ij = self.epsilon_ij_list[Ni]
- fij_xyz = potentials(self.potential_type, epsilon_ij,
- sigma_ij, rij, derivative = True)
+ fij_xyz = potentials(epsilon_ij, sigma_ij, rij, derivative = True)
if return_vector:
# Add the contribution to both Ni and its neighbors
force_vector[Ni] += np.sum((fij_xyz*rij_xyz.T/rij).T, axis=0)
@@ -346,23 +207,20 @@ let us create a new method that is dedicated solely to measuring forces:
# Add the contribution to the matrix
force_matrix[Ni][neighbor_of_i] += (fij_xyz*rij_xyz.T/rij).T
if return_vector:
- max_force = np.max(np.abs(force_vector))
- return force_vector, max_force
+ return force_vector
else:
return force_matrix
.. label:: end_Utilities_class
-Here, two types of outputs can
-be requested by the user: *force-vector*, and *force-matrix*.
-The *force-matrix* option will be useful for pressure calculation, see
-:ref:`chapter7-label`.
+The force can be returned as a vector (default) or as a matrix. The latter
+will be useful for pressure calculation; see the :ref:`chapter7-label` chapter.
-Wrap in box
+Wrap in Box
-----------
-Every time atoms are being displaced, one has to ensure that they remain in
-the box. This is done by the *wrap_in_box()* method that must be placed
+Every time atoms are displaced, one has to ensure that they remain within
+the box. This is done by the *wrap_in_box()* method, which must be placed
within the *Utilities* class:
.. label:: start_Utilities_class
@@ -370,7 +228,7 @@ within the *Utilities* class:
.. code-block:: python
def wrap_in_box(self):
- for dim in np.arange(self.dimensions):
+ for dim in np.arange(3):
out_ids = self.atoms_positions[:, dim] \
> self.box_boundaries[dim][1]
self.atoms_positions[:, dim][out_ids] \
@@ -394,15 +252,32 @@ typically negative.
.. code-block:: python
from MinimizeEnergy import MinimizeEnergy
-
- # Initialize the MinimizeEnergy object and run the minimization
+ from pint import UnitRegistry
+ ureg = UnitRegistry()
+
+ # Define atom number of each group
+ nmb_1, nmb_2= [2, 3]
+ # Define LJ parameters (sigma)
+ sig_1, sig_2 = [3, 4]*ureg.angstrom
+ # Define LJ parameters (epsilon)
+ eps_1, eps_2 = [0.2, 0.4]*ureg.kcal/ureg.mol
+ # Define atom mass
+ mss_1, mss_2 = [10, 20]*ureg.gram/ureg.mol
+ # Define box size
+ L = 20*ureg.angstrom
+ # Define a cut off
+ rc = 2.5*sig_1
+
+ # Initialize the prepare object
minimizer = MinimizeEnergy(
+ ureg = ureg,
maximum_steps=100,
- number_atoms=[2, 3],
- epsilon=[0.2, 0.4], # kcal/mol
- sigma=[3, 4], # A
- atom_mass=[10, 20], # g/mol
- box_dimensions=[20, 20, 20], # A
+ number_atoms=[nmb_1, nmb_2],
+ 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,
)
minimizer.run()
@@ -423,5 +298,6 @@ typically negative.
.. label:: end_test_4a_class
-For such as low density in particle, we can reasonably expect the energy to be always
-negative after 100 steps.
\ No newline at end of file
+For such a low particle density, we can reasonably expect that the potential
+energy will always be negative after 100 steps, and that the maximum force
+will be smaller than 10 (unitless).
diff --git a/docs/source/chapters/chapter5.rst b/docs/source/chapters/chapter5.rst
index b15b715..3c04c34 100644
--- a/docs/source/chapters/chapter5.rst
+++ b/docs/source/chapters/chapter5.rst
@@ -4,39 +4,19 @@ Output the simulation state
===========================
Here, the code is modified to allow us to follow the evolution of the system
-during the simulation. To do so, the *Output* class is modified.
-
-Update the MinimizeEnergy class
--------------------------------
-
-Let us start by calling two additional methods within the for loop of the
-*MinimizeEnergy* class, within the *run()* method.
-
-.. label:: start_MinimizeEnergy_class
-
-.. code-block:: python
-
- def run(self):
- (...)
- for self.step in range(0, self.maximum_steps+1):
- (...)
- self.displacement *= 0.2 # Multiply the displacement by a factor 0.2
- log_simulation_data(self)
- update_dump_file(self, "dump.min.lammpstrj")
-
-.. label:: end_MinimizeEnergy_class
-
-The two methods named *update_log_minimize()* and *update_dump_file()*, are used
-to print the information in the terminal and in a LAMMPS-type data file, respectively.
-These two methods will be written in the following.
+during the simulation. To do so, two new files will be created: one dedicated
+to logging information, and the other to printing the atom positions to a
+file, thus allowing for their visualization with software like
+VMD :cite:`humphrey1996vmd` or Ovito :cite:`stukowski2009visualization`.
Create logger
-------------
-Let us create functions named *log_simulation_data* to a file named *logger.py*.
+Let us create a function named *log_simulation_data()* to a file named *logger.py*.
With the logger, some output are being printed in a file, as well as in the terminal.
-The frequency of printing is set by *thermo_period*, see :ref:`chapter3-label`.
-All quantities are re-dimensionalized before getting outputed.
+The period of printing, which is a multiple of the simulation steps, will be set
+by a new parameter named *thermo_period* (integer). All quantities are
+converted into *real* units before getting outputed (see :ref:`chapter2-label`).
.. label:: start_logger_class
@@ -76,6 +56,7 @@ All quantities are re-dimensionalized before getting outputed.
return logger
def log_simulation_data(code):
+ # TOFIX: ceurrently, MaxF is returned dimensionless
# Setup the logger with the folder name, overwriting the log if code.step is 0
logger = setup_logger(code.data_folder, overwrite=(code.step == 0))
@@ -85,9 +66,9 @@ All quantities are re-dimensionalized before getting outputed.
if code.step % code.thermo_period == 0:
if code.step == 0:
Epot = code.compute_potential() \
- * code.reference_energy # kcal/mol
+ * code.ref_energy # kcal/mol
else:
- Epot = code.Epot * code.reference_energy # kcal/mol
+ Epot = code.Epot * code.ref_energy # kcal/mol
if code.step == 0:
if code.thermo_outputs == "Epot":
logger.info(f"step Epot")
@@ -96,24 +77,32 @@ All quantities are re-dimensionalized before getting outputed.
elif code.thermo_outputs == "Epot-press":
logger.info(f"step Epot press")
if code.thermo_outputs == "Epot":
- logger.info(f"{code.step} {Epot:.2f}")
+ logger.info(f"{code.step} {Epot.magnitude:.2f}")
elif code.thermo_outputs == "Epot-MaxF":
- logger.info(f"{code.step} {Epot:.2f} {code.MaxF:.2f}")
+ logger.info(f"{code.step} {Epot.magnitude:.2f} {code.MaxF:.2f}")
elif code.thermo_outputs == "Epot-press":
code.calculate_pressure()
- press = code.pressure * code.reference_pressure # Atm
- logger.info(f"{code.step} {Epot:.2f} {press:.2f}")
+ press = code.pressure * code.ref_pressure # Atm
+ logger.info(f"{code.step} {Epot.magnitude:.2f} {press.magnitude:.2f}")
.. label:: end_logger_class
-Create dumper
+For simplicify, the *logger* uses the |logging| library, which provides a
+flexible framework for emitting log messages from Python programs. Depending on
+the value of *thermo_outputs*, different informations are printed, such as
+*step*, *Epot*, or/and *MaxF*.
+
+.. |logging| raw:: html
+
+ logging
+
+Create Dumper
-------------
-Let us create a function named *update_dump_file* to a file named
-*dumper.py*. The dumper will print a *.lammpstrj file*, which contains the box
-dimensions and atom positions at every chosen frame (set by *dumping_period*,
-see :ref:`chapter3-label`). All quantities are dimensionalized before getting outputed, and the file follows
-a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
+Let us create a function named *update_dump_file()* in a file named
+*dumper.py*. The dumper will print a *.lammpstrj* file, which contains
+the box dimensions and atom positions at every chosen frame (set by
+*dumping_period*). All quantities are dimensionalized before being output.
.. label:: start_dumper_class
@@ -125,14 +114,12 @@ a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
if code.dumping_period is not None:
if code.step % code.dumping_period == 0:
# Convert units to the *real* dimensions
- box_boundaries = code.box_boundaries\
- * code.reference_distance # Angstrom
- atoms_positions = code.atoms_positions\
- * code.reference_distance # Angstrom
+ box_boundaries = code.box_boundaries*code.ref_length # Angstrom
+ atoms_positions = code.atoms_positions*code.ref_length # Angstrom
atoms_types = code.atoms_type
if velocity:
atoms_velocities = code.atoms_velocities \
- * code.reference_distance/code.reference_time # Angstrom/femtosecond
+ * code.ref_length/code.ref_time # Angstrom/femtosecond
# Start writting the file
if code.step == 0: # Create new file
f = open(code.data_folder + filename, "w")
@@ -141,18 +128,18 @@ a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
f.write("ITEM: TIMESTEP\n")
f.write(str(code.step) + "\n")
f.write("ITEM: NUMBER OF ATOMS\n")
- f.write(str(code.total_number_atoms) + "\n")
+ f.write(str(np.sum(code.number_atoms)) + "\n")
f.write("ITEM: BOX BOUNDS pp pp pp\n")
for dim in np.arange(3):
- f.write(str(box_boundaries[dim][0]) + " "
- + str(box_boundaries[dim][1]) + "\n")
+ f.write(str(box_boundaries[dim][0].magnitude) + " "
+ + str(box_boundaries[dim][1].magnitude) + "\n")
cpt = 1
if velocity:
f.write("ITEM: ATOMS id type x y z vx vy vz\n")
characters = "%d %d %.3f %.3f %.3f %.3f %.3f %.3f %s"
for type, xyz, vxyz in zip(atoms_types,
- atoms_positions,
- atoms_velocities):
+ atoms_positions.magnitude,
+ atoms_velocities.magnitude):
v = [cpt, type, xyz[0], xyz[1], xyz[2],
vxyz[0], vxyz[1], vxyz[2]]
f.write(characters % (v[0], v[1], v[2], v[3], v[4],
@@ -162,7 +149,7 @@ a LAMMPS dump format, and can be read by molecular dynamics softwares like VMD.
f.write("ITEM: ATOMS id type x y z\n")
characters = "%d %d %.3f %.3f %.3f %s"
for type, xyz in zip(atoms_types,
- atoms_positions):
+ atoms_positions.magnitude):
v = [cpt, type, xyz[0], xyz[1], xyz[2]]
f.write(characters % (v[0], v[1], v[2],
v[3], v[4], '\n'))
@@ -197,39 +184,115 @@ Add the same lines at the top of the *MinimizeEnergy.py* file:
.. label:: end_MinimizeEnergy_class
+Finally, let us make sure that *thermo_period*, *dumping_period*, and *thermo_outputs*
+parameters are passed the InitializeSimulation method:
+
+.. label:: start_InitializeSimulation_class
+
+.. code-block:: python
+
+ def __init__(self,
+ (...)
+ neighbor=1, # Integer
+ thermo_period = None,
+ dumping_period = None,
+ thermo_outputs = None,
+ data_folder="Outputs/",
+
+.. label:: end_InitializeSimulation_class
+
+.. label:: start_InitializeSimulation_class
+
+.. code-block:: python
+
+ def __init__(self,
+ (...)
+ self.initial_positions = initial_positions
+ self.thermo_period = thermo_period
+ self.dumping_period = dumping_period
+ self.thermo_outputs = thermo_outputs
+ self.data_folder = data_folder
+ if os.path.exists(self.data_folder) is False:
+ os.mkdir(self.data_folder)
+
+.. label:: end_InitializeSimulation_class
+
+Update the MinimizeEnergy class
+-------------------------------
+
+As a final step, let us call two functions *log_simulation_data* and
+*update_dump_file* within the *for* loop of the
+*MinimizeEnergy* class, within the *run()* method:
+
+.. label:: start_MinimizeEnergy_class
+
+.. code-block:: python
+
+ def run(self):
+ (...)
+ for self.step in range(0, self.maximum_steps+1):
+ (...)
+ self.displacement *= 0.2 # Multiply the displacement by a factor 0.2
+ log_simulation_data(self)
+ update_dump_file(self, "dump.min.lammpstrj")
+
+.. label:: end_MinimizeEnergy_class
+
+The name *dump.min.lammpstrj* will be used when printing the dump file
+during energy minimization.
+
Test the code
-------------
-One can use a test similar as the previous ones. Let us ask out code to print
-information in the dump and the log files, and then let us assert the
+One can use a test similar as the previous ones. Let us ask our code to print
+information in the *dump* and the *log* files, and then let us assert that the
files were indeed created without the *Outputs/* folder:
.. label:: start_test_5a_class
.. code-block:: python
- import os
from MinimizeEnergy import MinimizeEnergy
+ from pint import UnitRegistry
+ ureg = UnitRegistry()
+ import os
- # Initialize the MinimizeEnergy object and run the minimization
+ # Define atom number of each group
+ nmb_1, nmb_2= [2, 3]
+ # Define LJ parameters (sigma)
+ sig_1, sig_2 = [3, 4]*ureg.angstrom
+ # Define LJ parameters (epsilon)
+ eps_1, eps_2 = [0.2, 0.4]*ureg.kcal/ureg.mol
+ # Define atom mass
+ mss_1, mss_2 = [10, 20]*ureg.gram/ureg.mol
+ # Define box size
+ L = 20*ureg.angstrom
+ # Define a cut off
+ rc = 2.5*sig_1
+
+ # Initialize the prepare object
minimizer = MinimizeEnergy(
+ ureg = ureg,
maximum_steps=100,
thermo_period=25,
dumping_period=25,
- thermo_outputs="Epot-MaxF",
- number_atoms=[2, 3],
- epsilon=[0.1, 1.0], # kcal/mol
- sigma=[3, 4], # A
- atom_mass=[10, 20], # g/mol
- box_dimensions=[20, 20, 20], # A
+ number_atoms=[nmb_1, nmb_2],
+ 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,
data_folder="Outputs/",
+ thermo_outputs="Epot-MaxF",
)
minimizer.run()
# Test function using pytest
def test_output_files():
- assert os.path.exists("Outputs/dump.min.lammpstrj"), "Test failed: dump file was not created"
- assert os.path.exists("Outputs/simulation.log"), "Test failed: log file was not created"
+ assert os.path.exists("Outputs/dump.min.lammpstrj"), \
+ "Test failed: the dump file was not created"
+ assert os.path.exists("Outputs/simulation.log"), \
+ "Test failed: the log file was not created"
print("Test passed")
# If the script is run directly, execute the tests
@@ -240,8 +303,9 @@ files were indeed created without the *Outputs/* folder:
.. label:: end_test_5a_class
-I addition to the files getting created, information must be printed in the terminal
-during the simulation:
+In addition to making sure that both files were created, let us verify
+that the expected information was printed to the terminal during the
+simulation. The content of the *simulation.log* file should resemble:
.. code-block:: bw
@@ -254,8 +318,8 @@ during the simulation:
The data from the *simulation.log* can be used to generate plots using softwares
line XmGrace, GnuPlot, or Python/Pyplot. For the later, one can use a simple data
-reader to import the data from *Outputs/simulation.log* into Python. Copy the
-following lines in a file named *reader.py*:
+reader to import the data from *simulation.log* into Python. Copy the
+following lines in a new file named *reader.py*:
.. label:: start_reader_class
@@ -294,14 +358,28 @@ The *import_data* function from *reader.py* can simply be used as follows:
.. label:: start_test_5b_class
+.. code-block:: python
+
from reader import import_data
- file_path = "Outputs/simulation.log"
- header, data = import_data(file_path)
+ def test_file_not_empty():
+ # Import data from the file
+ file_path = "Outputs/simulation.log"
+ header, data = import_data(file_path)
+ # Check if the header and data are not empty
+ assert header, "Error, no header in simulation.log"
+ assert data, "Error, no data in simulation.log"
+ assert len(data) > 0, "Data should contain at least one row"
- print(header)
- for row in data:
- print(row)
+ print(header)
+ for row in data:
+ print(row)
+
+ # If the script is run directly, execute the tests
+ if __name__ == "__main__":
+ import pytest
+ # Run pytest programmatically
+ pytest.main(["-s", __file__])
.. label:: end_test_5b_class
@@ -314,4 +392,4 @@ Which must return:
[25.0, -2.12, 1.22]
[50.0, -2.19, 2.85]
[75.0, -2.64, 0.99]
- [100.0, -2.64, 0.51]
\ No newline at end of file
+ [100.0, -2.64, 0.51]
diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst
index 52fe593..6b7ea4d 100644
--- a/docs/source/chapters/chapter6.rst
+++ b/docs/source/chapters/chapter6.rst
@@ -47,22 +47,23 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
def monte_carlo_move(self):
"""Monte Carlo move trial."""
if self.displace_mc is not None: # only trigger if displace_mc was provided by the user
- try: # try using the last saved Epot, if it exists
- initial_Epot = self.Epot
- except: # If self.Epot does not exists yet, calculate it
- initial_Epot = self.compute_potential()
- # Make a copy of the initial atoms positions
+ # When needed, recalculate neighbor/coeff lists
+ self.update_neighbor_lists()
+ self.update_cross_coefficients()
+ # If self.Epot does not exists yet, calculate it
+ # It should only be necessary when step = 0
+ if hasattr(self, 'Epot') is False:
+ self.Epot = self.compute_potential()
+ # Make a copy of the initial atoms positions and initial energy
+ initial_Epot = self.Epot
initial_positions = copy.deepcopy(self.atoms_positions)
# Pick an atom id randomly
- atom_id = np.random.randint(self.total_number_atoms)
+ atom_id = np.random.randint(np.sum(self.number_atoms))
# Move the chosen atom in a random direction
# The maximum displacement is set by self.displace_mc
- if self.dimensions == 3:
- move = (np.random.random(3)-0.5)*self.displace_mc
- elif self.dimensions == 2: # the third value will be 0
- move = np.append((np.random.random(2) - 0.5) * self.displace_mc, 0)
+ move = (np.random.random(3)-0.5)*self.displace_mc
self.atoms_positions[atom_id] += move
- # Measure the optential energy of the new configuration
+ # 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
@@ -71,19 +72,23 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class:
acceptation_probability = np.min([1, np.exp(-beta*delta_E)])
if random_number <= acceptation_probability: # Accept new position
self.Epot = trial_Epot
+ self.successful_move += 1
else: # Reject new position
- self.Epot = initial_Epot
self.atoms_positions = initial_positions # Revert to initial positions
+ self.failed_move += 1
.. label:: end_MonteCarlo_class
+The counters *successful_move* and *failed_move* are incremented with each
+successful and failed attempt, respectively.
+
Parameters
----------
-The *monte_carlo_move* method requires a few parameters to be selected by the
-users, such as *displace_mc* (:math:`d_\text{mc}`), the maximum number of steps,
-and the desired temperature (:math:`T`). Let us add these parameters to the
-*__init__* method:
+The *monte_carlo_move* method requires a few parameters to be selected by
+the users, such as *displace_mc* (:math:`d_\text{mc}`, in Ångströms), the
+maximum number of steps, and the desired temperature (:math:`T`, in Kelvin).
+Let us add these parameters to the *__init__()* method:
.. label:: start_MonteCarlo_class
@@ -92,52 +97,25 @@ and the desired temperature (:math:`T`). Let us add these parameters to the
class MonteCarlo(Measurements):
def __init__(self,
maximum_steps,
- cut_off = 9,
+ desired_temperature,
displace_mc = None,
- neighbor = 1,
- desired_temperature = 300,
- thermo_outputs = "press",
- data_folder = None,
*args,
**kwargs):
self.maximum_steps = maximum_steps
- self.cut_off = cut_off
self.displace_mc = displace_mc
- self.neighbor = neighbor
self.desired_temperature = desired_temperature
- self.thermo_outputs = thermo_outputs
- self.data_folder = data_folder
- if self.data_folder is not None:
- if os.path.exists(self.data_folder) is False:
- os.mkdir(self.data_folder)
super().__init__(*args, **kwargs)
- self.nondimensionalize_units_3()
-
-.. label:: end_MonteCarlo_class
-
-Here, we anticipate that some of the parameters have to be nondimensionalized, which
-is done with the *nondimensionalize_units_3* method that must also be added to
-the *MonteCarlo* class:
-
-.. label:: start_MonteCarlo_class
-
-.. code-block:: python
-
- def nondimensionalize_units_3(self):
- """Use LJ prefactors to convert units into non-dimensional."""
- self.cut_off = self.cut_off/self.reference_distance
- self.desired_temperature = self.desired_temperature \
- /self.reference_temperature
- if self.displace_mc is not None:
- self.displace_mc /= self.reference_distance
+ self.nondimensionalize_units(["desired_temperature", "displace_mc"])
+ self.successful_move = 0
+ self.failed_move = 0
.. label:: end_MonteCarlo_class
-Run method
+Run Method
----------
-Finally, let us add a *run* method to the *MonteCarlo* class, that is used to
-perform a loop over the desired number of steps *maximum_steps*:
+Finally, let us add a *run* method to the *MonteCarlo* class that performs
+a loop over the desired number of steps, *maximum_steps*:
.. label:: start_MonteCarlo_class
@@ -146,20 +124,15 @@ perform a loop over the desired number of steps *maximum_steps*:
def run(self):
"""Perform the loop over time."""
for self.step in range(0, self.maximum_steps+1):
- self.update_neighbor_lists()
- self.update_cross_coefficients()
self.monte_carlo_move()
self.wrap_in_box()
.. label:: end_MonteCarlo_class
-At each step, the *monte_carlo_move* method is called. The previously defined
-methods *update_neighbor_lists* and *wrap_in_box* are also called to ensure that
-the neighbor lists are kept up to date despite the motion of the atoms, and that
-the atoms remain inside the box, respectively.
-
-Let us call *update_log_md_mc* from the run method of the MonteCarlo class.
-Let us add a dump too:
+At each step, the *monte_carlo_move()* method is called. The previously
+defined *wrap_in_box* method is also called to ensure that the atoms remain
+inside the box. Additionally, let us call *log_simulation_data()* and
+*update_dump_file()*:
.. label:: start_MonteCarlo_class
@@ -175,50 +148,73 @@ Let us add a dump too:
.. label:: end_MonteCarlo_class
-To output the density, let us add the following method to the *Utilities* class:
-# TOFIX: not used yet
-
-.. label:: start_Utilities_class
-
-.. code-block:: python
-
- def calculate_density(self):
- """Calculate the mass density."""
- volume = np.prod(self.box_size[:3]) # Unitless
- total_mass = np.sum(self.atoms_mass) # Unitless
- return total_mass/volume # Unitless
-
-.. label:: end_Utilities_class
-
-Test the code
+Test the Code
-------------
-One can use a similar test as previously. Let us use a displace distance of
-0.5 Angstrom, and make 1000 steps.
+Let us use a similar test as before. Set a displacement distance corresponding
+to a quarter of sigma, and perform a very small number of steps:
.. label:: start_test_6a_class
.. code-block:: python
- import os
from MonteCarlo import MonteCarlo
+ from pint import UnitRegistry
+ ureg = UnitRegistry()
+ import os
- mc = MonteCarlo(maximum_steps=1000,
- dumping_period=100,
- thermo_period=100,
- thermo_outputs = "Epot",
- displace_mc = 0.5,
- number_atoms=[50],
- epsilon=[0.1], # kcal/mol
- sigma=[3], # A
- atom_mass=[10], # g/mol
- box_dimensions=[20, 20, 20], # A
- data_folder="Outputs/",
- )
+ # Define atom number of each group
+ nmb_1= 50
+ # Define LJ parameters (sigma)
+ sig_1 = 3*ureg.angstrom
+ # Define LJ parameters (epsilon)
+ eps_1 = 0.1*ureg.kcal/ureg.mol
+ # Define atom mass
+ mss_1 = 10*ureg.gram/ureg.mol
+ # Define box size
+ L = 20*ureg.angstrom
+ # Define a cut off
+ rc = 2.5*sig_1
+ # Pick the desired temperature
+ T = 300*ureg.kelvin
+ # choose the displace_mc
+ displace_mc = sig_1/4
+
+ # Initialize the prepare object
+ mc = MonteCarlo(
+ 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
+ cut_off=rc,
+ thermo_outputs="Epot",
+ desired_temperature=T, # K
+ neighbor=20,
+ displace_mc = displace_mc,
+ )
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")
+
+ # If the script is run directly, execute the tests
+ if __name__ == "__main__":
+ import pytest
+ # Run pytest programmatically
+ pytest.main(["-s", __file__])
+
.. label:: end_test_6a_class
The evolution of the potential energy as a function of the number of steps
-are written in the *simulation.log* file. The data can be used to plot
-the evolution of the system with time.
+is recorded in the *simulation.log* file. The data in *simulation.log* can
+be used to plot the evolution of the system over time.
diff --git a/docs/source/chapters/chapter7.rst b/docs/source/chapters/chapter7.rst
index eb4d26c..0d0a8e8 100644
--- a/docs/source/chapters/chapter7.rst
+++ b/docs/source/chapters/chapter7.rst
@@ -1,38 +1,40 @@
.. _chapter7-label:
-Pressure measurement
+Pressure Measurement
====================
-In order to extract the equation of state in our simulation, we need to measure the
-pressure of the system, :math:`p`. The pressure in a molecular simulation can be
-calculated from the interactions between particles. The pressure can be measured as
-the sum of the ideal contribution, :math:`p_\text{ideal} = N_\text{DOF} k_\text{B} T / V d`,
-which comes from the ideal gas law, and a Virial term which accounts for the
-pressure contribution from the forces between particles,
+To extract the equation of state in our simulation, we need to measure the
+pressure of the system, :math:`p`. The pressure in a molecular simulation can
+be calculated from the interactions between particles. The pressure can be
+measured as the sum of the ideal contribution,
+:math:`p_\text{ideal} = N_\text{DOF} k_\text{B} T / V d`,
+which comes from the ideal gas law, and a Virial term that
+accounts for the pressure contribution from the forces between particles,
:math:`p_\text{non_ideal} = \left< \sum_i r_i \cdot F_i \right> / V d`. The final
-expression reads:
+expression reads :cite:`frenkel2023understanding`:
-.. math::
+.. math::
p = \dfrac{1}{V d} \left[ N_\text{DOF} k_\text{B} T + \left< \sum_i r_i \cdot F_i \right> \right]
:math:`N_\text{DOF}` is the number of degrees-of-freedom, which can be calculated
-from the number of particles, :math:`N`, and the dimension of the system, :math:`d`, as
-:math:`N_\text{DOF} = d N - d` :cite:`frenkel2023understanding`.
+from the number of particles, :math:`N`, and the dimension of the system,
+:math:`d = 3`, as :math:`N_\text{DOF} = d N - d` :cite:`frenkel2023understanding`.
-The calculation of :math:`p_\text{ideal}` is straighforward. For Monte Carlo simulation,
-as atoms do not have temperature, the *imposed* temperature will be used instead.
-The calculation of :math:`p_\text{non_ideal}` requires the measurement of all the
-force and distance between the atoms. The calculation of the forces was already
-implemented in a previous chapter, but a new function that returns all the
-vector direction between atoms pairs will have to be written here.
+The calculation of :math:`p_\text{ideal}` is straightforward. For Monte Carlo
+simulation, as atoms do not have a temperature, the *imposed* temperature will
+be used as ":math:`T`". The calculation of :math:`p_\text{non_ideal}` requires the
+measurement of all the forces and distances between atoms. The calculation of
+the forces was already implemented in a previous chapter (see *compute_force*
+the :ref:`chapter4-label` chapter), but a new function that returns all the
+vector directions between atom pairs will need to be written here.
Implement the Virial equation
-----------------------------
-Let us add the following method to the *Utilities* class.
+Let us add the following method to the *Measurements* class:
-.. label:: start_Utilities_class
+.. label:: start_Measurements_class
.. code-block:: python
@@ -40,21 +42,24 @@ Let us add the following method to the *Utilities* class.
"""Evaluate p based on the Virial equation (Eq. 4.4.2 in Frenkel-Smit,
Understanding molecular simulation: from algorithms to applications, 2002)"""
# Compute the ideal contribution
- Ndof = self.dimensions*self.total_number_atoms-self.dimensions
- volume = np.prod(self.box_size[:self.dimensions])
+ Nat = np.sum(self.number_atoms) # total number of atoms
+ dimension = 3 # 3D is the only possibility here
+ Ndof = dimension*Nat-dimension
+ volume = np.prod(self.box_size[:3]) # box volume
try:
self.calculate_temperature() # this is for later on, when velocities are computed
temperature = self.temperature
except:
temperature = self.desired_temperature # for MC, simply use the desired temperature
- p_ideal = Ndof*temperature/(volume*self.dimensions)
+ p_ideal = Ndof*temperature/(volume*dimension)
# Compute the non-ideal contribution
- distances_forces = np.sum(self.compute_force(return_vector = False)*self.evaluate_rij_matrix())
- p_nonideal = distances_forces/(volume*self.dimensions)
+ distances_forces = np.sum(self.compute_force(return_vector = False) \
+ *self.evaluate_rij_matrix())
+ p_nonideal = distances_forces/(volume*dimension)
# Final pressure
self.pressure = p_ideal+p_nonideal
-.. label:: end_Utilities_class
+.. label:: end_Measurements_class
To evaluate all the vectors between all the particles, let us also add the
*evaluate_rij_matrix* method to the *Utilities* class:
@@ -65,95 +70,77 @@ To evaluate all the vectors between all the particles, let us also add the
def evaluate_rij_matrix(self):
"""Matrix of vectors between all particles."""
- box_size = self.box_size[:3]
- half_box_size = self.box_size[:3]/2.0
- rij_matrix = np.zeros((self.total_number_atoms,self.total_number_atoms,3))
- positions_j = self.atoms_positions
- for Ni in range(self.total_number_atoms-1):
- position_i = self.atoms_positions[Ni]
- rij_xyz = (np.remainder(position_i - positions_j + half_box_size, box_size) - half_box_size)
+ Nat = np.sum(self.number_atoms)
+ Box = self.box_size[:3]
+ rij_matrix = np.zeros((Nat, Nat,3))
+ pos_j = self.atoms_positions
+ for Ni in range(Nat-1):
+ pos_i = self.atoms_positions[Ni]
+ rij_xyz = (np.remainder(pos_i - pos_j + Box/2.0, Box) - Box/2.0)
rij_matrix[Ni] = rij_xyz
- # use nan_to_num to avoid "nan" values in 2D
- return np.nan_to_num(rij_matrix)
+ return rij_matrix
.. label:: end_Utilities_class
Test the code
-------------
-Let us test the outputed pressure. An interesting test is to contront the output from
-our code with some data from the litterature. Let us used the same parameters
-as in Ref. :cite:`woodMonteCarloEquation1957`, where Monte Carlo simulations
-are used to simulate argon bulk phase. All we have to do is to apply our current
-code using their parameters, i.e. :math:`\sigma = 3.405~\text{Å}`, :math:`\epsilon = 0.238~\text{kcal/mol}`,
-and :math:`T = 55~^\circ\text{C}`. More details are given in the first illustration, :ref:`project1-label`.
-
-On the side note, a relatively small cut-off as well as a small number of atoms were
-chosen to make the calculation faster.
+Let us test the outputed pressure.
.. label:: start_test_7a_class
.. code-block:: python
- import numpy as np
from MonteCarlo import MonteCarlo
- from scipy import constants as cst
from pint import UnitRegistry
- from reader import import_data
-
- # Initialize the unit registry
ureg = UnitRegistry()
-
- # Constants
- kB = cst.Boltzmann * ureg.J / ureg.kelvin # Boltzmann constant
- Na = cst.Avogadro / ureg.mole # Avogadro's number
- R = kB * Na # Gas constant
-
- # Parameters taken from Wood1957
- tau = 2 # Ratio between volume / reduced volume
- epsilon = (119.76 * ureg.kelvin * kB * Na).to(ureg.kcal / ureg.mol) # kcal/mol
- r_star = 3.822 * ureg.angstrom # Angstrom
- sigma = r_star / 2**(1/6) # Angstrom
- N_atom = 50 # Number of atoms
- m_argon = 39.948 * ureg.gram / ureg.mol # g/mol
- T = 328.15 * ureg.degK # 328 K or 55°C
- volume_star = r_star**3 * Na * 2**(-0.5)
- cut_off = sigma * 2.5 # Angstrom
- displace_mc = sigma / 5 # Angstrom
- volume = N_atom * volume_star * tau / Na # Angstrom^3
- box_size = volume**(1/3) # Angstrom
-
- # Initialize and run the Monte Carlo simulation
+ import os
+
+ # Define atom number of each group
+ nmb_1= 50
+ # Define LJ parameters (sigma)
+ sig_1 = 3*ureg.angstrom
+ # Define LJ parameters (epsilon)
+ eps_1 = 0.1*ureg.kcal/ureg.mol
+ # Define atom mass
+ mss_1 = 10*ureg.gram/ureg.mol
+ # Define box size
+ L = 20*ureg.angstrom
+ # Define a cut off
+ rc = 2.5*sig_1
+ # Pick the desired temperature
+ T = 300*ureg.kelvin
+ # choose the displace_mc
+ displace_mc = sig_1/4
+
+ # Initialize the prepare object
mc = MonteCarlo(
- maximum_steps=15000,
- dumping_period=1000,
- thermo_period=1000,
+ 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
+ cut_off=rc,
thermo_outputs="Epot-press",
- neighbor=50,
- number_atoms=[N_atom],
- epsilon=[epsilon.magnitude],
- sigma=[sigma.magnitude],
- atom_mass=[m_argon.magnitude],
- box_dimensions=[box_size.magnitude, box_size.magnitude, box_size.magnitude],
- displace_mc=displace_mc.magnitude,
- desired_temperature=T.magnitude,
- cut_off=cut_off.magnitude,
- data_folder="Outputs/",
+ desired_temperature=T, # K
+ neighbor=20,
+ displace_mc = displace_mc,
)
+
+ # Run the Monte Carlo simulation
mc.run()
# Test function using pytest
- def test_pV_over_RT():
- # Import the data and calculate pV / RT
- pressure_data = np.array(import_data("Outputs/simulation.log")[1:])[0][:,2][10:] # Skip initial values
- pressure_atm = np.mean(pressure_data) # atm
- pressure_Pa = (pressure_atm * ureg.atm).to(ureg.pascal) # Pa
- volume = (volume_star * tau / Na).to(ureg.meter**3) # m3
- pV_over_RT = (pressure_Pa * volume / (R * T) * Na).magnitude
-
- # Assert that pV_over_RT is close to 1.5
- assert np.isclose(pV_over_RT, 1.5, atol=1.0), f"Test failed: pV/Rt = {pV_over_RT}, expected close to 1.5"
- print(f"pV/RT = {pV_over_RT:.2f} --- (The expected value from Wood1957 is 1.5)")
+ 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")
# If the script is run directly, execute the tests
if __name__ == "__main__":
@@ -163,12 +150,12 @@ chosen to make the calculation faster.
.. label:: end_test_7a_class
-Which should return a value for :math:`p V / R T` that is close to the expected value
-of 1.5 by Wood and Parker for :math:`\tau = V/V^* = 2` (see Fig. 4 in Ref. :cite:`woodMonteCarloEquation1957`):
+The pressure should be returned alongside the potential energy within *simulation.log*:
.. code-block:: bw
- (...)
- p v / R T = 1.56 --- (The expected value from Wood1957 is 1.5)
-
-The exact value will vary from one simulation to the other due to noise.
\ No newline at end of file
+ step Epot press
+ 0 134248.72 4608379.27
+ 10 124905.76 4287242.75
+ 20 124858.91 4285584.40
+ (...)
\ No newline at end of file
diff --git a/docs/source/chapters/chapter8.rst b/docs/source/chapters/chapter8.rst
index ba2926b..920c899 100644
--- a/docs/source/chapters/chapter8.rst
+++ b/docs/source/chapters/chapter8.rst
@@ -3,64 +3,138 @@
Monte Carlo insert
==================
-In the Grand Canonical ensemble, the number of atom in the
-simulation box is not constant.
-
-Let us add the following method to the MonteCarlo class:
+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 resemble the Monte Carlo move from :ref:`chapter6-label`:
+
+- 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 tried to be inserted, or an existing particle will tried to be deleted.
+- 3) The energy of the system after the insersion/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_insert_delete(self):
- if self.desired_mu is not None: # trigger method if desired_mu was entered
- initial_Epot = self.compute_potential(output="potential")
- initial_positions = copy.deepcopy(self.atoms_positions)
- initial_number_atoms = copy.deepcopy(self.number_atoms)
- initial_neighbor_lists = copy.deepcopy(self.neighbor_lists)
- volume = np.prod(np.diff(self.box_boundaries))
- if np.random.random() < 0.5:
- # Try adding an atom
- self.number_atoms[self.inserted_type] += 1
- atom_position = np.zeros((1, self.dimensions))
- for dim in np.arange(self.dimensions):
- atom_position[:, dim] = np.random.random(1)*np.diff(self.box_boundaries[dim]) - np.diff(self.box_boundaries[dim])/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], atom_position, self.atoms_positions[shift_id:]])
- self.update_neighbor_lists()
- self.calculate_cross_coefficients()
- trial_Epot = self.compute_potential(output="potential")
- Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type])
- beta = 1/self.desired_temperature
- acceptation_probability = np.min([1, volume/(Lambda**self.dimensions*(self.total_number_atoms))*np.exp(beta*(self.desired_mu-trial_Epot+initial_Epot))])
+ 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 np.random.random() < insert_or_delete:
+ self.successful_insert += 1
else:
- # 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.calculate_cross_coefficients()
- trial_Epot = self.compute_potential(output="potential")
- Lambda = self.calculate_Lambda(self.atom_mass[self.inserted_type])
- beta = 1/self.desired_temperature
- acceptation_probability = np.min([1, (Lambda**self.dimensions*(self.total_number_atoms-1)/volume)*np.exp(-beta*(self.desired_mu+trial_Epot-initial_Epot))])
- else:
- acceptation_probability = 0
- if np.random.random() < acceptation_probability:
- # Accept the new position
- pass
+ 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:
- # Reject the new position
- self.neighbor_lists = initial_neighbor_lists
- self.atoms_positions = initial_positions
- self.number_atoms = initial_number_atoms
- self.calculate_cross_coefficients()
+ 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:
+
+.. 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 relly 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 relly 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
@@ -76,8 +150,6 @@ Complete the *__init__* method as follows:
displace_mc = None,
desired_mu = None,
inserted_type = 0,
- desired_temperature = 300,
- (...)
.. label:: end_MonteCarlo_class
@@ -93,27 +165,30 @@ and
self.displace_mc = displace_mc
self.desired_mu = desired_mu
self.inserted_type = inserted_type
- self.desired_temperature = desired_temperature
- (...)
.. label:: end_MonteCarlo_class
-Let us non-dimentionalize desired_mu by adding:
+Let us also normalised the "desired_mu":
.. label:: start_MonteCarlo_class
.. code-block:: python
- def nondimensionalize_units_3(self):
- (...)
- self.displace_mc /= self.reference_distance
- if self.desired_mu is not None:
- self.desired_mu /= self.reference_energy
- (...)
+ class MonteCarlo(Outputs):
+ 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
-Finally, *monte_carlo_insert_delete* must be included in the run:
+Finally, the *monte_carlo_exchange()* method must be included in the run:
.. label:: start_MonteCarlo_class
@@ -121,14 +196,14 @@ Finally, *monte_carlo_insert_delete* must be included in the run:
def run(self):
(...)
- self.monte_carlo_displacement()
- self.monte_carlo_insert_delete()
+ self.monte_carlo_move()
+ self.monte_carlo_exchange()
self.wrap_in_box()
(...)
.. label:: end_MonteCarlo_class
-We need to calculate Lambda:
+We need to calculate :math:`\Lambda`:
.. label:: start_MonteCarlo_class
@@ -136,47 +211,92 @@ We need to calculate Lambda:
def calculate_Lambda(self, mass):
"""Estimate the de Broglie wavelength."""
- # Is it normal that mass is unitless ???
- m = mass/cst.Avogadro*cst.milli # kg
- kB = cst.Boltzmann # J/K
- Na = cst.Avogadro
- kB *= Na/cst.calorie/cst.kilo # kCal/mol/K
T = self.desired_temperature # N
- T_K = T*self.reference_energy/kB # K
- Lambda = cst.h/np.sqrt(2*np.pi*kB*m*T_K) # m
- Lambda /= cst.angstrom # Angstrom
- return Lambda / self.reference_distance # dimensionless
+ 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. Let us use a displace distance of
-0.5 Angstrom, and make 1000 steps.
+One can use a similar test as previously, but with an imposed chemical
+potential *desired_mu*:
-.. label:: start_test_MonteCarlo_class
+.. label:: start_test_8a_class
.. code-block:: python
- import os
from MonteCarlo import MonteCarlo
+ from pint import UnitRegistry
+ ureg = UnitRegistry()
+ import os
- mc = MonteCarlo(maximum_steps=1000,
- dumping_period=100,
- thermo_period=100,
- desired_mu = -3,
- inserted_type = 0,
- number_atoms=[50],
- epsilon=[0.1], # kcal/mol
- sigma=[3], # A
- atom_mass=[1], # g/mol
- box_dimensions=[20, 20, 20], # A
- )
+ # Define atom number of each group
+ nmb_1= 50
+ # Define LJ parameters (sigma)
+ sig_1 = 3*ureg.angstrom
+ # Define LJ parameters (epsilon)
+ eps_1 = 0.1*ureg.kcal/ureg.mol
+ # Define atom mass
+ mss_1 = 10*ureg.gram/ureg.mol
+ # Define box size
+ L = 20*ureg.angstrom
+ # Define a cut off
+ rc = 2.5*sig_1
+ # Pick the desired temperature
+ T = 300*ureg.kelvin
+ # choose the desired_mu
+ desired_mu = -3*ureg.kcal/ureg.mol
+
+ # Initialize the prepare object
+ mc = MonteCarlo(
+ 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
+ cut_off=rc,
+ thermo_outputs="Epot-press",
+ desired_temperature=T, # K
+ neighbor=1,
+ desired_mu = desired_mu,
+ )
mc.run()
-.. label:: end_test_MonteCarlo_class
+ # 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")
+
+ # 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
The evolution of the potential energy as a function of the
-number of steps are written in the *Outputs/Epot.dat* file
-and can be plotted.
+number of steps are written in the *Outputs/Epot.dat*.
diff --git a/docs/source/index.rst b/docs/source/index.rst
index cf64b19..f094c45 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -57,6 +57,7 @@ Learn Molecular Simulations with Python can be followed by anyone with a compute
chapters/chapter5
chapters/chapter6
chapters/chapter7
+ chapters/chapter8
.. toctree::
:maxdepth: 2
diff --git a/docs/source/journal-article.bib b/docs/source/journal-article.bib
index 7148bbc..6616f3b 100644
--- a/docs/source/journal-article.bib
+++ b/docs/source/journal-article.bib
@@ -61,4 +61,60 @@ @book{Rossum2009Python3
isbn = {1441412697},
publisher = {CreateSpace},
address = {Scotts Valley, CA}
-}
\ No newline at end of file
+}
+
+@article{ harris2020array,
+ title = {Array programming with {NumPy}},
+ author = {Charles R. Harris and K. Jarrod Millman and St{\'{e}}fan J.
+ van der Walt and Ralf Gommers and Pauli Virtanen and David
+ Cournapeau and Eric Wieser and Julian Taylor and Sebastian
+ Berg and Nathaniel J. Smith and Robert Kern and Matti Picus
+ and Stephan Hoyer and Marten H. van Kerkwijk and Matthew
+ Brett and Allan Haldane and Jaime Fern{\'{a}}ndez del
+ R{\'{i}}o and Mark Wiebe and Pearu Peterson and Pierre
+ G{\'{e}}rard-Marchant and Kevin Sheppard and Tyler Reddy and
+ Warren Weckesser and Hameer Abbasi and Christoph Gohlke and
+ Travis E. Oliphant},
+ year = {2020},
+ month = sep,
+ journal = {Nature},
+ volume = {585},
+ number = {7825},
+ pages = {357--362},
+ doi = {10.1038/s41586-020-2649-2},
+ publisher = {Springer Science and Business Media {LLC}},
+ url = {https://doi.org/10.1038/s41586-020-2649-2}
+}
+
+@article{michaud2011mdanalysis,
+ title={MDAnalysis: a toolkit for the analysis of molecular dynamics simulations},
+ author={Michaud-Agrawal, Naveen and Denning, Elizabeth J and Woolf, Thomas B and Beckstein, Oliver},
+ journal={Journal of computational chemistry},
+ volume={32},
+ number={10},
+ pages={2319--2327},
+ year={2011},
+ publisher={Wiley Online Library}
+}
+
+@article{humphrey1996vmd,
+ title={VMD: visual molecular dynamics},
+ author={Humphrey, William and Dalke, Andrew and Schulten, Klaus},
+ journal={Journal of molecular graphics},
+ volume={14},
+ number={1},
+ pages={33--38},
+ year={1996},
+ publisher={Elsevier}
+}
+
+@article{stukowski2009visualization,
+ title={Visualization and analysis of atomistic simulation data with OVITO--the Open Visualization Tool},
+ author={Stukowski, Alexander},
+ journal={Modelling and simulation in materials science and engineering},
+ volume={18},
+ number={1},
+ pages={015012},
+ year={2009},
+ publisher={IOP Publishing}
+}
diff --git a/tests/build-documentation.py b/tests/build-documentation.py
index 247c50a..b6c5395 100644
--- a/tests/build-documentation.py
+++ b/tests/build-documentation.py
@@ -18,7 +18,7 @@
os.mkdir("generated-codes/")
# loop on the different chapter
-for chapter_id in [1, 2, 3, 4, 5, 6, 7]:
+for chapter_id in [1, 2, 3, 4, 5, 6, 7, 8]:
# for each chapter, create the corresponding code
RST_EXISTS, created_tests, folder = sphinx_to_python(path_to_docs, chapter_id)
if RST_EXISTS: