From 7eade5434054d6e1cbfc10264600b9ff0d2a9ddc Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Thu, 29 Aug 2024 18:48:08 +0200 Subject: [PATCH 1/4] start adapting the code to 3D --- docs/source/chapters/chapter3.rst | 19 ++++++++++++------- docs/source/chapters/chapter4.rst | 5 ++++- docs/source/chapters/chapter6.rst | 5 ++++- 3 files changed, 20 insertions(+), 9 deletions(-) diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst index 2952ef3..2009b91 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,11 +138,16 @@ 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_size = np.diff(self.box_boundaries).reshape(self.dimension) box_geometry = np.array([90, 90, 90]) self.box_size = np.array(box_size.tolist()+box_geometry.tolist()) @@ -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..974d682 100644 --- a/docs/source/chapters/chapter4.rst +++ b/docs/source/chapters/chapter4.rst @@ -291,7 +291,10 @@ 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": 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") From 4d39713453819ec268dbbff6e3ba00e3160da775 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Thu, 29 Aug 2024 18:52:02 +0200 Subject: [PATCH 2/4] fix broken line --- docs/source/chapters/chapter3.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/chapters/chapter3.rst b/docs/source/chapters/chapter3.rst index 2009b91..4731e76 100644 --- a/docs/source/chapters/chapter3.rst +++ b/docs/source/chapters/chapter3.rst @@ -147,7 +147,7 @@ method to the *InitializeSimulation* class: box_boundaries[dim] = -L/2, L/2 dim += 1 self.box_boundaries = box_boundaries - box_size = np.diff(self.box_boundaries).reshape(self.dimension) + box_size = np.diff(self.box_boundaries).reshape(3) box_geometry = np.array([90, 90, 90]) self.box_size = np.array(box_size.tolist()+box_geometry.tolist()) From 88a706473a566acefc3b9d49bec7d30a09299f1b Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Thu, 29 Aug 2024 19:03:55 +0200 Subject: [PATCH 3/4] improved pressure calculation in 2D --- docs/source/chapters/chapter7.rst | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) 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 From 13132aeb721016794f145aa20d0e334bc99e7548 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Thu, 29 Aug 2024 22:42:26 +0200 Subject: [PATCH 4/4] implemented hard sphere potential --- docs/source/chapters/chapter1.rst | 40 ++++++++++++++++++++++--------- docs/source/chapters/chapter2.rst | 9 +++++-- docs/source/chapters/chapter4.rst | 4 ++-- 3 files changed, 38 insertions(+), 15 deletions(-) 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/chapter4.rst b/docs/source/chapters/chapter4.rst index 974d682..2a2df5e 100644 --- a/docs/source/chapters/chapter4.rst +++ b/docs/source/chapters/chapter4.rst @@ -298,9 +298,9 @@ class. 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