diff --git a/docs/source/chapters/chapter1.rst b/docs/source/chapters/chapter1.rst index a130331..599d46e 100644 --- a/docs/source/chapters/chapter1.rst +++ b/docs/source/chapters/chapter1.rst @@ -66,33 +66,51 @@ between atoms. Although more complicated 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 named *LJ_potential* for the Lennard-Jones -potential (LJ). Copy the following into *potentials.py*: +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. + +Copy the following lines into *potentials.py*: .. label:: start_potentials_class .. code-block:: python - def LJ_potential(epsilon, sigma, r, derivative = False): - if derivative: - return 48*epsilon*((sigma/r)**12-0.5*(sigma/r)**6)/r + 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: + raise ValueError("Derivative is not defined for Hard-Sphere potential.") + else: + return 1000 if r < sigma else 0 else: - return 4*epsilon*((sigma/r)**12-(sigma/r)**6) + raise ValueError(f"Unknown potential type: {potential_type}") .. label:: end_potentials_class -Depending on the value of the optional argument *derivative*, which can be -either *False* or *True*, the *LJ_potential* function will return the force: +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: .. 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], 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]. Create the Classes ------------------ @@ -124,7 +142,7 @@ copy the following lines: .. code-block:: python - from potentials import LJ_potential + from potentials import potentials class Utilities: diff --git a/docs/source/chapters/chapter2.rst b/docs/source/chapters/chapter2.rst index 7d2f9d6..f297513 100644 --- a/docs/source/chapters/chapter2.rst +++ b/docs/source/chapters/chapter2.rst @@ -99,12 +99,14 @@ Modify the *Prepare* class as follows: epsilon=[0.1], # List - Kcal/mol sigma=[3], # List - Angstrom atom_mass=[10], # List - g/mol + potential_type="Lennard-Jones", *args, **kwargs): 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 @@ -114,7 +116,10 @@ Here, the four lists *number_atoms* :math:`N`, *epsilon* :math:`\epsilon`, :math:`10`, :math:`0.1~\text{[Kcal/mol]}`, :math:`3~\text{[Å]}`, and :math:`10~\text{[g/mol]}`, respectively. -All four parameters are assigned to *self*, allowing other methods to access +The type of potential is also specified, with Lennard-Jones being closen as +the default option. + +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 of positional and keyword arguments, respectively. @@ -130,7 +135,7 @@ unit system to the *LJ* unit system: .. code-block:: python def calculate_LJunits_prefactors(self): - """Calculate the Lennard-Jones units prefacors.""" + """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 diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst index 2952ef3..4731e76 100644 --- a/docs/source/chapters/chapter3.rst +++ b/docs/source/chapters/chapter3.rst @@ -35,7 +35,8 @@ Let us improve the previously created *InitializeSimulation* class: ): super().__init__(*args, **kwargs) self.box_dimensions = box_dimensions - self.dimensions = len(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.initial_positions = initial_positions self.thermo_period = thermo_period @@ -137,9 +138,14 @@ method to the *InitializeSimulation* class: .. code-block:: python def define_box(self): - box_boundaries = np.zeros((self.dimensions, 2)) - for dim, L in zip(range(self.dimensions), self.box_dimensions): + """Define the simulation box. + For 2D simulations, the third dimensions only contains 0. + """ + box_boundaries = np.zeros((3, 2)) + dim = 0 + for L in 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]) @@ -183,9 +189,8 @@ case, the array must be of size 'number of atoms' times 'number of dimensions'. def populate_box(self): if self.initial_positions is None: - atoms_positions = np.zeros((self.total_number_atoms, - self.dimensions)) - for dim in np.arange(self.dimensions): + atoms_positions = np.zeros((self.total_number_atoms, 3)) + for dim in np.arange(3): diff_box = np.diff(self.box_boundaries[dim]) random_pos = np.random.random(self.total_number_atoms) atoms_positions[:, dim] = random_pos*diff_box-diff_box/2 diff --git a/docs/source/chapters/chapter4.rst b/docs/source/chapters/chapter4.rst index e8fc0f0..2a2df5e 100644 --- a/docs/source/chapters/chapter4.rst +++ b/docs/source/chapters/chapter4.rst @@ -291,13 +291,16 @@ class. sigma_ij = self.sigma_ij_list[Ni] epsilon_ij = self.epsilon_ij_list[Ni] # Measure distances - rij_xyz = (np.remainder(position_i - positions_j + half_box_size, box_size) - half_box_size) + # Measure distances + # The nan_to_num is crutial in 2D to avoid nan value along third dimension + rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j + + half_box_size, box_size) - half_box_size) rij = np.linalg.norm(rij_xyz, axis=1) # Measure potential if output == "potential": - energy_potential += np.sum(LJ_potential(epsilon_ij, sigma_ij, rij)) + energy_potential += np.sum(potentials(self.potential_type, epsilon_ij, sigma_ij, rij)) else: - derivative_potential = LJ_potential(epsilon_ij, sigma_ij, rij, derivative = True) + derivative_potential = potentials(self.potential_type, epsilon_ij, sigma_ij, rij, derivative = True) if output == "force-vector": forces[Ni] += np.sum((derivative_potential*rij_xyz.T/rij).T, axis=0) forces[neighbor_of_i] -= (derivative_potential*rij_xyz.T/rij).T diff --git a/docs/source/chapters/chapter6.rst b/docs/source/chapters/chapter6.rst index 4c3b93a..49b18ca 100644 --- a/docs/source/chapters/chapter6.rst +++ b/docs/source/chapters/chapter6.rst @@ -57,7 +57,10 @@ Let us add a method named *monte_carlo_move* to the *MonteCarlo* class: atom_id = np.random.randint(self.total_number_atoms) # Move the chosen atom in a random direction # The maximum displacement is set by self.displace_mc - move = (np.random.random(self.dimensions)-0.5)*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) self.atoms_positions[atom_id] += move # Measure the optential energy of the new configuration trial_Epot = self.compute_potential(output="potential") diff --git a/docs/source/chapters/chapter7.rst b/docs/source/chapters/chapter7.rst index a430494..68cdfa2 100644 --- a/docs/source/chapters/chapter7.rst +++ b/docs/source/chapters/chapter7.rst @@ -39,16 +39,16 @@ Let us add the following method to the *Utilities* class. def calculate_pressure(self): """Evaluate p based on the Virial equation (Eq. 4.4.2 in Frenkel-Smit, Understanding molecular simulation: from algorithms to applications, 2002)""" - # Ideal contribution + # Compute the ideal contribution Ndof = self.dimensions*self.total_number_atoms-self.dimensions - volume = np.prod(self.box_size[:3]) + volume = np.prod(self.box_size[:self.dimensions]) 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) - # Non-ideal contribution + # Compute the non-ideal contribution distances_forces = np.sum(self.compute_potential(output="force-matrix")*self.evaluate_rij_matrix()) p_nonideal = distances_forces/(volume*self.dimensions) # Final pressure @@ -73,7 +73,8 @@ To evaluate all the vectors between all the particles, let us also add the position_i = self.atoms_positions[Ni] rij_xyz = (np.remainder(position_i - positions_j + half_box_size, box_size) - half_box_size) rij_matrix[Ni] = rij_xyz - return rij_matrix + # use nan_to_num to avoid "nan" values in 2D + return np.nan_to_num(rij_matrix) .. label:: end_Utilities_class