@@ -32,8 +32,10 @@ The steepest descent method is used for energy minimization, following these ste
32
32
- 3) Move the atoms in the opposite direction of the maximum
33
33
force to minimize the potential energy by a displacement step.
34
34
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}`.
37
39
38
40
- If :math: `\Delta E < 0 `, the new configuration is accepted as it results in
39
41
lower energy, and the step size is increased.
@@ -63,13 +65,8 @@ Let us fill the *__init__()* method:
63
65
64
66
.. label :: end_MinimizeEnergy_class
65
67
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.
73
70
74
71
Energy minimizer
75
72
----------------
@@ -89,8 +86,13 @@ following *run()* method to the *MinimizeEnergy* class:
89
86
self .update_neighbor_lists() # Rebuild neighbor list, if necessary
90
87
self .update_cross_coefficients() # Recalculate the cross coefficients, if necessary
91
88
# 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
94
96
# Save the current atom positions
95
97
init_positions = copy.deepcopy(self .atoms_positions)
96
98
# Move the atoms in the opposite direction of the maximum force
@@ -101,10 +103,9 @@ following *run()* method to the *MinimizeEnergy* class:
101
103
# Keep the more favorable energy
102
104
if trial_Epot < init_Epot: # accept new position
103
105
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
106
107
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
108
109
self .displacement *= 1.2 # Multiply the displacement by a factor 1.2
109
110
else : # reject new position
110
111
self .Epot = init_Epot # Revert to old energy
@@ -113,17 +114,19 @@ following *run()* method to the *MinimizeEnergy* class:
113
114
114
115
.. label :: end_MinimizeEnergy_class
115
116
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.
119
122
120
- Compute_potential
121
- -----------------
123
+ Compute the Potential
124
+ ---------------------
122
125
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:
127
130
128
131
.. label :: start_Utilities_class
129
132
@@ -139,16 +142,16 @@ class:
139
142
rij = self .compute_distance(self .atoms_positions[Ni],
140
143
self .atoms_positions[neighbor_of_i],
141
144
self .box_size)
142
- # Measure potential using information about cross coefficients
145
+ # Measure potential using pre-calculated cross coefficients
143
146
sigma_ij = self .sigma_ij_list[Ni]
144
147
epsilon_ij = self .epsilon_ij_list[Ni]
145
148
energy_potential += np.sum(potentials(epsilon_ij, sigma_ij, rij))
146
149
return energy_potential
147
150
148
151
.. label :: end_Utilities_class
149
152
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 *
152
155
class as well:
153
156
154
157
.. label :: start_Utilities_class
@@ -158,7 +161,7 @@ class as well:
158
161
def compute_distance (self ,position_i , positions_j , box_size , only_norm = True ):
159
162
"""
160
163
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?
162
165
"""
163
166
rij_xyz = np.nan_to_num(np.remainder(position_i - positions_j
164
167
+ box_size[:3 ]/ 2.0 , box_size[:3 ]) - box_size[:3 ]/ 2.0 )
@@ -171,7 +174,9 @@ class as well:
171
174
172
175
Finally, the energy minimization requires the computation of the minimum
173
176
force 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):
175
180
176
181
.. label :: start_Utilities_class
177
182
@@ -202,23 +207,20 @@ let us create a new method that is dedicated solely to measuring forces:
202
207
# Add the contribution to the matrix
203
208
force_matrix[Ni][neighbor_of_i] += (fij_xyz* rij_xyz.T/ rij).T
204
209
if return_vector:
205
- max_force = np.max(np.abs(force_vector))
206
- return force_vector, max_force
210
+ return force_vector
207
211
else :
208
212
return force_matrix
209
213
210
214
.. label :: end_Utilities_class
211
215
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.
216
218
217
- Wrap in box
219
+ Wrap in Box
218
220
-----------
219
221
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
222
224
within the *Utilities * class:
223
225
224
226
.. label :: start_Utilities_class
@@ -296,5 +298,6 @@ typically negative.
296
298
297
299
.. label :: end_test_4a_class
298
300
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