diff --git a/files/tutorial3/PEG-GROMOS.mol b/files/tutorial3/mergePEGH2O/PEG-GROMOS.mol
similarity index 100%
rename from files/tutorial3/PEG-GROMOS.mol
rename to files/tutorial3/mergePEGH2O/PEG-GROMOS.mol
diff --git a/files/tutorial3/mergePEGH2O/input.lmp b/files/tutorial3/mergePEGH2O/input.lmp
new file mode 100644
index 0000000..48f2c0a
--- /dev/null
+++ b/files/tutorial3/mergePEGH2O/input.lmp
@@ -0,0 +1,45 @@
+
+units real
+atom_style full
+bond_style harmonic
+angle_style harmonic
+dihedral_style harmonic
+pair_style lj/cut/coul/long 10
+kspace_style pppm 1e-5
+special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 1.0 &
+angle yes dihedral yes
+
+read_data ../pureH2O/H2O.data &
+    extra/bond/per/atom 3 &
+    extra/angle/per/atom 6 &
+    extra/dihedral/per/atom 10 &
+    extra/special/per/atom 14
+include ../PARM.lmp
+
+molecule pegmol PEG-GROMOS.mol
+create_atoms 0 single 0 0 0 mol pegmol 454756
+
+group H2O type OW HW
+group PEG type C CPos H HC OAlc OE
+
+delete_atoms overlap 2.0 H2O PEG mol yes
+
+fix mynpt all npt temp 300 300 100 x 1 1 1000
+timestep 1.0
+
+dump mydmp all image 1000 dump.*.ppn type &
+ type shiny 0.1 box no 0.01 &
+ view 0 90 zoom 1.8  fsaa yes bond atom 0.8 
+dump_modify mydmp backcolor white &
+ acolor OW red acolor HW white &
+ acolor OE red acolor OAlc red &
+ acolor C gray acolor CPos gray &
+ acolor H white acolor HC white &
+ adiam OW 0.2 adiam HW 0.2 &
+ adiam C 3 adiam CPos 3 adiam OAlc 2.8 &
+ adiam H 1 adiam HC 1 adiam OE 2.8
+
+thermo 1000
+
+run 30000
+write_data mix.data
diff --git a/files/tutorial3/pullonPEG/input.lmp b/files/tutorial3/pullonPEG/input.lmp
new file mode 100644
index 0000000..6b64b1c
--- /dev/null
+++ b/files/tutorial3/pullonPEG/input.lmp
@@ -0,0 +1,59 @@
+variable f0 equal 5
+
+units real
+atom_style full
+bond_style harmonic
+angle_style harmonic
+dihedral_style harmonic
+pair_style lj/cut/coul/long 10
+kspace_style pppm 1e-5
+special_bonds lj 0.0 0.0 0.5 &
+coul 0.0 0.0 1.0 &
+angle yes dihedral yes
+
+read_data ../mergePEGH2O/mix.data
+include ../PARM.lmp
+
+group H2O type OW HW
+group PEG type C CPos H HC OAlc OE
+group ends type OAlc
+variable xcm equal xcm(ends,x)
+variable oxies atom type==label2type(atom,OAlc)
+variable end1 atom v_oxies*(x>v_xcm)
+variable end2 atom v_oxies*(x<v_xcm)
+
+group topull1 variable end1
+group topull2 variable end2
+
+dump mydmp all image 1000 dump.*.ppn type &
+type shiny 0.1 box no 0.01 &
+view 0 90 zoom 1.8 fsaa yes bond atom 0.8
+dump_modify mydmp backcolor white &
+acolor OW red acolor HW white &
+acolor OE red acolor OAlc red &
+acolor C gray acolor CPos gray &
+acolor H white acolor HC white &
+adiam OW 0.2 adiam HW 0.2 &
+adiam C 3 adiam CPos 3 adiam OAlc 2.8 &
+adiam H 1 adiam HC 1 adiam OE 2.8
+
+variable x1 equal xcm(topull1,x)
+variable x2 equal xcm(topull2,x)
+variable y1 equal xcm(topull1,y)
+variable y2 equal xcm(topull2,y)
+variable z1 equal xcm(topull1,z)
+variable z2 equal xcm(topull2,z)
+variable delta_r equal sqrt((v_x1-v_x2)^2+(v_y1-v_y2)^2+(v_z1-v_z2)^2) 
+compute rgyr PEG gyration
+thermo_style custom step temp etotal v_delta_r c_rgyr
+thermo 1000
+
+timestep 1.0
+fix mynvt all nvt temp 300 300 100
+
+run 30000
+
+fix myaf1 topull1 addforce ${f0} 0 0
+fix myaf2 topull2 addforce -${f0} 0 0
+
+run 30000
diff --git a/files/tutorial3/H2O-SPCFw.mol b/files/tutorial3/pureH2O/H2O-SPCFw.mol
similarity index 100%
rename from files/tutorial3/H2O-SPCFw.mol
rename to files/tutorial3/pureH2O/H2O-SPCFw.mol
diff --git a/files/tutorial3/pureH2O/input.lmp b/files/tutorial3/pureH2O/input.lmp
new file mode 100644
index 0000000..034181e
--- /dev/null
+++ b/files/tutorial3/pureH2O/input.lmp
@@ -0,0 +1,49 @@
+
+units real
+atom_style full
+bond_style harmonic
+angle_style harmonic
+dihedral_style harmonic
+pair_style lj/cut/coul/long 10
+kspace_style pppm 1e-5
+special_bonds lj 0.0 0.0 0.5 &
+coul 0.0 0.0 1.0 angle yes
+
+region box block -45 45 -15 15 -15 15
+  create_box 8 box &
+  bond/types 7 &
+  angle/types 8 &
+  dihedral/types 4 &
+  extra/bond/per/atom 3 &
+  extra/angle/per/atom 6 &
+  extra/dihedral/per/atom 10 &
+  extra/special/per/atom 14
+
+include ../PARM.lmp
+
+molecule h2omol H2O-SPCFw.mol
+create_atoms 0 random 1050 87910 NULL mol &
+h2omol 454756 overlap 1.0 maxtry 50
+
+group H2O type OW HW
+minimize 1.0e-4 1.0e-6 100 1000
+reset_timestep 0
+
+fix mynpt all npt temp 300 300 100 &
+  iso 1 1 1000
+
+dump mydmp all image 500 dump.*.jpg type type &
+  shiny 0.1 box no 0.01 view 0 90 zoom 1.8
+dump_modify mydmp backcolor white &
+  acolor OW red acolor HW white &
+  adiam OW 3 adiam HW 1.5
+
+variable myvol equal vol
+variable myoxy equal count(H2O)/3
+variable rho equal ${myoxy}/v_myvol
+thermo 1000
+thermo_style custom step temp etotal v_myvol v_rho
+
+timestep 1.0
+run 20000
+write_data H2O.data
diff --git a/lammps-tutorials.tex b/lammps-tutorials.tex
index 3fd2cdb..068cfb9 100644
--- a/lammps-tutorials.tex
+++ b/lammps-tutorials.tex
@@ -1856,7 +1856,7 @@ \subsubsection{Preparing the water reservoir}
 \end{verbatim}}
 
 Download and save the
-\href{https://raw.githubusercontent.com/lammpstutorials/lammpstutorials-article/main/files/tutorial3/PARM.lmp}{\textit{PARM.lmp}}
+\href{\filepath tutorial3/PARM.lmp}{\textit{PARM.lmp}}
 file next to the \textit{pureH2O/} folder. This example uses type labels, which
 are strings that map to each of the numeric atom types, bond types, angles types,
 etc. in the simulation. Type labels can be any string defined by the user, and we
@@ -1888,7 +1888,7 @@ \subsubsection{Preparing the water reservoir}
 \end{verbatim}}
 When LAMMPS fails to create the desired number of molecules, a WARNING appears
 in the \textit{log} file. The molecule template called
-\href{https://raw.githubusercontent.com/lammpstutorials/lammpstutorials-article/main/files/tutorial3/H2O-SPCFw.mol}{\textit{H2O-SPCFw.mol}}
+\href{\filepath tutorial3/pureH2O/H2O-SPCFw.mol}{\textit{H2O-SPCFw.mol}}
 must be downloaded and saved in the \textit{pureH2O/} folder. This template contains
 the necessary structural information of a water molecule, such as the number of
 atoms, or the IDs of the atoms that are connected by bonds, angles, etc.
@@ -1921,8 +1921,9 @@ \subsubsection{Preparing the water reservoir}
 
 Let us output the system into images by adding the following commands to \textit{input.lmp}:
 {\normalsize \begin{verbatim}
- dump mydmp all image 1000 dump.*.ppn type type &
- shiny 0.1 box no 0.01 view 0 90 zoom 1.8
+ dump mydmp all image 1000 dump.*.ppn type &
+     type shiny 0.1 box no 0.01 &
+     view 0 90 zoom 1.8
  dump_modify mydmp backcolor white &
    acolor OW red acolor HW white &
    adiam OW 3 adiam HW 1.5
@@ -1933,7 +1934,8 @@ \subsubsection{Preparing the water reservoir}
  variable myoxy equal count(H2O)/3
  variable myrho equal ${myoxy}/v_myvol
  thermo 1000
- thermo_style custom step temp etotal v_myvol v_myrho
+ thermo_style custom step temp etotal &
+     v_myvol v_myrho
 \end{verbatim}}
 The variable \textit{myoxy} corresponds to the number of atoms divided by 3, i.e.
 the number of molecules, and the variable \textit{myrho} is the number of molecules
@@ -1964,8 +1966,8 @@ \subsubsection{Preparing the water reservoir}
 \begin{figure}
 \centering
 \includegraphics[width=\linewidth]{PEG-density}
-\caption{Evolution of the density $\rho$ of water with time $t$. The density $\rho$
-reaches a plateau after $t \approx 10\,\text{ps}$.}
+\caption{Evolution of the density $\rho$ of water with time $t$. The density
+quickly reaches a plateau.}
 \label{fig:PEG-density}
 \end{figure}
 
@@ -1988,81 +1990,80 @@ \subsubsection{Solvating the PEG in water}
 Create a new blank file in it, call it \textit{input.lmp}. Within \textit{input.lmp},
 copy the same first lines as previously:
 {\normalsize \begin{verbatim}
-units real
-atom_style full
-bond_style harmonic
-angle_style harmonic
-dihedral_style harmonic
-pair_style lj/cut/coul/long 10
-kspace_style pppm 1e-5
-special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 1.0 &
-    angle yes dihedral yes
+ units real
+ atom_style full
+ bond_style harmonic
+ angle_style harmonic
+ dihedral_style harmonic
+ pair_style lj/cut/coul/long 10
+ kspace_style pppm 1e-5
+ special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 1.0 &
+     angle yes dihedral yes
 \end{verbatim}}
 Then, import the previously generated data file \textit{H2O.data} as well as the \textit{PARM.lmp} file:
 {\normalsize \begin{verbatim}
-read_data ../pureH2O/H2O.data &
-    extra/bond/per/atom 3 &
-    extra/angle/per/atom 6 &
-    extra/dihedral/per/atom 10 &
-    extra/special/per/atom 14
-include ../PARM.lmp
+ read_data ../pureH2O/H2O.data &
+     extra/bond/per/atom 3 &
+     extra/angle/per/atom 6 &
+     extra/dihedral/per/atom 10 &
+     extra/special/per/atom 14
+ include ../PARM.lmp
 \end{verbatim}}
 Download the template called
-\href{https://raw.githubusercontent.com/lammpstutorials/lammpstutorials-article/main/files/tutorial3/PEG-GROMOS.mol}{\textit{PEG-GROMOS.mol}}
+\href{\filepath tutorial3/mergePEGH2O/PEG-GROMOS.mol}{\textit{PEG-GROMOS.mol}}
 for the PEG molecule, and then create a single molecule in the middle of the box:
 {\normalsize \begin{verbatim}
-molecule pegmol PEG-GROMOS.mol
-create_atoms 0 single 0 0 0 mol pegmol 454756
+ molecule pegmol PEG-GROMOS.mol
+ create_atoms 0 single 0 0 0 mol pegmol 454756
 \end{verbatim}}
 Let us create 2 groups to differentiate the PEG from the H2O, by adding the following
 lines into \textit{input.lmp}:
 {\normalsize \begin{verbatim}
-group H2O type OW HW
-group PEG type C CPos H HC OAlc OE
+ group H2O type OW HW
+ group PEG type C CPos H HC OAlc OE
 \end{verbatim}}
 Water molecules that are overlapping with the PEG must be deleted to avoid future crashing.
 Add the following line into \textit{input.lmp}:
 {\normalsize \begin{verbatim}
-delete_atoms overlap 2.0 H2O PEG mol yes
+ delete_atoms overlap 2.0 H2O PEG mol yes
 \end{verbatim}}
 Here the value of 2 Ångstroms for the overlap cutoff was fixed arbitrarily and can
 be chosen through trial and error. If the cutoff is too small, the simulation will
 crash. If the cutoff is too large, too many water molecules will unnecessarily be
-deleted. Finally, let us use the \textit{fix npt} to control the temperature, as
+deleted. 
+
+Finally, let us use the \textit{fix npt} to control the temperature, as
 well as the pressure by allowing the box size to be rescaled along the \textit{x} axis:
 {\normalsize \begin{verbatim}
-fix mynpt all npt temp 300 300 100 x 1 1 1000
-timestep 1.0
-\end{verbatim}}
-Once more, let us dump the atom positions as well as the system temperature and volume:
-{\normalsize \begin{verbatim}
-dump mydmp all atom 100 dump.lammpstrj
-thermo 100
-variable mytemp equal temp
-variable myvol equal vol
-fix myat1 all ave/time 10 10 100 &
-    v_mytemp file temperature.dat
-fix myat2 all ave/time 10 10 100 &
-    v_myvol file volume.dat
+ fix mynpt all npt temp 300 300 100 x 1 1 1000
+ timestep 1.0
 \end{verbatim}}
-Let us also print the total enthalpy:
+Once more, let us create images of the systems:
 {\normalsize \begin{verbatim}
-variable myenthalpy equal enthalpy
-fix myat3 all ave/time 10 10 100 &
-    v_myenthalpy file enthalpy.dat
+ dump mydmp all image 1000 dump.*.ppn type &
+     type shiny 0.1 box no 0.01 &
+     view 0 90 zoom 1.8  fsaa yes bond atom 0.8 
+ dump_modify mydmp backcolor white &
+     acolor OW red acolor HW white &
+     acolor OE red acolor OAlc red &
+     acolor C gray acolor CPos gray &
+     acolor H white acolor HC white &
+     adiam OW 0.2 adiam HW 0.2 &
+     adiam C 3 adiam CPos 3 adiam OAlc 2.8 &
+     adiam H 1 adiam HC 1 adiam OE 2.8
+ thermo 1000
 \end{verbatim}}
 Finally, let us perform a short equilibration and print the
 final state in a data file. Add the following lines to the data file:
 {\normalsize \begin{verbatim}
-run 30000
-write_data mix.data
+ run 30000
+ write_data mix.data
 \end{verbatim}}
-If you open the \textit{dump.lammpstrj} file using VMD or have a look at the
-evolution of the volume in \textit{volume.dat}, you should see that the box
-dimensions slightly evolve along \textit{x} to accommodate the new configuration.
-In addition, the temperature remains close to the target value of $300~\text{K}$
-throughout the entire simulation, and the enthalpy is almost constant, suggesting
-that the system was close to equilibrium from the start.
+From the outputs, you can make sure that the temperature remains close to the
+target value of $300~\text{K}$ throughout the entire simulation, and the that 
+the volume and total energy are almost constant, suggesting
+that the system was close to equilibrium from the start. See a snapshot of the
+system in Fig.\,\ref{fig:PEG-solvated}.
 
 \begin{figure}
 \centering
@@ -2086,22 +2087,22 @@ \subsubsection{Stretching the PEG molecule}
 contribution from both water and PEG molecules (it was chosen by trial and error).
 Then, copy the same lines as previously:
 {\normalsize \begin{verbatim}
-units real
-atom_style full
-bond_style harmonic
-angle_style harmonic
-dihedral_style harmonic
-pair_style lj/cut/coul/long 10
-kspace_style pppm 1e-5
-special_bonds lj 0.0 0.0 0.5 &
-    coul 0.0 0.0 1.0 &
-    angle yes dihedral yes
+ units real
+ atom_style full
+ bond_style harmonic
+ angle_style harmonic
+ dihedral_style harmonic
+ pair_style lj/cut/coul/long 10
+ kspace_style pppm 1e-5
+ special_bonds lj 0.0 0.0 0.5 &
+     coul 0.0 0.0 1.0 &
+     angle yes dihedral yes
 \end{verbatim}}
 Start the simulation from the equilibrated PEG-water system and include again the
 parameter file by adding the following lines to the \textit{input.lmp}:
 {\normalsize \begin{verbatim}
-read_data ../mergePEGH2O/mix.data
-include ../PARM.lmp
+ read_data ../mergePEGH2O/mix.data
+ include ../PARM.lmp
 \end{verbatim}}
 Next, let us create some useful atom groups, including H2O and PEG (as previously),
 as well as 2 groups containing a single atom each, the oxygen atoms at the ends
@@ -2109,20 +2110,29 @@ \subsubsection{Stretching the PEG molecule}
 atoms located at the ends of the PEG molecule, which we are going to use to pull
 on the PEG molecule. Add the following lines to the \textit{input.lmp}:
 {\normalsize \begin{verbatim}
-group H2O type OW HW
-group PEG type C CPos H HC OAlc OE
-group ends type OAlc
-variable xcm equal xcm(ends,x)
-variable oxies atom type==label2type(atom,OAlc)
-variable end1 atom v_oxies*(x>v_xcm)
-variable end2 atom v_oxies*(x<v_xcm)
-group topull1 variable end1
-group topull2 variable end2
-\end{verbatim}}
-Add the following \textit{dump} command to the input to print the atom positions
-every 1000 steps:
-{\normalsize \begin{verbatim}
-dump mydmp all atom 1000 dump.lammpstrj
+ group H2O type OW HW
+ group PEG type C CPos H HC OAlc OE
+ group ends type OAlc
+ variable xcm equal xcm(ends,x)
+ variable oxies atom type==label2type(atom,OAlc)
+ variable end1 atom v_oxies*(x>v_xcm)
+ variable end2 atom v_oxies*(x<v_xcm)
+ group topull1 variable end1
+ group topull2 variable end2
+\end{verbatim}}
+Add the following \textit{dump} command to create images of the system:
+{\normalsize \begin{verbatim}
+ dump mydmp all image 1000 dump.*.ppn type &
+     type shiny 0.1 box no 0.01 &
+     view 0 90 zoom 1.8 fsaa yes bond atom 0.8
+ dump_modify mydmp backcolor white &
+     acolor OW red acolor HW white &
+     acolor OE red acolor OAlc red &
+     acolor C gray acolor CPos gray &
+     acolor H white acolor HC white &
+     adiam OW 0.2 adiam HW 0.2 &
+     adiam C 3 adiam CPos 3 adiam OAlc 2.8 &
+     adiam H 1 adiam HC 1 adiam OE 2.8
 \end{verbatim}}
 Let us use a single Nosé-Hoover thermostat applied to all the atoms by adding the
 following lines to \textit{input.lmp}:
@@ -2132,44 +2142,37 @@ \subsubsection{Stretching the PEG molecule}
 \end{verbatim}}
 Let us also print the end-to-end distance of the PEG,
 here defined as the distance between the groups \textit{topull1}
-and \textit{topull2}, as well as the temperature of the system and the gyration
+and \textit{topull2}, as well as the gyration
 radius of the molecule \cite{fixmanRadiusGyrationPolymer1962a}
 by adding the following lines to \textit{input.lmp}:
-% S.G.: should we keep the data file, or just use thermo custom as its better for LAMMPS-GUI?
-{\normalsize \begin{verbatim}
-variable mytemp equal temp
-fix myat1 all ave/time 10 10 100 &
-    v_mytemp file temperature.dat
-variable x1 equal xcm(topull1,x)
-variable x2 equal xcm(topull2,x)
-variable y1 equal xcm(topull1,y)
+{\normalsize \begin{verbatim}
+ variable x1 equal xcm(topull1,x)
+ variable x2 equal xcm(topull2,x)
+ variable y1 equal xcm(topull1,y)
 variable y2 equal xcm(topull2,y)
-variable z1 equal xcm(topull1,z)
-variable z2 equal xcm(topull2,z)
-variable delta_r equal &
-    sqrt((v_x1-v_x2)^2+(v_y1-v_y2)^2 &
-    +(v_z1-v_z2)^2)
-fix myat2 all ave/time 10 10 100 v_delta_r &
-    file end-to-end-distance.dat
-compute rgyr PEG gyration
-fix myat3 all ave/time 10 10 100 c_rgyr &
-    file gyration-radius.dat
-thermo 1000
+ variable z1 equal xcm(topull1,z)
+ variable z2 equal xcm(topull2,z)
+ variable delta_r equal &
+     sqrt((v_x1-v_x2)^2+(v_y1-v_y2)^2 &
+     +(v_z1-v_z2)^2) 
+ compute rgyr PEG gyration
+ thermo_style custom step temp etotal v_delta_r c_rgyr
+ thermo 1000
 \end{verbatim}}
 Finally, let us simulate 30 picoseconds without any external forcing:
 {\normalsize \begin{verbatim}
-run 30000
+ run 30000
 \end{verbatim}}
 This first run will serve as a benchmark to quantify later the changes induced
-by the forcing. Then, let us apply a forcing on the 2 oxygen atoms using two
+by the forcing. Then, let us apply a forcing on the 2 selected oxygen atoms using two
 \textit{add\_force} commands, and run for an extra 30 ps:
 {\normalsize \begin{verbatim}
-fix myaf1 topull1 addforce ${f0} 0 0
-fix myaf2 topull2 addforce -${f0} 0 0
-run 30000
+ fix myaf1 topull1 addforce ${f0} 0 0
+ fix myaf2 topull2 addforce -${f0} 0 0
+ run 30000
 \end{verbatim}}
-Run the \textit{input.lmp} file using LAMMPS. If you open the \textit{dump.lammpstrj}
-file using \textit{VMD}, you should see that the PEG molecule eventually aligns
+Run the \textit{input.lmp} file using LAMMPS. From the generated images of the system,
+you should see that the PEG molecule eventually aligns
 in the direction of the force (as seen in Fig.\,\ref{fig:PEG-in-water}).
 The evolutions of the end-to-end distance and of the gyration radius over
 time show that the PEG is adjusting to the external forcing in less than