-
Notifications
You must be signed in to change notification settings - Fork 24
Tutorial: writing .ff input files
under construction
In this tutorial we will learn how to write vermouth (.ff) input files, which are used by polyply to generate polymer itp-files. This tutorial covers standard cases and concepts. For an in-depth view see the documentation on the file format. As a hands on example we will generate input files for atomistic polyethylene glycol in the GROMOS 2016H66 force-field.
Polyply is a residue oriented program. This means it is based upon the fact that polymers can be split into residues (i.e. monomers) and those residues are then to be connected. The input files reflect this notion in the sense that they exists two simple directives. Here we call a directive anything that has the following format [ directive ]. If you are familiar with the GROMACS itp files you will know that to define the a molecule you need to start by writing a [molecultype] directive followed by a [atoms] directive, a [bonds] directive and so forth. The vermouth input file format is very similar to this style of writing input files. Coming back to polymers: We need in general two directive types. The [moleculetype] directive is used to define all atoms and bonded interactions that are within a single residue. On the other hand the [link] directive specifies all bonded interactions that link two residues together. Both of them together will enable you to write custom polymer itp-files in a very efficient manner. First we see how to write the [moleculetype] directive.
Residues (i.e. monomers in a polymer) are defined by writing the [moleculetype] directive. It follows the general format of a GROMACS molecule itp file. So writing a input file for a residue is like writing the itp-file for that residue. The following itp file is the PEG repeat unit file for atomistic PEO. Note that the repeat unit here is taken to be -[CH2-O-CH2]-.
[ moleculetype ]
PEO 3
[ atoms ]
1 CH3 1 PEO EC1 1 0.29 14.0
2 OE 1 PEO O1 2 -0.58 16.0
3 CH3 1 PEO EC2 3 0.29 14.0
[ bonds ]
1 2 2 gb_18
3 2 2 gb_18
[ angles ]
3 2 1 2 ga_15
There are few important things to note: (1) all atom-names need to be unique within one residue; (2) residue name and molecule name should be the same; (3) you only need to provide the atom-indices with the bonded directives. For example, you can write the bond between atom 1 and atom 2 as 1 2 2 gb_18 or alternatively also as 1 2 2 1.4300000e-01 8.1800000e+06, which is the GROMOS definition for that bond or you can also write it as a blank (i.e. 1 2 2). This gives full freedom to write parameters and is common to both link and residue definitions. Another thing to note about this PEG residue is that we have written atom 1 and 3 as CH3, meaning the itp file is equivalent to dimethyl ether. This is important because if we call polyply to ask for a single residue we will automatically get it as CH3 terminated. Later when writing the links we will account for the fact that the atom-type in the polymer should be CH2. For the residue definition we are already done so we move on to the links.
Next we will have to define a number of links which specify, how we string together the PEO residues. The link definition begins by writing the the link directive as follows:
[ link ]
resname PEO
Note that if your link can apply to multiple residues you can write resname "<resname1>|<resname2>|<resname3>" etc., but in this case we only concern ourselves with PEO. In general it is recommended to write every bonded directive into it's own link definition. That is writing a seperate link directive for all the bonds, angles, dihedrals. For PEO we will need bonds, angles, dihedrals and pairs. Let us start with the bond. When two PEOs connect we expect there to be a bond between the atoms EC1 and EC2 of the next residue. So we write the following link entry:
[ link ]
resname "PEO"
[ atoms ]
EC2 { }
+EC1 { }
[ bonds ]
+EC1 EC2 2 gb_27
In this link entry we note the atom names involved in the interaction in the atoms directive followed by empty curly braces. Then we define the bonds directive in the same was as in an itp file, however, the atom names (e.g. EC2) are used instead of the numbers. In addition the relative residue ID needs to be defined. This link is supposed to apply between the atom EC2 of one residue and the to the atom EC1, which belongs to the next residue. The next residue has a resid +1 higher so we use the prefix '+'. However, we still have the atomtypes CH3. To take care of that we add the keyword "replace" in the atoms directive as follows: [ link ] resname "PEO" [ atoms ] EC2 {"replace": {"atype": "CH2"}} +EC1 {"replace": {"atype": "CH2"}} [ bonds ] +EC1 EC2 2 gb_27 The "replace" directive will replace the atom-types of EC2 and +EC1, whenever this link is applied. In this way when we get a bond between to residues we also change the atom-types. Next we define the links for angles in the same way as for bonds. Nothing new about those. Note that we can omit the atoms directive in this case, because we don't use it.
[ link ]
resname "PEO"
[ angles ]
EC2 +EC1 +O1 2 ga_15
[ link ]
resname "PEO"
[ angles ]
+EC1 EC2 O1 2 ga_15
Finally PEG also requires dihedrals. In GROMOS there are three sets of three dihedrals defined for the polyether bonds. Each set of dihedrals describes one torsion and becomes it's own link as shown below.
[ link ]
resname "PEO"
[ dihedrals ]
EC1 O1 EC2 +EC1 1 gd_42 {"version":"1"}
EC1 O1 EC2 +EC1 1 gd_43 {"version":"2"}
EC1 O1 EC2 +EC1 1 gd_44 {"version":"3"}
[ link ]
resname "PEO"
[ dihedrals ]
O1 EC2 +EC1 +O1 1 gd_45 {"version":"1"}
O1 EC2 +EC1 +O1 1 gd_46 {"version":"2"}
O1 EC2 +EC1 +O1 1 gd_47 {"version":"3"}
[ link ]
resname "PEO"
[ dihedrals ]
EC2 +EC1 +O1 +EC2 1 gd_42 {"version":"1"}
EC2 +EC1 +O1 +EC2 1 gd_43 {"version":"2"}
EC2 +EC1 +O1 +EC2 1 gd_44 {"version":"3"}
In general interactions overwrite each other. For example, taking the last set of dihedral we would usually only get gd_44, because the other two dihedral interactions are overwritten. There can only be one interaction of the same type for each set of indices. To circumvent this problem the "version" keyword can be used between curly braces. This way we tell polyply not to overwrite those dihedrals but keep all three of them. This is a general feature and can be used for any link directive. Last but not least we need to introduce pair interactions. Pairs are almost the same as the interactions we have seen before.
[ link ]
resname "PEO"
[ atoms ]
EC1 { }
O1 { }
EC2 { }
+EC1 { }
+O1 { }
+EC2 { }
[ pairs ]
EC1 +EC1 1
O1 +O1 1
EC2 +EC2 1
[ edges ]
EC1 O1
O1 EC2
EC2 +EC1
+EC1 +O1
+O1 +EC2