-
Notifications
You must be signed in to change notification settings - Fork 24
Tutorial: Martini Polymers
Polyply facilitates the easy generation of polymer melts / amorphous blends, polymer solutions, and even liquid-liquid phase separated systems. In this tutorial, we are going to use polyply to generate itp files and starting configurations for these three cases. In particular, you will generate an amorphous polymer simulation using pure PEO and a blend of the two other polymers. Subsequently, we turn polymer solvent systems and last but not least liquid-liquid phase separated systems of Dextran and PEO. Note that all basic steps apply in the same manner to all-atom models, but make sure to also check out the GROMOS tutorial for all-atom simulations.
First, we need to generate an itp file for the polymers that we want to study. Polyply is shipped with a range of polymer input files collected in libraries, which can be used to generate itp files for polymers. Note that you can also implement your own polymer definition. To do so check out the tutorial on generating your own polymer library. If you have any trouble always feel free to raise an issue. To inspect which libraries are available within polyply run:
polyply -list-lib
You should see that a library for martini3 polymers is implemented. Next, you can check that the building block you require is present in the library by running:
polyply -list-blocks martini3
To check any other library simply replace martini3 with the library name. Running the command tells us that there are parameters for both PS and PEO present. Note that there are not always links to combine different polymers with each other. Next, we want to invoke the gen_params sub-program to generate an itp file. For a single polymer, the command has the following basic structure. Run polyply gen_params -h for more options.
polyply gen_params -lib <library_name> -o <filename>.itp -name <name of polymer> -seq [<monomer-name>:<number of monomers>, ...]
To generate PEO of length 200 we simply need to provide the library name using -lib, a name for the itp file using -o, a name using -name and the sequence of the polymer using -seq. In this case, the sequence is just 200 residues of PEO (i.e. -seq PEO:200). The full command would be:
polyply gen_params -lib martini3 -o PEO200.itp -name PEO -seq PEO:200
This generates the required itp file. To do the same for PS simply change PEO to PS. If you need more complex linear polymer-sequences you can provide a text file with a list of residues and provide it using the -seqf option. Branched polymers are also possible to generate itp files for but you have to provide a networkx .json describing the graph of the branched polymer. This discussion offers an example. The file also has to be provided using the -seqf option.
To generate starting coordinates all you need is the target topology file. First, we need the interaction parameters. Those you can download here. Subsequently, write the topology file. In our case let's generate an amorphous system of 200 polymers. In this case, the topology file looks as follows:
#include "martini_v3.0.0.itp"
#include "PEO200.itp"
[ system ]
; name
Melt of PEO 200
[ molecules ]
; name number
PEO 200
The martini interaction itp file (i.e. martini_v3.0.0.itp) can be downloaded from here. To create starting coordinates we call the gen_coords sub-program with the general command looking as follows:
polyply gen_coords -p <topfilename> -o <filename>.gro -name <name> -dens <density>
In order to generate coordinates you have two options for setting the system size: (1) you set the box size (i.e. -box <x> <y> <z>); keep in mind only rectangular boxes are currently supported; (2) you provide the overall system density using -dens. Most of the time it is convenient to provide the target system density. Note that sometimes it is convenient to set the target density a little lower than what is expected in the simulation to gain a speed up in the structure generation. To further speed up the generation you can install the numba package. For our example, we will generate PEO at a target density of 1000 kg/cm3, which is the Martini system density. The overall command for our polymer amorphous system becomes:
polyply gen_coords -p system.top -o melt.gro -name PEO -dens 1000
This should take less than 60 seconds on a single CPU provided you installed numba. Once you have generated your melt you will need to first run an energy minimization, followed by a short NpT equilibration. You can download the files for minimization and equilibration here. The GROMACS commands are as follows:
gmx grompp -f min.mdp -c melt.gro -p system.top -o min.tpr
gmx mdrun -v -s min.tpr -deffnm min
gmx grompp -f pre_NpT.mdp -c min.gro -p system.top -o eq.tpr
gmx mdrun -v -s eq.tpr -deffnm eq
Note to run the system as a polymer melt you need to increase the temperature above the PEO melting temperature (~335K).
To generate the polymer blend you simply have to generate an itp file for PS as well. Subsequently, simply add it to your topology file and run the command above again. This will generate the polymer blend. Edit the topology file to include PS:
#include "martini_v3.0.0.itp"
#include "PEO200.itp"
#include "PS200.itp"
[ system ]
; name
A blend of PEO and PS
[ molecules ]
; name number
PEO 200
PS 200
Now run the coordinate generation again:
polyply gen_coords -p system.top -o blend.gro -name PEO_PS_blend -dens 1000
Again energy minimization and equilibration are needed to get the system ready for production.
If you already generated an itp file using the gen_params tool described in Section 2, a polymer solution is generated as simply as the polymer blend. First, however, a solvent has to be selected. Apart from water, Martini3 offers a whole host of solvents. A comprehensive list can be found in this file, which you will also need to download and include in the topology. Additional mostly aromatic solvents can be found here.
For example, PS is soluble in acetone but insoluble in water. To generate and compare both systems, we can simply write the following two topology files and generate the coordinates as before. Note that the concentration is adjusted by specifying an appropriate number of solvent molecules. For this example, we use 8000 acetone molecules, a 1:4 ratio of monomers to solvent. In subsequent commands we assume the file name to be sol_sys.top
#include "martini_v3.0.0.itp"
#include "martini_v3.0.0_solvents_v1.itp"
#include "PS200.itp"
[ system ]
; name
A blend of PEO and PS
[ molecules ]
; name number
PS 10
PPN 8000
The command for generating the system is pretty much the same as before. Note that the Martini Actone density is about 970
polyply gen_coords -p sol_sys.top -o sol_sys.gro -name PS_PPN -dens 970
To generate the water system, just replace PPN with W and change the density to 1000
Dextran and PEO form an LLPS in when solvated in water above a critical concentration. Martini can capture these phase effects. This particular system has already been published. To set up the system we first have to generate the Dextran itp file. Dextran with a number average molecular weight of 65 is hardly branched and forms an LLPS with PEO. Thus we generate linear Dextran:
polyply gen_params -seq DEX:65 -o dextran.itp -name DEX -lib martini3
Subsequently, we need to build the system above the critical concentration. To speed up phase separation and keep the system size small, we will choose a rather high concentration in a small system. Note that for production the system sizes should not be too small. Write the topology file called llps.top:
#include "martini_v3.0.0.itp"
#include "martini_v3.0_solvents.itp"
#include "dextran.itp"
#include "PEO.itp"
[ system ]
LLPS of dextran PEO in water
[ molecules ]
DEX 18
PEO 20
W 13693
Next, we generate the coordinates for the mixed system as done in the other examples. Only in this case, we lower the maximum allowed force in the random-walk process by setting -mf 5000. It is not strictly needed but helps with avoiding the occasional overlap.
polyply gen_coords -p llps.top -o llps.gro -name llps -dens 1000 -mf 5000
After energy minimization, you should see a largely mixed system of PEO and Dextran in water. To facilitate phase separation you can change the barostat settings to semi-isotropic. The phase separation takes about 500ns before it becomes visible.
- run energy minimizations with the -DFLEXIBLE option to use flexible bonds
- energy minimizing using GROMACS double precision can help to relax close contact configurations better
- run equilibration with Berendsen barostat to relax to the target volume or Crescale barostat for isotropic systems
- the solvent density in Martini is not always the same as the real density, because Martini beads have a fixed mass. Compare the solvent mass to the real molecule mass to see if the density has to be scaled