@@ -32,8 +32,10 @@ The steepest descent method is used for energy minimization, following these ste
3232- 3) Move the atoms in the opposite direction of the maximum
3333 force to minimize the potential energy by a displacement step.
3434 The size of the step is adjusted iteratively based on the reduction in energy.
35- - 4) Compute the new potential energy after the displacement, :math: `E_\text {pot}^\text {trial}`.
36- - 5) Evaluate the change in energy: :math: `\Delta E = E_\text {pot}^\text {trial} - E_\text {pot}^\text {initial}`.
35+ - 4) Compute the new potential energy after the displacement,
36+ :math: `E_\text {pot}^\text {trial}`.
37+ - 5) Evaluate the change in energy:
38+ :math: `\Delta E = E_\text {pot}^\text {trial} - E_\text {pot}^\text {initial}`.
3739
3840 - If :math: `\Delta E < 0 `, the new configuration is accepted as it results in
3941 lower energy, and the step size is increased.
@@ -63,13 +65,8 @@ Let us fill the *__init__()* method:
6365
6466 .. label :: end_MinimizeEnergy_class
6567
66- An important parameter is *maximum_steps *, which sets the maximum number
67- of steps for the energy minimization process. The *displacement *
68- parameter, with a default value of 0.01 Ångström, sets the initial atom
69- displacement value.
70-
71- The *thermo_outputs * and *data_folder * parameters are used for printing data
72- to files. These two parameters will be useful in the next chapter, :ref: `chapter5-label `.
68+ The parameter *maximum_steps * sets the maximum number
69+ of steps for the energy minimization process.
7370
7471Energy minimizer
7572----------------
@@ -89,8 +86,13 @@ following *run()* method to the *MinimizeEnergy* class:
8986 self .update_neighbor_lists() # Rebuild neighbor list, if necessary
9087 self .update_cross_coefficients() # Recalculate the cross coefficients, if necessary
9188 # Compute Epot/MaxF/force
92- init_Epot = self .compute_potential()
93- forces, init_MaxF = self .compute_force()
89+ if hasattr (self , ' Epot' ) is False : # If self.Epot does not exists yet, calculate it
90+ self .Epot = self .compute_potential()
91+ if hasattr (self , ' MaxF' ) is False : # If self.MaxF does not exists yet, calculate it
92+ forces = self .compute_force()
93+ self .MaxF = np.max(np.abs(forces))
94+ init_Epot = self .Epot
95+ init_MaxF = self .MaxF
9496 # Save the current atom positions
9597 init_positions = copy.deepcopy(self .atoms_positions)
9698 # Move the atoms in the opposite direction of the maximum force
@@ -101,10 +103,9 @@ following *run()* method to the *MinimizeEnergy* class:
101103 # Keep the more favorable energy
102104 if trial_Epot < init_Epot: # accept new position
103105 self .Epot = trial_Epot
104- # calculate the new max force and save it
105- forces, init_MaxF = self .compute_force()
106+ forces = self .compute_force() # calculate the new max force and save it
106107 self .MaxF = np.max(np.abs(forces))
107- self .wrap_in_box() # Wrap atoms in the box, if necessary
108+ self .wrap_in_box() # Wrap atoms in the box
108109 self .displacement *= 1.2 # Multiply the displacement by a factor 1.2
109110 else : # reject new position
110111 self .Epot = init_Epot # Revert to old energy
@@ -113,17 +114,19 @@ following *run()* method to the *MinimizeEnergy* class:
113114
114115 .. label :: end_MinimizeEnergy_class
115116
116- The displacement, which has an initial value of 0.01, is adjusted through energy
117- minimization. When the trial is successful, its value is multiplied by 1.2. When
118- the trial is rejected, its value is multiplied by 0.2.
117+ The displacement, which has an initial value of 0.01, is adjusted through
118+ energy minimization. When the trial is successful, its value is multiplied
119+ by 1.2. When the trial is rejected, its value is multiplied by 0.2. Thus,
120+ as the minimization progresses, the displacements that the particles
121+ perform increase.
119122
120- Compute_potential
121- -----------------
123+ Compute the Potential
124+ ---------------------
122125
123- Computing the potential energy of the system is central to the energy minimizer,
124- as the value of the potential is used to decide if the trial is accepted or
125- rejected. Add the following method called *compute_potential() * to the * Utilities *
126- class:
126+ Computing the potential energy of the system is central to the energy
127+ minimizer, as the value of the potential is used to decide if the trial is
128+ accepted or rejected. Add the following method called *compute_potential() *
129+ to the * Utilities * class:
127130
128131.. label :: start_Utilities_class
129132
@@ -139,16 +142,16 @@ class:
139142 rij = self .compute_distance(self .atoms_positions[Ni],
140143 self .atoms_positions[neighbor_of_i],
141144 self .box_size)
142- # Measure potential using information about cross coefficients
145+ # Measure potential using pre-calculated cross coefficients
143146 sigma_ij = self .sigma_ij_list[Ni]
144147 epsilon_ij = self .epsilon_ij_list[Ni]
145148 energy_potential += np.sum(potentials(epsilon_ij, sigma_ij, rij))
146149 return energy_potential
147150
148151 .. label :: end_Utilities_class
149152
150- Measuring the distance is an important step of computing the potential. Let us
151- do it using a dedicated method. Add the following method to the *Utilities *
153+ Measuring the distance is a central step in computing the potential. Let us
154+ do this using a dedicated method. Add the following method to the *Utilities *
152155class as well:
153156
154157.. label :: start_Utilities_class
@@ -158,7 +161,7 @@ class as well:
158161 def compute_distance (self ,position_i , positions_j , box_size , only_norm = True ):
159162 """
160163 Measure the distances between two particles.
161- # TOFIX: Move as function instead of a method?
164+ # TOFIX: Move as a function instead of a method?
162165 """
163166 rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j
164167 + box_size[:3 ]/ 2.0 , box_size[:3 ]) - box_size[:3 ]/ 2.0 )
@@ -171,7 +174,9 @@ class as well:
171174
172175Finally, the energy minimization requires the computation of the minimum
173176force in the system. Although not very different from the potential measurement,
174- let us create a new method that is dedicated solely to measuring forces:
177+ let us create a new method that is dedicated solely to measuring forces.
178+ The force can be returned as a vector (default), or as a matrix (which will be
179+ useful for pressure calculation, see the :ref: `chapter7-label ` chapter):
175180
176181.. label :: start_Utilities_class
177182
@@ -202,23 +207,20 @@ let us create a new method that is dedicated solely to measuring forces:
202207 # Add the contribution to the matrix
203208 force_matrix[Ni][neighbor_of_i] += (fij_xyz* rij_xyz.T/ rij).T
204209 if return_vector:
205- max_force = np.max(np.abs(force_vector))
206- return force_vector, max_force
210+ return force_vector
207211 else :
208212 return force_matrix
209213
210214 .. label :: end_Utilities_class
211215
212- Here, two types of outputs can
213- be requested by the user: *force-vector *, and *force-matrix *.
214- The *force-matrix * option will be useful for pressure calculation, see
215- :ref: `chapter7-label `.
216+ The force can be returned as a vector (default) or as a matrix. The latter
217+ will be useful for pressure calculation; see the :ref: `chapter7-label ` chapter.
216218
217- Wrap in box
219+ Wrap in Box
218220-----------
219221
220- Every time atoms are being displaced, one has to ensure that they remain in
221- the box. This is done by the *wrap_in_box() * method that must be placed
222+ Every time atoms are displaced, one has to ensure that they remain within
223+ the box. This is done by the *wrap_in_box() * method, which must be placed
222224within the *Utilities * class:
223225
224226.. label :: start_Utilities_class
@@ -296,5 +298,6 @@ typically negative.
296298
297299 .. label :: end_test_4a_class
298300
299- For such as low density in particle, we can reasonably expect the energy to be always
300- negative after 100 steps.
301+ For such a low particle density, we can reasonably expect that the potential
302+ energy will always be negative after 100 steps, and that the maximum force
303+ will be smaller than 10 (unitless).
0 commit comments