.. _stretching-polymer-label: Stretching a polymer ******************** .. container:: hatnote Solvating and stretching a small polymer molecule .. figure:: ../figures/level2/stretching-a-polymer/video-PEG-dark.webp :alt: Movie of a peg polymer molecule in water as simulated with GROMACS :height: 250 :align: right :class: only-dark .. figure:: ../figures/level2/stretching-a-polymer/video-PEG-light.webp :alt: Movie of a peg polymer molecule in water as simulated with GROMACS :height: 250 :align: right :class: only-light .. container:: justify The goal of this tutorial is to use GROMACS and solvate a small hydrophilic polymer in a reservoir of water. .. container:: justify An all-atom description is used for both polymer and water. The polymer is PolyEthylene Glycol (PEG). Once the system is properly equilibrated at the desired temperature and pressure, a force is applied to both ends of the polymer. The evolution of the polymer length is measured, and the energetics of the system is measured. .. (GROMOS 54A7 force field :cite:`schmid2011definition`) (SPC flexible model :cite:`wu2006flexible`) .. container:: justify This tutorial was inspired by a |Liese2017| by Liese and coworkers, in which molecular dynamics simulations are compared with force spectroscopy experiments :cite:`liese2017hydration`. .. |Liese2017| raw:: html publication .. include:: ../../non-tutorials/recommand-salt.rst .. include:: ../../non-tutorials/needhelp.rst .. include:: ../../non-tutorials/GROMACS2024.2.rst Prepare the PEG molecule ======================== .. container:: justify Download the *peg.gro* file for the PEG molecule by clicking |download_H2O.data|. The *peg.gro* file can be visualized using vmd, by typing in a terminal: .. |download_H2O.data| raw:: html here .. code-block:: bash vmd peg.gro .. figure:: ../figures/level2/stretching-a-polymer/light-PEG.png :alt: PEG polymer for molecular dynamics simulation in GROMACS :class: only-light :width: 500 .. figure:: ../figures/level2/stretching-a-polymer/dark-PEG.png :alt: PEG polymer for molecular dynamics simulation in GROMACS :class: only-dark :width: 500 .. container:: figurelegend Figure: The PEG molecule is a polymer chain made of carbon atoms (in gray), oxygen atoms (in red), and hydrogen atoms (in white). See the corresponding |video_peg_youtube|. .. |video_peg_youtube| raw:: html video .. container:: justify Save *peg.gro* in a new folder. Next to *peg.gro* create an empty file named *topol.top*, and copy the following lines in it: .. code-block:: bw [ defaults ] ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 1 no 1.0 1.0 ; Include forcefield parameters #include "ff/charmm35r.itp" #include "ff/peg.itp" #include "ff/tip3p.itp [ system ] ; Name PEG [ molecules ] ; Compound #mols PEG 1 .. container:: justify Next to *conf.gro* and *topol.top*, create a folder named *ff/*, and copy the following 3 *.itp* files into it: |download_charmm35r.itp|, |download_peg.itp|, and |download_tip3p.itp|. .. |download_charmm35r.itp| raw:: html charmm35r.itp .. |download_peg.itp| raw:: html peg.itp .. |download_tip3p.itp| raw:: html tip3p.itp .. container:: justify These 3 files contain the parameters for the PEG and the water molecules with oxygen (OW) and hydrogen (HW) atoms. .. container:: justify Create an *inputs/* folder next to *ff/*, and create a new empty file called *em.mdp*. Copy the following lines into it: .. code-block:: bw integrator = steep emtol = 10 emstep = 0.0001 nsteps = 5000 nstenergy = 1000 nstxout = 100 cutoff-scheme = Verlet coulombtype = PME rcoulomb = 1 rvdw = 1 pbc = xyz .. container:: justify Most of these commands have been seen in previous tutorials. The most important command is *integrator = steep*, which set the algorithm used by GROMACS as the steepest-descent, which moves the atoms following the direction of the largest forces until one of the stopping criteria is reached :cite:`debye1909naherungsformeln`. .. container:: justify Run the energy minimization using GROMACS by typing in a terminal: .. code-block:: bash gmx grompp -f inputs/em.mdp -c peg.gro -p topol.top -o em-peg gmx mdrun -deffnm em-peg -v -nt 8 .. container:: justify The *-nt 8* option limits the number of threads that GROMACS uses. Adjust the number to your computer. .. container:: justify After the simulation is over, open the trajectory file with VMD by typing in a terminal: .. code-block:: bash vmd peg.gro em-peg.trr .. container:: justify From VMD, the PEG molecule can be seen moving a little by the steepest-descent algorithm. .. container:: justify Before adding the water, let us reshape the box and recenter the PEG molecule in the box. As a first step, let us use a cubic box of lateral size :math:`2.6~\text{nm}`. .. code-block:: bash gmx trjconv -f em-peg.gro -s em-peg.tpr -o peg-recentered.gro -center -pbc mol -box 2.6 2.6 2.6 .. container:: justify Select *system* for both centering and output. The newly created *gro* file named *peg-recentered.gro* will be used as a starting point for the next step of the tutorial. Solvate the PEG molecule ======================== .. container:: justify Let us add the water molecules to the system by using *gmx solvate*: .. code-block:: bash gmx solvate -cp peg-recentered.gro -cs spc216.gro -o peg-solvated.gro -p topol.top .. container:: justify Here *spc216.gro* is a default GROMACS file containing a pre-equilibrated water reservoir. .. container:: justify The newly created file *peg-solvated.gro* contains the water molecules, and a a new line in was added to the topology file *topol.top*: .. code-block:: bw [ molecules ] ; Compound #mols PEG 1 SOL 546 .. container:: justify We can apply the same energy minimization to the newly created solvated system. Simply add the following line to *em.mdp*: .. code-block:: bw define = -DFLEXIBLE .. container:: justify And then launch the energy minimization again using: .. code-block:: bash gmx grompp -f inputs/em.mdp -c peg-solvated.gro -p topol.top -o em gmx mdrun -deffnm em -v -nt 8 .. container:: justify The *define = -DFLEXIBLE* option triggers the following *if* condition within the *tip3p.itp* file: .. code-block:: bw #ifdef FLEXIBLE [ bonds ] ; i j funct length force.c. 1 2 1 0.09572 502416.0 0.09572 502416.0 1 3 1 0.09572 502416.0 0.09572 502416.0 [ angles ] ; i j k funct angle force.c. 2 1 3 1 104.52 628.02 104.52 628.02 .. container:: justify With this *if* condition the water molecules behave as flexible. This is better because rigid molecules and energy minimization usually don't go along well. For the next molecular dynamics steps, rigid water molecules will be used by not including the *define = -DFLEXIBLE* command in the inputs. Equilibrate the PEG-water system ================================ .. container:: justify Let use equilibrate the system in two steps: first a NVT simulation, with constant number of particles, constant volume, and imposed temperature, and second a NPT simulation with imposed pressure. .. container:: justify Within the *inputs/* folder, create a new input named *nvt-peg-h2o.mdp*, and copy the following lines into it: .. code-block:: bw integrator = md dt = 0.002 nsteps = 10000 nstenergy = 500 nstlog = 500 nstxout-compressed = 500 constraint-algorithm = lincs constraints = hbonds continuation = no coulombtype = pme rcoulomb = 1.0 rlist = 1.0 vdwtype = Cut-off rvdw = 1.0 tcoupl = v-rescale tau_t = 0.1 0.1 ref_t = 300 300 tc_grps = PEG Water gen-vel = yes gen-temp = 300 gen-seed = 65823 comm-mode = linear comm-grps = PEG .. container:: justify Most of these commands have already been seen. In addition to the conventional *md* leap-frog algorithm integrator, long-range Coulomb and short-range van der Waals interactions, the LINCS constraint algorithm is used to maintain the hydrogen bonds as rigid. An initial temperature of :math:`300~K` is given to the system by the *gen-* commands, and the PEG is maintained in the center of the box by the *comm-mode* and *comm-grps* commands. .. container:: justify Launch the NVT simulation using: .. code-block:: bash gmx grompp -f inputs/nvt-peg-h2o.mdp -c em.gro -p topol.top -o nvt -maxwarn 1 gmx mdrun -deffnm nvt -v -nt 8 .. container:: justify The *maxwarn 1* option is used to avoid a GROMACS WARNING related to the centering of the PEG in the box. .. container:: justify Let us follow-up with the NPT equilibration. Duplicate the *nvt-peg-h2o.mdp* file into a new input file named *npt-peg-h2o.mdp*. Within *npt-peg-h2o.mdp*, Within the *npt-peg-h2o.mdp*, delete the lines related to the creation of velocity as its better to keep the velocities generated during the NVT run: .. code-block:: bw gen_vel = yes gen-temp = 300 gen-seed = 65823 .. container:: justify In addition to the removal the previous 3 lines, add the following lines to *npt-peg-h2o.mdp* to specify the isotropic barostat with imposed pressure of :math:`1~\text{bar}`: .. code-block:: bw pcoupl = c-rescale pcoupltype = isotropic tau-p = 0.5 ref-p = 1.0 compressibility = 4.5e-5 .. container:: justify Run the NPT simulation, using the final state of the NVT simulation *nvt.gro* as starting configuration: .. code-block:: bash ${gmx} grompp -f inputs/npt-peg-h2o.mdp -c nvt.gro -p topol.top -o npt -maxwarn 1 ${gmx} mdrun -deffnm npt -v -nt 8 .. container:: justify Let us observe the evolution of the potential energy of the system during the 3 successive equilibration steps, i.e. the *em*, *nvt*, and *npt* steps, using the *gmx energy* command as follow: .. code-block:: bash gmx energy -f em.edr -o energy-em.xvg gmx energy -f nvt.edr -o energy-nvt.xvg gmx energy -f npt.edr -o energy-npt.xvg .. container:: justify For each of the 3 *gmx energy* commands, select *potential*. .. figure:: ../figures/level2/stretching-a-polymer/potential-energy-light.png :alt: Potential energy from molecular dynamics simulation in GROMACS :class: only-light .. figure:: ../figures/level2/stretching-a-polymer/potential-energy-dark.png :alt: Potential energy from molecular dynamics simulation in GROMACS :class: only-dark .. container:: figurelegend Figure: Evolution of the potential energy during the 3 equilibration steps, respectively the energy minimization (a), the NVT step (b), and the NPT step (c). .. container:: justify Let us launch a longer simulation, and extract the angle distribution between the different atoms of the PEG molecules. This angle distribution will be used later as a benchmark to probe the effect of of the stretching on the PEG structure. .. container:: justify Create a new input named *production-peg-h2o.mdp*, and copy the following lines into it: .. code-block:: bw integrator = md dt = 0.002 nsteps = 50000 nstenergy = 100 nstlog = 100 nstxout-compressed = 100 constraint-algorithm = lincs constraints = hbonds continuation = no coulombtype = pme rcoulomb = 1.0 rlist = 1.0 vdwtype = Cut-off rvdw = 1.0 tcoupl = v-rescale tau_t = 0.1 0.1 ref_t = 300 300 tc_grps = PEG Water comm-mode = linear comm-grps = PEG .. container:: justify This script resembles the *nvt-peg-h2o.mdp* input, but the duration and output frequency is different, and without the *gen-vel* commands. .. container:: justify Run it using: .. code-block:: bash gmx grompp -f inputs/production-peg-h2o.mdp -c npt.gro -p topol.top -o production -maxwarn 1 gmx mdrun -deffnm production -v -nt 8 .. container:: justify First, create an index file called *angle.ndx* using the *gmx mk_angndx* command: .. code-block:: bash gmx mk_angndx -s production.tpr -hyd no .. container:: justify The *angle.ndx* file generated contains groups with all the atoms involved by an angle constraint, with the exception of the hydrogen atoms due to the use of *-hyd no*. The atom ids selected in the groups can be seen from the *index.ndx* file: .. code-block:: bw [ Theta=109.7_795.49 ] 2 5 7 10 12 14 17 19 21 24 26 28 31 33 35 38 40 42 45 47 49 52 54 56 59 61 63 66 68 70 73 75 77 80 82 84 .. container:: justify Here, each number corresponds to the atom index, as can be seen from the initial *peg.gro* file. For instance, the atom of *id 2* is a carbon atom, and the atom with *id 5* is an oxygen: .. code-block:: bw PEG in water 86 1PEG H 1 2.032 1.593 1.545 0.6568 2.5734 1.2192 1PEG C 2 1.929 1.614 1.508 0.1558 -0.2184 0.8547 1PEG H1 3 1.902 1.721 1.523 -3.6848 -0.3932 -3.0658 1PEG H2 4 1.921 1.588 1.400 -1.5891 1.4960 0.5057 1PEG O 5 1.831 1.544 1.576 0.0564 -0.5300 -0.6094 1PEG H3 6 1.676 1.665 1.494 -2.6585 -0.5997 0.3128 1PEG C1 7 1.699 1.559 1.519 0.6996 0.0066 0.2900 1PEG H4 8 1.699 1.500 1.425 4.2893 1.6837 -0.9462 (...) .. container:: justify Then, extract the angle distribution from the *production.xtc* file using *gmx angle*: .. code-block:: bash gmx angle -n angle.ndx -f production.xtc -od angle-distribution.xvg -binwidth 0.25 .. container:: justify Select 1 for the O-C-C-O dihedral. .. figure:: ../figures/level2/stretching-a-polymer/dihedral-distribution-light.png :alt: Angular distribution from molecular dynamics simulation in GROMACS :class: only-light .. figure:: ../figures/level2/stretching-a-polymer/dihedral-distribution-dark.png :alt: Angular distribution from molecular dynamics simulation in GROMACS :class: only-dark .. container:: figurelegend Figure: Angular distribution for the O-C-C-O dihedral of the PEG molecules. Stretch on the polymer ====================== .. container:: justify Create a new folder named *elongated-box/* next to *cubic-box/*, and copy *ff/*, *inputs/*, *em-peg.gro*, and em-peg.tpr from *cubic-box/* into *elongated-box/*: .. container:: justify To leave space for the stretched PEG molecule, let us create an elongated box of length :math:`6~\text{nm}` along the *x* direction: .. code-block:: bash gmx trjconv -f em-peg.gro -s em-peg.tpr -o peg-elongated.gro -center -pbc mol -box 6 2.6 2.6 .. container:: justify Select *system* for both centering and output. .. container:: justify Then, follow the exact same steps as previously to solvate and equilibrate the system: .. code-block:: bash gmx solvate -cp peg-elongated.gro -cs spc216.gro -o peg-solvated.gro -p topol.top gmx grompp -f inputs/em.mdp -c peg-solvated.gro -p topol.top -o em -maxwarn 1 gmx mdrun -deffnm em -v -nt 8 gmx grompp -f inputs/nvt-peg-h2o.mdp -c em.gro -p topol.top -o nvt -maxwarn 1 gmx mdrun -deffnm nvt -v -nt 8 gmx grompp -f inputs/npt-peg-h2o.mdp -c nvt.gro -p topol.top -o npt -maxwarn 1 gmx mdrun -deffnm npt -v -nt 8 The index file -------------- .. container:: justify To apply a forcing to the ends of the PEG, one needs to create atom groups. Specificaly, we want to create two groups, each containing a single oxygen atom from the edges of the PEG molecules (with ID 82 and 5). In GROMACS, this can be done using and index file *.ndx*. Create a new index file named *index.ndx* using the *gmx make_ndx* command: .. code-block:: bash gmx make_ndx -f nvt.gro -o index.ndx .. container:: justify When prompted, type the following 4 lines to create 2 additional groups: .. code-block:: bash a 82 a 5 name 6 End1 name 7 End2 .. container:: justify Then, type *q* for quitting. The index file *index.ndx* contains 2 additional groups named *End1* and *End2*: .. code-block:: bw (...) [ PEG ] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 [ End1 ] 82 [ End2 ] 5 The input file -------------- .. container:: justify Let us create an input file for the stretching of the PEG molecule. .. container:: justify Create a new input file named *stretching-peg-h2o.mdp* within *inputs/*, and copy the following lines in it: .. code-block:: bw integrator = md dt = 0.002 nsteps = 50000 nstenergy = 100 nstlog = 100 nstxout-compressed = 100 constraint-algorithm = lincs constraints = hbonds continuation = no coulombtype = pme rcoulomb = 1.0 rlist = 1.0 vdwtype = Cut-off rvdw = 1.0 tcoupl = v-rescale tau_t = 0.1 0.1 ref_t = 300 300 tc_grps = PEG Water .. container:: justify So far, the script is similar to the previously created *production-peg-h2o.mdp* file, but without the *comm-mode* commands. To apply the constant forcing to the *End1* and *End2* groups, add the following lines to *production-peg-h2o.mdp*: .. code-block:: bw pull = yes pull-coord1-type = constant-force pull-ncoords = 1 pull-ngroups = 2 pull-group1-name = End1 pull-group2-name = End2 pull-coord1-groups = 1 2 pull-coord1-geometry = direction-periodic pull-coord1-dim = Y N N pull-coord1-vec = 1 0 0 pull-coord1-k = 200 pull-coord1-start = yes pull-print-com = yes .. container:: justify The force constant is requested along the *x* direction only (Y N N), with a force constant :math:`k = 200~\text{kJ}~\text{mol}^{-1}~\text{nm}^{-1}`. .. container:: justify Launch the simulation using the *-n index.ndx* option for the *gmx grompp* command to refer to the previously created index file, so that GROMACS finds the *End1* and *End2* groups. .. code-block:: bash gmx grompp -f inputs/stretching-peg-h2o.mdp -c npt.gro -p topol.top -o stretching -n index.ndx gmx mdrun -deffnm stretching -v -nt 8 .. container:: justify Two data files named *stretching_pullf.xvg* and *stretching_pullx.xvg* are created during the simulation, and contain respectively the force and distance between the 2 groups *End1* and *End2* as a function of time. .. figure:: ../figures/level2/stretching-a-polymer/pull-position-light.png :alt: Pull position from molecular dynamics simulation in GROMACS :class: only-light .. figure:: ../figures/level2/stretching-a-polymer/pull-position-dark.png :alt: Pull position from molecular dynamics simulation in GROMACS :class: only-dark .. container:: figurelegend Figure: Distance between the two pulled groups *End1* and *End2* along the *x* direction, :math:`\Delta x`, as a function of time :math:`t`. .. container:: justify It can be seen from the evolution of the distance between the groups, :math:`\Delta x`, that the system reaches its equilibrium state after approximately 20 pico-seconds. .. container:: justify Let us probe the effect of the stretching on the structure of the PEG by remeasuring the dihedral angle values: .. code-block:: bash gmx mk_angndx -s stretching.tpr -hyd no -type dihedral gmx angle -n angle.ndx -f stretching-centered.xtc -od dihedral-distribution.xvg -binwidth 0.25 -type dihedral -b 20 .. container:: justify Select 1 for the O-C-C-O dihedral. Here the option *-b 20* is used to disregard the first 20 pico-seconds of the simulation during which the PEG has not reach is final length. .. figure:: ../figures/level2/stretching-a-polymer/comparison-dihedral-distribution-light.png :alt: Angular distribution from molecular dynamics simulation in GROMACS :class: only-light .. figure:: ../figures/level2/stretching-a-polymer/comparison-dihedral-distribution-dark.png :alt: Angular distribution from molecular dynamics simulation in GROMACS :class: only-dark .. container:: figurelegend Figure: Angular distribution for the O-C-C-O dihedral of the PEG molecules, comparing the unstretched (cyan) and stretched case (orange). .. container:: justify The change in dihedral angles disribution reveals a configurational change of the polymer induced by the forcing. This transition is called gauche-trans, where *gauche* and *trans* refer to possible states for the PEG monomer :cite:`binder1995monte, liese2017hydration`. .. figure:: ../figures/level2/stretching-a-polymer/gauche-trans.png :alt: Illustration of the gauche (left) and trans (right) states of the PEG polymer. .. container:: figurelegend Figure: Illustration of the gauche (left) and trans (right) states of the PEG polymer. .. include:: ../../non-tutorials/accessfile.rst