diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst index 5540209..632d616 100644 --- a/docs/source/chapters/chapter1.rst +++ b/docs/source/chapters/chapter1.rst @@ -39,32 +39,32 @@ 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.py - *Utilities* class: General-purpose methods, inherited by all other classes - * - *InitializeSimulation.py* + * - 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 - * - *logger.py* + * - logger.py - Functions for outputting data into text files - * - *dumper.py* + * - dumper.py - Functions for outputting data into trajectory files for visualization - * - *reader.py* + * - reader.py - Functions for importing data from text files Some of these files are created in this chapter; others will be created later @@ -73,14 +73,13 @@ on. All of these files must be created within the same folder. 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. -Create the first file named *potentials.py*. This file will contain a function -called *potentials*. For the moment, the only potential that can be returned by -this function is the Lennard-Jones potential (LJ). This may change in the -future. +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*: @@ -88,16 +87,11 @@ 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) + 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 @@ -108,14 +102,19 @@ 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], + - \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]. + - \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 ------------------ @@ -183,8 +182,14 @@ Within the *InitializeSimulation.py* file, copy the following lines: .. label:: end_InitializeSimulation_class The *InitializeSimulation* class inherits from the previously created -*Prepare* and Utilities classes. 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: @@ -204,13 +209,13 @@ 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, named *MinimizeEnergy*, -*MonteCarlo*, and *MolecularDynamics*. Each of these three classes inherits -from the *Measurements* class, and thus from the classes inherited by -*Measurements*. +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). Within the *MinimizeEnergy.py* file, copy the following lines: @@ -219,6 +224,8 @@ Within the *MinimizeEnergy.py* file, copy the following lines: .. code-block:: python from Measurements import Measurements + import numpy as np + import copy import os @@ -230,8 +237,8 @@ Within the *MinimizeEnergy.py* file, copy the following lines: .. 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: @@ -239,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 @@ -257,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: @@ -331,9 +334,9 @@ Alternatively, this test can also be launched using *Pytest* by typing in a term 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 f6682d4..2714e19 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`, @@ -100,7 +96,6 @@ Modify the *Prepare* class as follows: epsilon, # List - Kcal/mol sigma, # List - Angstrom atom_mass, # List - g/mol - potential_type="Lennard-Jones", *args, **kwargs): self.ureg = ureg @@ -108,18 +103,19 @@ Modify the *Prepare* class as follows: 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. + +.. |Pint| raw:: html -The type of potential is also specified, with Lennard-Jones being chosen 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 diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst index 0f1a439..0746a69 100644 --- a/docs/source/chapters/chapter4.rst +++ b/docs/source/chapters/chapter4.rst @@ -47,19 +47,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 @@ -154,8 +142,7 @@ class: # Measure potential using information about 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 @@ -207,8 +194,7 @@ let us create a new method that is dedicated solely to measuring forces: # 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) diff --git a/docs/source/chapters/chapter8.rst b/docs/source/chapters/chapter8.rst index e1bf6cc..e4c9b8c 100644 --- a/docs/source/chapters/chapter8.rst +++ b/docs/source/chapters/chapter8.rst @@ -215,7 +215,6 @@ We need to calculate Lambda: def calculate_Lambda(self, mass): """Estimate the de Broglie wavelength.""" - m = mass/cst.Avogadro*cst.milli # kg T = self.desired_temperature # N return 1/np.sqrt(2*np.pi*mass*T) diff --git a/docs/source/journal-article.bib b/docs/source/journal-article.bib index 7148bbc..8869996 100644 --- a/docs/source/journal-article.bib +++ b/docs/source/journal-article.bib @@ -61,4 +61,27 @@ @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} +} \ No newline at end of file