@@ -140,7 +140,9 @@ lines to **free-sampling.lmp**:
140140 variable F atom ${U0}/((x-${x0})^2/${dlt}^2+1)/${dlt}-${U0}/((x+${x0})^2/${dlt}^2+1)/${dlt}
141141 fix myadf all addforce v_F 0.0 0.0 energy v_U
142142
143- Next, we combine the ``fix nve `` with a ``fix langevin `` thermostat:
143+ Next, we use the Newtonian equations of motion with
144+ a Langevin thermostat by combining the ``fix nve `` with a
145+ ``fix langevin `` command:
144146
145147.. code-block :: lammps
146148
@@ -152,6 +154,13 @@ in the NVT ensemble, maintaining a constant number of
152154atoms :math: `N`, constant volume :math: `V`, and a temperature :math: `T` that
153155fluctuates around a target value.
154156
157+ .. admonition :: Note
158+ :class: non-title-info
159+
160+ LAMMPS documentation suggests using damping constants for thermostats that
161+ are approximately 100 times the timestep value. In this case, a value of
162+ 500 is used, resulting in a relatively weak coupling to the thermostat.
163+
155164To ensure that the equilibration time is sufficient, we will track the evolution of
156165the number of atoms in the central - energetically unfavorable - region,
157166referred to as ``mymes ``, using the ``n_center `` variable:
@@ -224,7 +233,7 @@ Add the following line to **free-sampling.lmp**:
224233 Here, the ``chunk/atom `` command discretizes the simulation
225234domain into spatial bins of size 2~\A A{} along the :math: `x` direction,
226235and the ``ave/chunk `` command computes and outputs the number density of
227- atoms within each bin to the file **free-sampling.dat **.}
236+ atoms within each bin to the file **free-sampling.dat **.
228237The step count is reset to 0 using ``reset_timestep `` to synchronize it
229238with the output times of ``fix density/number ``. Run the simulation using
230239LAMMPS.
@@ -253,7 +262,7 @@ Next, we plot :math:`-R T \ln(\rho/\rho_\mathrm{bulk})`, where :math:`\rho/\rho_
253262is the the density ratio, and compare it
254263with the imposed potential :math: `U` from Eq. :eq: `eq_U `.
255264The reference density, :math: `\rho _\text {bulk} = 0.0009 ~\text {Å}^{-3 }`,
256- was estimated by measuring the density of the reservoir from the raw density
265+ was estimated by measuring the density of the reservoir from the density
257266profiles. The agreement between the MD results and the imposed energy profile
258267is excellent, despite some noise in the central part, where fewer data points
259268are available due to the repulsive potential.
@@ -278,8 +287,8 @@ The limits of free sampling
278287---------------------------
279288
280289Increasing the value of :math: `U_0 ` reduces the average number of atoms in the central
281- region, making it difficult to achieve a high-resolution free energy profile.
282- For example, running the same simulation with :math: `U_0 = 10 \epsilon `,
290+ region, making it difficult to achieve a high-resolution free energy profile
291+ within reasonable simulation times. For example, running the same simulation with :math: `U_0 = 10 \epsilon `,
283292corresponding to :math: `U_0 \approx 10 k_\text {B} T`, results in no atoms exploring
284293the central part of the simulation box during the simulation.
285294In such a case, employing an enhanced sampling method is recommended, as done in the next section.
@@ -381,7 +390,7 @@ So far, our code resembles that of Method 1, except for the additional particle
381390of type 2. Particles of types 1 and 2 are identical, with the same mass
382391and LJ parameters. However, the particle of type 2 will also
383392be exposed to the biasing potential :math: `V`, which forces it to explore the
384- central part of the box.
393+ central part of the box, thus justifying the definition of two atom types .
385394
386395..
387396 TOFIX: Add a figure with one single particle exploring the central part of the system.
@@ -412,9 +421,13 @@ bias potential by increments of 0.4 nm. Add the following lines to **umbrella-s
412421 next a
413422 jump SELF loop
414423
424+ The definition of a variable of loop style serves the same purpose as in
425+ :ref: `reactive-silicon-dioxide-label `, and we highlight here the particular
426+ utility of using its value to distinguish the files written by the
427+ fix ``ave_time `` command for the different bias potentials.
415428The ``spring `` command imposes the additional harmonic potential :math: `V` with
416- the previously defined spring constant :math: `k`. The center of the harmonic
417- potential, :math: `x_\text {des}`, successively takes values
429+ the previously defined spring constant :math: `k` to the atoms in the group `` pull ``.
430+ The center of the harmonic potential, :math: `x_\text {des}`, successively takes values
418431from :math: `-28 \,\text {Å}` to :math: `28 \,\text {Å}`. For each value of :math: `x_\text {des}`,
419432an equilibration step of 40 ps is performed, followed by a step
420433of 400 ps during which the position of the particle of
0 commit comments