@@ -102,7 +102,7 @@ for both the PEG (C, OE, H, OT, and HT atoms) and the water molecules (OW and HW
102102 angles, and dihedrals contraints.
103103
104104Create an **inputs/ ** folder next to **ff/ **, and create a new empty file
105- called **em.mdp **. Copy the following lines into it :
105+ called **em.mdp ** into it . Copy the following lines into ** em.mdp ** :
106106
107107.. code-block :: bw
108108
@@ -125,6 +125,15 @@ most important command is ``integrator = steep``, which sets the algorithm
125125used by GROMACS as the steepest-descent method. This algorithm moves the
126126atoms following the direction of the largest forces until one of the stopping
127127criteria is reached :cite: `debyeNaeherungsformelnFuerZylinderfunktionen1909 `.
128+ The stopping criteria include reaching the maximum force tolerance set by
129+ ``emtol ``, or completing the maximum number of steps specified by
130+ ``nsteps ``. The initial step size for energy minimization is controlled by
131+ ``emstep ``. Other parameters define the treatment of nonbonded interactions,
132+ such as ``cutoff-scheme = Verlet ``, electrostatics handled by Particle-Mesh
133+ Ewald with ``coulombtype = PME ``, cutoff distances ``rcoulomb `` and
134+ ``rvdw `` set to 1 nm. Periodic boundary conditions are applied in all
135+ directions with ``pbc = xyz ``. Output frequencies for energy and coordinate
136+ writing are set by ``nstenergy `` and ``nstxout ``, respectively.
128137
129138Run the energy minimization using GROMACS by typing in a terminal:
130139
@@ -136,43 +145,40 @@ Run the energy minimization using GROMACS by typing in a terminal:
136145 The ``-nt 8 `` option limits the number of threads that GROMACS uses. Adjust
137146the number to your computer.
138147
139- After the simulation is over , open the trajectory file with VMD by typing
140- in a terminal:
148+ After the simulation is complete , open the trajectory file with VMD by typing
149+ the following command in a terminal:
141150
142151.. code-block :: bash
143152
144153 vmd peg.gro em-peg.trr
145154
146- From VMD, the PEG molecule can be seen moving a little by the
147- steepest-descent algorithm .
155+ In VMD, you can observe the PEG molecule moving slightly as a result of the
156+ steepest-descent energy minimization .
148157
149- Before adding the water, let us reshape the box and recenter the PEG
150- molecule in the box. As a first step, let us use a cubic box of lateral
151- size :math: `2.6 ~\text {nm}`.
158+ Before adding water, let us reshape the box and recenter the PEG molecule within
159+ it. Let us also place it in a cubic box with a lateral size of :math: `2.6 ~\text {nm}`.
152160
153161.. code-block :: bash
154162
155163 gmx trjconv -f em-peg.gro -s em-peg.tpr -o peg-recentered.gro -center
156164 -pbc mol -box 2.6 2.6 2.6
157165
158- Select ``system `` for both centering and output. The newly created ** .gro **
159- file named **peg-recentered.gro ** will be used as a starting point for the
160- next step of the tutorial.
166+ Select ``system `` for both the centering and output prompts . The newly created
167+ **peg-recentered.gro ** file will be used as the starting point for the next step
168+ of the tutorial.
161169
162170Solvate the PEG molecule
163171========================
164172
165- Let us add the water molecules to the system by using * gmx solvate * :
173+ Let us add water molecules to the system using `` gmx solvate `` :
166174
167175.. code-block :: bash
168176
169177 gmx solvate -cp peg-recentered.gro -cs spc216.gro -o peg-solvated.gro -p topol.top
170178
171- Here *spc216.gro * is a default GROMACS file containing a pre-equilibrated
172- water reservoir.
173-
174- The newly created file *peg-solvated.gro * contains the water molecules,
175- and a a new line in was added to the topology file *topol.top *:
179+ Here, **spc216.gro ** is a default GROMACS file containing a pre-equilibrated
180+ water reservoir. The newly created file **peg-solvated.gro ** contains the
181+ water molecules, and a new line has been added to the topology file **topol.top **:
176182
177183.. code-block :: bw
178184
@@ -181,22 +187,17 @@ and a a new line in was added to the topology file *topol.top*:
181187 PEG 1
182188 SOL 546
183189
184- We can apply the same energy minimization to the newly created solvated
185- system. Simply add the following line to *em.mdp *:
190+ We can apply the same energy minimization as before to the newly created
191+ solvated system. Simply add the following line to ** em.mdp * *:
186192
187193.. code-block :: bw
188194
189195 define = -DFLEXIBLE
190196
191- And then launch the energy minimization again using:
192-
193- .. code-block :: bash
194-
195- gmx grompp -f inputs/em.mdp -c peg-solvated.gro -p topol.top -o em
196- gmx mdrun -deffnm em -v -nt 8
197-
198- The ``define = -DFLEXIBLE `` option triggers the following **if ** condition
199- within the **tip3p.itp ** file:
197+ With the ``define = -DFLEXIBLE `` option, the water molecules are treated as
198+ flexible during energy minimization, enabling bond stretching and angle bending
199+ in water. The ``define = -DFLEXIBLE `` option triggers the following **if **
200+ condition within the **tip3p.itp ** file:
200201
201202.. code-block :: bw
202203
@@ -209,22 +210,33 @@ within the **tip3p.itp** file:
209210 [ angles ]
210211 ; i j k funct angle force.c.
211212 2 1 3 1 104.52 628.02 104.52 628.02
212-
213- With this **if ** condition the water molecules
214- behave as flexible. This is better because rigid molecules and
215- energy minimization usually don't go along well. For the next molecular
216- dynamics steps, rigid water molecules will be used by not including
217- the ``define = -DFLEXIBLE `` command in the inputs.
213+
214+ With this **if ** condition, the water molecules behave as flexible, which is
215+ preferable because rigid molecules and energy minimization usually do not work
216+ well together.
217+
218+ .. admonition :: Note
219+ :class: non-title-info
220+
221+ For the subsequent molecular dynamics steps, rigid water
222+ molecules will be used by removing the ``define = -DFLEXIBLE `` command from the
223+ inputs.
224+
225+ Finally, launch the energy minimization again using:
226+
227+ .. code-block :: bash
228+
229+ gmx grompp -f inputs/em.mdp -c peg-solvated.gro -p topol.top -o em
230+ gmx mdrun -deffnm em -v -nt 8
218231
219232 Equilibrate the PEG-water system
220233================================
221234
222- Let use equilibrate the system in two steps: first a NVT simulation,
223- with constant number of particles, constant volume, and imposed temperature,
224- and second a NPT simulation with imposed pressure.
225-
226- Within the **inputs/ ** folder, create a new input named **nvt-peg-h2o.mdp **,
227- and copy the following lines into it:
235+ Let us further equilibrate the system in two steps: first, an *NVT * simulation
236+ with constant number of particles, constant volume, and imposed temperature;
237+ and second, an *NpT * simulation with imposed pressure. Within the **inputs/ **
238+ folder, create a new input file named **nvt-peg-h2o.mdp **, and copy the
239+ following lines into it:
228240
229241.. code-block :: bw
230242
@@ -260,11 +272,15 @@ and copy the following lines into it:
260272 comm-grps = PEG
261273
262274 Most of these commands have already been seen. In addition to the conventional
263- *md * leap-frog algorithm integrator, long-range Coulomb and short-range
264- van der Waals interactions, the LINCS constraint algorithm is used to maintain
265- the hydrogen bonds as rigid. An initial temperature of :math: `300 ~K` is given
266- to the system by the ``gen- `` commands, and the PEG is maintained in the center
267- of the box by the ``comm-mode `` and ``comm-grps `` commands.
275+ ``md `` leap-frog algorithm integrator, with long-range Coulomb and short-range
276+ van der Waals interactions, the LINCS constraint algorithm is used to keep the
277+ hydrogen bonds rigid. Temperature coupling at :math: `300 ~K` is imposed on both
278+ the PEG and water groups using velocity rescaling with a stochastic term
279+ (``tcoupl = v-rescale ``). Initial velocities at :math: `300 ~K` are generated
280+ by the ``gen- `` commands, with a specified random seed.
281+ The ``comm-mode `` and ``comm-grps `` commands ensure that the PEG molecule
282+ remains centered in the box by removing its center-of-mass motion separately
283+ from the solvent.
268284
269285Launch the *NVT * simulation using:
270286
@@ -273,23 +289,23 @@ Launch the *NVT* simulation using:
273289 gmx grompp -f inputs/nvt-peg-h2o.mdp -c em.gro -p topol.top -o nvt -maxwarn 1
274290 gmx mdrun -deffnm nvt -v -nt 8
275291
276- The ``maxwarn 1 `` option is used to avoid a GROMACS WARNING related to the
277- centering of the PEG in the box.
292+ The ``- maxwarn 1 `` option is used to bypass a GROMACS warning related to the
293+ centering of the PEG molecule in the box.
278294
279- Let us follow- up with the NPT equilibration. Duplicate the ** nvt-peg-h2o.mdp **
280- file into a new input file named **npt-peg-h2o.mdp **. Within ** npt-peg-h2o.mdp **,
281- Within the **npt-peg-h2o.mdp **, delete the lines related to the creation
282- of velocity as its better to keep the velocities generated during the
283- *NVT * run:
295+ Let us follow up with the * NPT * equilibration. Duplicate the
296+ ** nvt-peg-h2o.mdp ** file into a new input file named **npt-peg-h2o.mdp **.
297+ Within **npt-peg-h2o.mdp **, delete the lines related to the creation of
298+ velocities, as it is better to preserve the velocities generated during the
299+ previous *NVT * run:
284300
285301.. code-block :: bw
286302
287303 gen_vel = yes
288304 gen-temp = 300
289305 gen-seed = 65823
290306
291- In addition to the removal the previous 3 lines, add the following lines
292- to **npt-peg-h2o.mdp ** to specify the isotropic barostat with imposed pressure
307+ In addition to removing these three lines, add the following lines to
308+ **npt-peg-h2o.mdp ** to enable an isotropic barostat with an imposed pressure
293309of :math: `1 ~\text {bar}`:
294310
295311.. code-block :: bw
@@ -300,6 +316,16 @@ of :math:`1~\text{bar}`:
300316 ref-p = 1.0
301317 compressibility = 4.5e-5
302318
319+ These last five lines configure the pressure coupling. The ``pcoupl `` option selects
320+ the ``c-rescale `` barostat, a stochastic barostat suitable for equilibration. The
321+ pressure coupling type is set to isotropic, meaning the box dimensions scale
322+ uniformly in all directions. The ``tau-p `` parameter defines the coupling time
323+ constant in picoseconds, determining how quickly the system adjusts to the
324+ target pressure. The ``ref-p `` sets the reference pressure at
325+ :math: `1 ~\text {bar}`, and ``compressibility `` specifies the isothermal
326+ compressibility of water, which is required for volume fluctuations to properly
327+ reflect the solvent properties.
328+
303329Run the *NpT * simulation, using the final state of the *NVT * simulation
304330**nvt.gro ** as starting configuration:
305331
@@ -309,16 +335,16 @@ Run the *NpT* simulation, using the final state of the *NVT* simulation
309335 ${gmx} mdrun -deffnm npt -v -nt 8
310336
311337 Let us observe the evolution of the potential energy of the system during the
312- 3 successive equilibration steps, i.e. the energy minimization, *NVT *, and * NpT * steps,
313- using the ``gmx energy `` command as follow :
338+ three successive equilibration steps, i.e. the energy minimization, *NVT *, and
339+ * NpT * steps, using the ``gmx energy `` command as follows :
314340
315341.. code-block :: bash
316342
317343 gmx energy -f em.edr -o energy-em.xvg
318344 gmx energy -f nvt.edr -o energy-nvt.xvg
319345 gmx energy -f npt.edr -o energy-npt.xvg
320346
321- For each of the 3 ``gmx energy `` commands, select ``potential `` .
347+ For each of the three ``gmx energy `` commands, select ``Potential `` when prompted .
322348
323349.. figure :: ../figures/level2/stretching-a-polymer/potential-energy-light.png
324350 :alt: Potential energy from molecular dynamics simulation in GROMACS
@@ -334,13 +360,11 @@ For each of the 3 ``gmx energy`` commands, select ``potential``.
334360 respectively the energy minimization (a), the NVT step (b), and the NPT
335361 step (c).
336362
337- Let us launch a longer simulation, and extract the angle distribution
338- between the different atoms of the PEG molecules. This angle
339- distribution will be used later as a benchmark to probe the effect of
340- of the stretching on the PEG structure.
341-
342- Create a new input named **production-peg-h2o.mdp **, and copy the following
343- lines into it:
363+ Let us launch a longer simulation and extract the angle distribution
364+ between the different atoms of the PEG molecule. This angle distribution
365+ will later serve as a benchmark to assess the effect of stretching on the
366+ PEG structure. Create a new input file named **production-peg-h2o.mdp **, and
367+ copy the following lines into it:
344368
345369.. code-block :: bw
346370
@@ -371,8 +395,9 @@ lines into it:
371395 comm-mode = linear
372396 comm-grps = PEG
373397
374- This script resembles the **nvt-peg-h2o.mdp ** input, but the duration and
375- output frequency is different, and without the ``gen-vel `` commands.
398+ This input file is similar to **nvt-peg-h2o.mdp **, but with a longer simulation
399+ time and more frequent output, and without the ``gen-vel `` commands to preserve
400+ the velocities from the previous equilibration.
376401
377402Run it using:
378403
@@ -388,10 +413,10 @@ command:
388413
389414 gmx mk_angndx -s production.tpr -hyd no
390415
391- The **angle.ndx ** file generated contains groups with all the atoms
392- involved by an angle constraint , with the exception of the hydrogen
393- atoms due to the use of ``-hyd no ``. The atom ids selected in the groups
394- can be seen from the **index.ndx ** file:
416+ The **angle.ndx ** file generated contains groups of all atoms involved
417+ in angle constraints , with hydrogen atoms excluded thanks to the
418+ ``-hyd no `` option . The atom indices included in the groups can be
419+ verified in the **index.ndx ** file:
395420
396421.. code-block :: bw
397422
@@ -400,9 +425,9 @@ can be seen from the **index.ndx** file:
400425 31 33 35 38 40 42 45 47 49 52 54 56
401426 59 61 63 66 68 70 73 75 77 80 82 84
402427
403- Here, each number corresponds to the atom index, as can be seen from the
404- initial **peg.gro ** file. For instance, the atom of ``id 2 `` is a carbon atom,
405- and the atom with ``id 5 `` is an oxygen:
428+ Here, each number corresponds to an atom index, as listed in the
429+ initial **peg.gro ** file. For example, atom ``id 2 `` is a carbon atom,
430+ and atom ``id 5 `` is an oxygen atom :
406431
407432.. code-block :: bw
408433
@@ -419,13 +444,13 @@ and the atom with ``id 5`` is an oxygen:
419444 (...)
420445
421446 Then, extract the angle distribution from the **production.xtc **
422- file using ``gmx angle ``:
447+ trajectory using ``gmx angle ``:
423448
424449.. code-block :: bash
425450
426451 gmx angle -n angle.ndx -f production.xtc -od angle-distribution.xvg -binwidth 0.25
427452
428- Select 1 for the O-C-C-O dihedral.
453+ When prompted, select group 1 corresponding to the O-C-C-O dihedral.
429454
430455.. figure :: ../figures/level2/stretching-a-polymer/dihedral-distribution-light.png
431456 :alt: Angular distribution from molecular dynamics simulation in GROMACS
0 commit comments