@@ -2468,11 +2468,11 @@ \subsubsection{System preparation}
24682468\end {lstlisting }
24692469These lines are used to define the most basic parameters, including the
24702470{\color {blue}atom style, the form of the non-bonded, bond and angle potentials as well as
2471- other specifics about the non-boned interactions}. Here, \lmpcmd {lj/cut/tip4p/long}
2472- imposes a Lennard-Jones potential with a cut-off at $ 12 \, \text {\AA {}}$ and a long-range
2473- Coulomb potential. {\color {blue} The parameters \lmpcmd {O}, \lmpcmd {H}, \lmpcmd {O-H},
2474- and \lmpcmd {H-O-H} correspond respectively to the oxygens, hydrogens, O-H bonds, and
2475- H-O-H angle constraints of the water molecules; their definitions, provided by the
2471+ other specifics about the non-boned interactions}. Here, \lmpcmd {lj/cut/tip4p/long}
2472+ imposes a Lennard-Jones potential with a cut-off at $ 12 \, \text {\AA {}}$ and a long-range
2473+ Coulomb potential. {\color {blue} The parameters \lmpcmd {O}, \lmpcmd {H}, \lmpcmd {O-H},
2474+ and \lmpcmd {H-O-H} correspond respectively to the oxygens, hydrogens, O-H bonds, and
2475+ H-O-H angle constraints of the water molecules; their definitions, provided by the
24762476\lmpcmd {labelmap} commands, will be clarified below.}
24772477
24782478So far, the commands are relatively similar to those in the previous tutorial,
@@ -2618,15 +2618,17 @@ \subsubsection{System preparation}
26182618bond_coeff O-H 0 0.9572
26192619angle_coeff H-O-H 0 104.52
26202620\end {lstlisting }
2621- The \lmpcmd {bond\_ coeff} command, used here for the O-H bond of the water
2622- molecule, sets both the spring constant of the harmonic potential and the
2623- equilibrium bond distance of $ 0.9572 ~\text {\AA {}}$ . The constant can be 0 for a
2624- rigid water molecule because the SHAKE algorithm, {\color {blue} which will be summoned
2625- later in the script using a command line,} will {\color {blue} fix the intramolecular
2626- structure} of the water {\color {blue} molecules} (see below)~\cite {ryckaert1977numerical , andersen1983rattle }.
2627- Similarly, the \lmpcmd {angle\_ coeff} command for the H-O-H angle of the water molecule sets
2628- the force constant of the angular harmonic potential to 0 and the equilibrium
2629- angle to $ 104.52 ^\circ $ .
2621+ The \lmpcmd {bond\_ coeff} command, used here for the O-H bond of the
2622+ water molecule, sets both the spring constant of the harmonic potential
2623+ and the equilibrium bond distance of $ 0.9572 ~\text {\AA {}}$ . The force
2624+ constant can be 0 for a rigid water molecule because the SHAKE
2625+ algorithm, {\color {blue} which will be used in the input at a later
2626+ step,} will {\color {blue} constrain the intramolecular structure} of
2627+ the water {\color {blue} molecules} (see
2628+ below)~\cite {ryckaert1977numerical , andersen1983rattle }. Similarly, the
2629+ \lmpcmd {angle\_ coeff} command for the H-O-H angle of the water molecule
2630+ sets the force constant of the angular harmonic potential to 0 and the
2631+ equilibrium angle to $ 104.52 ^\circ $ .
26302632
26312633Alongside \flecmd {parameters.inc}, the \flecmd {groups.inc} file contains
26322634several \lmpcmd {group} commands to \textcolor {blue}{define groups of atoms based
@@ -3105,7 +3107,7 @@ \subsubsection{Prepare and relax}
31053107and a \lmpcmd {.data} file is imported by the \lmpcmd {read\_ data} command.
31063108
31073109The initial topology given by \href {\filepath tutorial5/silica.data}{\dwlcmd {silica.data}}
3108- {\color {blue} corresponds to} a small amorphous silica structure.
3110+ {\color {blue} corresponds to} a small amorphous silica structure.
31093111{\color {blue}This structure was generated in a prior
31103112simulation using the Vashishta force field~\cite {vashishta1990interaction }.}
31113113If you open the \flecmd {silica.data} file, you will find in the \lmpcmd {Atoms}
@@ -3125,7 +3127,7 @@ \subsubsection{Prepare and relax}
31253127pair_coeff * * ffield.reax.CHOFe Si O
31263128fix myqeq all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff maxiter 400
31273129\end {lstlisting }
3128- In this case, the \lmpcmd {pair\_ style reaxff} is used without a control file
3130+ In this case, the \lmpcmd {pair\_ style reaxff} is used without a control file
31293131{\color {blue}(see note below)}. The
31303132\lmpcmd {safezone} and \lmpcmd {mincap} keywords are added to prevent
31313133allocation issues, which sometimes can trigger segmentation faults and
@@ -3159,7 +3161,7 @@ \subsubsection{Prepare and relax}
31593161variable vq atom q
31603162\end {lstlisting }
31613163{\color {blue} The definition of the equal style variables qSi and qO
3162- make use of functions pre-defined within LAMMPS that allow calculating
3164+ make use of functions pre-defined within LAMMPS that allow calculating
31633165the total charge of atoms belonging to a group (charge()) and the total
31643166number of atoms in the group (count()). }
31653167To print the averaged charges \lmpcmd {qSi} and \lmpcmd {qO} using the
@@ -3445,7 +3447,7 @@ \subsubsection{Decorate the surface}
34453447fix myspec all reaxff/species 5 1 5 decorate.species &
34463448 element Si O H
34473449\end {lstlisting }
3448- {\color {blue} The commands above are, once again, similar to the ones of the previous script.}
3450+ {\color {blue} The commands above are, once again, similar to the ones of the previous script.}
34493451Here, the $ +1 \mathrm {e}{-10}$ was added to the denominator of the \lmpcmd {variable qH}
34503452to avoid dividing by 0 at the beginning of the simulation{\color {blue}, as no hydrogen
34513453atoms exists in the simulation domain yet}. Finally, let us
@@ -3529,7 +3531,7 @@ \subsubsection{Generation of the silica block}
35293531create_atoms Si random 240 5802 box overlap 2.0 maxtry 500
35303532create_atoms O random 480 1072 box overlap 2.0 maxtry 500
35313533\end {lstlisting }
3532- {\color {blue}In line with what is done in previous tutorials, the}
3534+ {\color {blue}In line with what is done in previous tutorials, the}
35333535\lmpcmd {create\_ atoms} commands are used to place
35343536240 Si atoms and 480 O atoms, respectively. This corresponds to
35353537an initial density of approximately $ 2 $ \, g/cm$ ^3 $ , which is close
@@ -3585,9 +3587,9 @@ \subsubsection{Generation of the silica block}
35853587fix mynvt all nvt temp 6000 300 0.1
35863588run 30000
35873589\end {lstlisting }
3588- {\color {blue} In this case, the initial and final target temperatures
3590+ {\color {blue} In this case, the initial and final target temperatures
35893591set for the Nose-Hoover thermostat is different, causing it to evolve
3590- linearly within the number of timesteps evoked in the \lmpcmd {run}
3592+ linearly within the number of timesteps evoked in the \lmpcmd {run}
35913593command. }
35923594In the third step, the system is equilibrated at the final desired
35933595conditions, $ T = 300 \, \text {K}$ and $ p = 1 \, \text {atm}$ ,
@@ -3601,7 +3603,7 @@ \subsubsection{Generation of the silica block}
36013603write_data generate.data
36023604\end {lstlisting }
36033605Here, an anisotropic barostat is used.
3604- {\color {blue}As previously mentioned, a}nisotropic
3606+ {\color {blue}As previously mentioned, a}nisotropic
36053607barostats adjust the dimensions independently, which is
36063608generally suitable for a solid phase.
36073609
@@ -3611,7 +3613,7 @@ \subsubsection{Generation of the silica block}
36113613{\color {blue}is deforming} during the last stage of the simulation
36123614(Fig.~\ref {fig:GCMC-dimension }\, b). After the simulation completes, the final
36133615{\color {blue} microstate attained during the dynamics and the system topology
3614- will be written to a} LAMMPS {\color {blue} data }file called \flecmd {generate.data}
3616+ will be written to a} LAMMPS {\color {blue} data }file called \flecmd {generate.data}
36153617{\color {blue}which} will be located next to \flecmd {generate.lmp} (Fig.~\ref {fig:GCMC-snapshot }).
36163618
36173619\begin {figure }
@@ -3647,7 +3649,6 @@ \subsubsection{Cracking the silica}
36473649thermo 250
36483650thermo_style custom step temp etotal vol density
36493651\end {lstlisting }
3650-
36513652Let us progressively increase the size of the box in the $ x$ direction,
36523653forcing the silica to deform and eventually crack. To achive this,
36533654the \lmpcmd {fix deform} command is used, with a rate
@@ -3662,8 +3663,7 @@ \subsubsection{Cracking the silica}
36623663write_data cracking.data
36633664\end {lstlisting }
36643665{\color {blue}As discussed, the \lmpcmd {fix nvt} command integrates the Nosé-Hoover equations
3665- of motion as originally derived to sample the NVT ensemble,
3666- which allows controlling the temperature of the system.}
3666+ of motion to sample the NVT ensemble, which allows controlling the temperature of the system.}
36673667As observed from the generated images, the atoms
36683668progressively adjust to the changing box dimensions. At some point,
36693669bonds begin to break, leading to the appearance of
@@ -3688,13 +3688,13 @@ \subsubsection{Adding water}
36883688
36893689To add the water molecules to the silica, we will employ the Monte Carlo
36903690method in the grand canonical ensemble (GCMC). In short, the system is
3691- placed into contact with a virtual reservoir {\color {blue} containing pure
3692- water at a given thermodynamic state}, and multiple attempts to insert
3693- water molecules at random positions are made. {\color {blue} In the grand
3694- canonical ensemble, each} attempt is either accepted or rejected based on
3695- {\color {blue} internal} energy {\color {blue} and chemical potential, $ \mu $ }
3691+ placed into contact with a virtual reservoir {\color {blue} containing pure
3692+ water at a given thermodynamic state}, and multiple attempts to insert
3693+ water molecules at random positions are made. {\color {blue} In the grand
3694+ canonical ensemble, each} attempt is either accepted or rejected based on
3695+ {\color {blue} internal} energy {\color {blue} and chemical potential, $ \mu $ }
36963696considerations. For further details, please refer to
3697- classical textbooks like Ref.~\citenum {frenkel2023understanding}.
3697+ classical textbooks like Ref.~\citenum {frenkel2023understanding}.
36983698
36993699% \paragraph{Using hydrid potentials}
37003700
@@ -3706,11 +3706,11 @@ \subsubsection{Adding water}
37063706% $\text{SiO}_2$), and TIP4P (for water). Here, the TIP4P/2005 model is
37073707% employed for the water~\cite{abascal2005general}. Open the
37083708% \flecmd{gcmc.lmp} file, which should contain the following lines:
3709- For this next step, we need to {\color {blue} specify the force field used
3709+ For this next step, we need to {\color {blue} specify the force field used
37103710to model the interactions in the system}. The TIP4P/2005 model is employed
37113711for the water~\cite {abascal2005general }, while no {\color {blue} interaction
3712- within silica is defined, as it will be seen farther below}.
3713- {\color {blue} This is because the} atoms of the silica will
3712+ within silica is defined, as it will be seen farther below}.
3713+ {\color {blue} This is because the} atoms of the silica will
37143714remain frozen during this part {\color {blue}of the simulation}.
37153715Only the cross-interactions between water and silica need
37163716to be defined. Open the \flecmd {gcmc.lmp} file, which should contain the following lines:
@@ -3731,7 +3731,7 @@ \subsubsection{Adding water}
37313731The PPPM solver~\cite {luty1996calculating } is specified with the \lmpcmd {kspace}
37323732command, and is used to compute the long-range Coulomb interactions associated
37333733with \lmpcmd {tip4p/long}. Finally, the {\color {blue} form of} the bond
3734- and angle{\color {blue} potentials} of the water molecules are defined; however,
3734+ and angle{\color {blue} potentials} of the water molecules are defined; however,
37353735{\color {blue} as previously discussed,} these specifications are
37363736not critical since TIP4P/2005 is a rigid water model.}
37373737
@@ -3812,7 +3812,7 @@ \subsubsection{Adding water}
38123812file, the \lmpcmd {create\_ atoms} command is used to include three water molecules
38133813in the system. Then, add the following \lmpcmd {pair\_ coeff} (and
38143814\lmpcmd {bond\_ coeff} and \lmpcmd {angle\_ coeff}) commands
3815- to \flecmd {gcmc.lmp} {\color {blue}in order to specify the potential parameters}:
3815+ to \flecmd {gcmc.lmp} {\color {blue}in order to set the potential parameters}:
38163816\ begin{lstlisting}
38173817pair_coeff * * 0 0
38183818pair_coeff Si OW 0.0057 4.42
@@ -3825,8 +3825,8 @@ \subsubsection{Adding water}
38253825% The force field Vashishta applies only to \lmpcmd{Si} and \lmpcmd{O} of $\text{SiO}_2$,
38263826% and not to the \lmpcmd{OW} and \lmpcmd{HW} of $\text{H}_2\text{O}$, thanks to the \lmpcmd{NULL} parameters
38273827% used for atoms of types \lmpcmd{OW} and \lmpcmd{HW}.
3828- Pair coefficients for the { \color {blue}potentials underlying the} \ lmpcmd {lj/cut/tip4p/long}
3829- { \color {blue}style} are defined between O($ \text {H}_2 \text {O}$ ) and between H($ \text {H}_2 \text {O}$ )
3828+ Pair coefficients for the \ lmpcmd {lj/cut/tip4p/long} { \color {blue}pair style} are
3829+ defined between O($ \text {H}_2 \text {O}$ ) and between H($ \text {H}_2 \text {O}$ )
38303830atoms, as well as between O($ \text {SiO}_2 $ )-O($ \text {H}_2 \text {O}$ ) and
38313831Si($ \text {SiO}_2 $ )-O($ \text {H}_2 \text {O}$ ). Thus, the fluid-fluid and the
38323832fluid-solid interactions will be adressed with by the \lmpcmd {lj/cut/tip4p/long} potential.
@@ -4116,7 +4116,7 @@ \subsubsection{Method 1: Free sampling}
41164116fix myadf all addforce v_F 0.0 0.0 energy v_U
41174117\end {lstlisting }
41184118Next, we {\color {blue}use the Newtonian equations of motion with a
4119- Langevin thermostat by combining }the \lmpcmd {fix nve} with a
4119+ Langevin thermostat by combining }the \lmpcmd {fix nve} with a
41204120\lmpcmd {fix langevin} {\color {blue}command}:
41214121\ begin{lstlisting}
41224122fix mynve all nve
@@ -4128,7 +4128,7 @@ \subsubsection{Method 1: Free sampling}
41284128fluctuates around a target value.
41294129% SG: may be discuss the choice of constant "500" -> chosen for a relatiely weak coupling with thermostat (add a box?)
41304130% CA: I can propose something that can be put in a yellow box. Feel free to uncomment the lines below and modify it if you want to.
4131- % {\color{blue}You can find proposed in the LAMMPS documentation
4131+ % {\color{blue}You can find proposed in the LAMMPS documentation
41324132% values of 100x and 1000x the value of the timestep for the damping
41334133% constants of the thermostats and barostats. Here, a smaller
41344134% value is used in order to have the temperature of the system
@@ -4147,7 +4147,7 @@ \subsubsection{Method 1: Free sampling}
41474147
41484148To ensure that the equilibration time is sufficient, we will track the evolution of
41494149the number of atoms in the central - energetically unfavorable - region,
4150- {\color {blue} defined below under the name} \lmpcmd {mymes}, using
4150+ {\color {blue} defined below under the name} \lmpcmd {mymes}, using
41514151the \lmpcmd {n\_ center} variable:
41524152\ begin{lstlisting}
41534153region mymes block -${x0} ${x0} INF INF INF INF
@@ -4161,7 +4161,7 @@ \subsubsection{Method 1: Free sampling}
41614161dump_modify viz backcolor white acolor 1 cyan &
41624162 adiam 1 3 boxcolor black
41634163\end {lstlisting }
4164- A \lmpcmd {dump image} command was also added for system visualization.
4164+ A \lmpcmd {dump image} command was also added for system visualization.
41654165{\color {blue} The other commands should also be familiar from previous tutorials.}
41664166
41674167Finally, let us perform an equilibration of 50000 steps,
@@ -4328,7 +4328,7 @@ \subsubsection{Method 2: Umbrella sampling}
43284328\end {lstlisting }
43294329
43304330So far, our code resembles that of Method 1, except for the additional particle
4331- of type 2. Particles of types 1 and 2 are identical.
4331+ of type 2. Particles of types 1 and 2 are identical.
43324332However, the particle of type 2 will also
43334333be exposed to the biasing potential $ V$ , which forces it to explore the
43344334central part of the box (Fig.~\ref {fig:US-system-biased }), {\color {blue} thus
@@ -4366,9 +4366,9 @@ \subsubsection{Method 2: Umbrella sampling}
43664366jump SELF loop
43674367\end {lstlisting }
43684368{\color {blue} The definition of a variable of loop style serves the same purpose
4369- as in (\hyperref [reactive-silicon-dioxide-label]{Tutorial 5}), and we highlight
4369+ as in (\hyperref [reactive-silicon-dioxide-label]{Tutorial 5}), and we highlight
43704370here the particular utility of using its value
4371- to distinguish the files written by the \lmpcmd {fix ave_/time } command for the
4371+ to distinguish the files written by the \lmpcmd {fix ave \_ time } command for the
43724372different bias potentials.}
43734373The \lmpcmd {spring} command imposes the additional harmonic potential $ V$ with
43744374the previously defined spring constant $ k$ {\color {blue}to the atoms in the group \lmpcmd {topull}}. The center of the harmonic
@@ -4476,8 +4476,8 @@ \subsubsection{Creating the system}
44764476\end {lstlisting }
44774477The \lmpcmd {class2} styles compute a 6/9 Lennard-Jones potential~\cite {sun1998compass }.
44784478The \textit {class2 } bond, angle, dihedral, and improper styles are used as
4479- well, see the documentation for a description of the respective {\color {blue} potential
4480- form they, each, prescribe}.
4479+ well, see the documentation for a description of the respective {\color {blue} potential
4480+ form they, each, prescribe}.
44814481{\color {blue}The \lmpcmd {tail yes} option adds long-range van der Waals tail
44824482corrections to the energy and pressure.}
44834483The \lmpcmd {mix sixthpower} imposes the following mixing rule for the calculation
@@ -4510,10 +4510,10 @@ \subsubsection{Creating the system}
45104510minimize 1.0e-4 1.0e-6 100 1000
45114511reset_timestep 0
45124512\end {lstlisting }
4513- {\color {blue}These commands should be familiar from previous tutorials.}
4513+ {\color {blue}These commands should be familiar from previous tutorials.}
45144514
45154515Then, let us densify the system to a target value of $ 0.9 ~\text {g/cm}^3 $
4516- by {\color {blue}imposing the} shrinking {\color {blue}of} the simulation box at a constant rate.
4516+ by {\color {blue}imposing the} shrinking {\color {blue}of} the simulation box at a constant rate.
45174517The dimension parallel to the CNT axis is maintained fixed because the CNT is periodic in that direction.
45184518Add the following commands to \flecmd {mixing.lmp}:
45194519% SG: I removed the loop local, unless its important? But if it is, we have to
@@ -4565,7 +4565,7 @@ \subsubsection{Creating the system}
45654565
45664566write_data mixing.data
45674567\end {lstlisting }
4568- For visualization purposes, the atoms {\color {blue} of}
4568+ For visualization purposes, the atoms {\color {blue} of}
45694569the CNT \lmpcmd {group} {\color {blue}are} moved
45704570to the center of the box using \lmpcmd {fix recenter}.
45714571As the time progresses, the system density,
@@ -4704,25 +4704,18 @@ \subsubsection{Simulating the reaction}
47044704% as I literally just figured it out to write them. (PS: I just realized that some
47054705% of it may be redundant with the note, but I would still leave it as the reader
47064706% will have the full information.
4707- With the \lmpcmd {stabilization} keyword, the \lmpcmd {bond/react} command will
4707+ With the \lmpcmd {stabilization} keyword, the \lmpcmd {fix bond/react} command will
47084708stabilize the atoms involved in the reaction using the \lmpcmd {nve/limit}
4709- command with a maximum displacement of $ 0.03 \, \mathrm {\AA {}}$ .
4710- {\color {blue}The \lmpcmd {nve/limit} is a command that functions similar to the
4711- \lmpcmd {nve} in terms of equation of motion implied, but has a constraint on
4712- how close atoms will be. In this case, it will be summoned internally within
4713- LAMMS as the \lmpcmd {bond/react} appears.}
4709+ command with a maximum displacement of $ 0.03 \, \mathrm {\AA {}}$ .
4710+ {\color {blue}The \lmpcmd {fix nve/limit} command functions similar to
4711+ \lmpcmd {fix nve}, but restricts how far atoms can move in a single time step, even with
4712+ very large forces.}
47144713By default, each reaction is stabilized for 60 time steps. Each \lmpcmd {react} keyword
47154714corresponds to a reaction, e.g.,~a transformation of \lmpcmd {mol1} into \lmpcmd {mol2}
47164715based on the atom map \lmpcmd {M-M.rxnmap}. Implementation details about each reaction,
47174716such as the reaction distance cutoffs and the frequency with which to search for
47184717reaction sites, are also specified in this command.
47194718
4720- \begin {note }
4721- {\color {blue}The \lmpcmd {fix nve/limit} integrates command Newton's equations of motion
4722- while limiting the maximum displacement of atoms per timestep. This is
4723- useful for preventing atoms from moving too far due to large forces.}
4724- \end {note }
4725-
47264719\begin {figure }
47274720\centering
47284721\includegraphics [width=\linewidth ]{REACT-final.png}
@@ -4741,7 +4734,7 @@ \subsubsection{Simulating the reaction}
47414734undergoing stabilization is named by appending `\_ REACT' to the user-provided prefix.
47424735\end {note }
47434736
4744- Add the following commands to \flecmd {polymerize.lmp} to
4737+ Add the following commands to \flecmd {polymerize.lmp} to
47454738{\color {blue} carry out the dynamics using a Nosé-Hoover thermostat}
47464739while ensuring that the CNT remains centered in the simulation box:
47474740\ begin{lstlisting}
0 commit comments