This task is devoted to introduce GROMACS notation for force fields. We build the model of hydrazine molecule, introducing lone electron pair as an artificial atom. We set molecular conformations with pairs or by defining dihedral potential. Varying dihedral potentials or pairs we obtain different equilibrium conformations of hydrazine.
1. Draw hydrazine molecule in PyMOL. Add +1 charge to nitrogen atoms and additional protons will appears. The resulting PDB file should be converted into GRO/G96 fromat:
gmx editconf -f hydrazine.pdb -o hydrazine.gro -box 3 3 3
-box sets cell dimenensions. We will not use periodic conditions but the dimensions should be set up.
2. Rename atoms and reshuffle atom strings in GRO file:
1HZN N1 1 1.555 1.546 1.490
1HZN H11 2 1.571 1.599 1.574
1HZN H12 3 1.533 1.609 1.414
...
1HZN - residue number (1) and name (HZN). N1,H11,.. - atom names, they should be different and well-knowable. 1,2,3 — atom numbers. Make them consequential after reshuffling strings. Keep position of nubers and words as GRO is positional format (like PDB). One hydrogen atom near every nitrogen must be renamed (e.g. X1/X2) and will serve as a lone pair.
3. Create stub topology file hydrazine.itp with the following blocks:
; semicolon means comments [defaults] ; common settings of the force field [atomtypes] ; list of atom types with default VdW [bodtypes] ; list of bond types with bond potentials [angletypes] ; list of angle types with angle potentials [dihedraltypes] ; list of dihedral types with dihedral potentials [moleculetype] ; name nrexcl [atoms] ; list of atoms in molecule [bonds] ; list of bonds in molecule [angles] ;list of angles in molecule [dihedrals] ;list of dihedrals in molecule [pairs] ;list of 1-4 pairs
Now let study the content of every section. Do not rewrite example lines to your file without editing. All these lines are stub and don’t work.
4. Define defaults section.
[defaults] 1 2 no 0.0 0.0
1 — type of VdW potential function, 2 — type of combination rules, yes — automatic generation of 1-4 pairs, 0.8 — fudgeLJ, 0.5 — fudgeQQ.
Automatic pair generation is used when pdb2gmx is executed. It does not affect in our case.
Type of VdW potential may switch constant definition between C6/C12 and σ/ε. It is easier to use the second notation (number 1) because σ corresponds double atomic radius and ε — to the deepness of the energetic well.
It is hard to observe an effect of combination rules in this task.
fudgeLJ/fudgeQQ scales VdW and coulomb interaction defined for 1-4 pairs. It is important for us if we will use pairs.
5. Define atomtypes section.
You must define three atom types of your molecule nitrogen, hydrogen and lone pair.
[atomtypes] NZ 7 14.0 0.0 A 0.1 0.2 ...
NZ — name of atomtype, 7 — atomic number, 14 — atomic mass, 0.0 — default parital charge, A — atom or dummy, 0.1 — first VdW parameter (C6 or σ), 0.2 — second VdW parameter (C12 or ε). Set first to double radius in nm and second to allowable energy of interaction (in kJ/mol).
Note, do not use zero mass for lone pair. Give some small mass to it.
6. Define bondtype section
[bondtypes] NZ HZ 1 0.5 150 ...
NZ, HZ — atomtypes of atoms constituting the bond; 1 — type of bond potential, the unity means parabolic potential; 0.5 — length of the equilibrium bond in parabolic potential (nm); 150 — spring constant of the bond in parabolic potential (kJ/mol/nm2).
Note that spirng constant with hydrogens is much higher. The lower spring constant allow the bond fluctuate with the bigger amplitude. The allowable order of bond spring constants is about 105 kJ/mol/nm2.
7. Define angletype section.
[bondtypes] HZ NZ HZ 1 109 150 ...
Everything is the same as we said previously for bondtype. We give equilibrium angle in degrees and spring constant in kJ/mol/rad2 (it is not a misprint! degrees and radians!). The allowable order of angle spring constants is about 102-103 kJ/mol/rad2.
7. Define angletype section.
[dihedraltype] HZ NZ NZ HZ 3 1.0 -1.0 2.0 -2.0 3.0 -3.0 ...
There are different ways to define the dihedral potential but here we use only one: Ryckaert-Bellemans notation.
All constants are measured in kJ/mol. For our molecule we can write 9 quads of atoms because we have 3 atoms at one nitrogen and 3 atoms at the another one. There are two ways in defining dihedral potential. First way is to use one dihedral potential and to define it for one fourth of atoms. The second way (mostly used) is to define dihedrals for every quad. I’m not going to write the exact numbers you should put into the potential. But I give you plot U(φ) for the following sets of C-s: -1, 1, 0, 0, 0, 0; -1, 0, 1, 0, 0, 0; 1, 3, 0, -4, 0, 0.
The allowable range of dihedral barriers is 1-100 kJ/mol.
8. Define moleculetype section.
[moleculetype] HZN 3
Just write your residue name, space and number 3. There are no variants for large molecules. 3 means that 1-2, 1-3 and 1-4 interactions are normally excluded from non-bonded potential calculations. We may add 1-4 pairs back but it is completely different story and these interaction will be automatically scaled by fudge coefficients. If number 3 is written here 5th, 6th atom, and further neighbourhood interact the atom in the same way that atoms from another molecules do.
9. Define atoms section.
[atoms] 1 NZ 8 HZN N1 2 0.0 ...
All atoms of the molecule should be listed in a correct order. 1 - first number is order number; NZ - atomtype name defined in [atomtypes]; 8 - second number is residue number. When the molecule consists of one residue, just leave the unity in all lines; HZN — name of your residue from [moleculetype]; N1 — name of your atom from PDB/GRO; 2 — charge group. Mark small chemical structures with summary charge close to zero as one charge group (e.g. NH2-, CH3-, -CH2-, >C=O, etc); 0.0 - partial charge of the atom.
Partial charges should be accurately computed but here you can guess them. 1 = electron charge. Mostly, partial charge absolute values are much less than 1.
You can take partial charges from ab initio calculation of the hydrazin molecule (Hydrazin 631*/B3LYP [ORCA] output). You can find charges in «MULLIKEN ATOMIC CHARGE» and «LOWDIN ATOMIC CHARGE» sections.
10. Define bonds section.
[bonds] 1 2 1 ...
1 — index of the first atom in bond, 2 - index from the second atom in bond. The order doesn’t matter. Indexes should be the same as in [atoms] group; the last 1 — the type of bond potential. If the potential for bond between these 2 atomtypes and defining with this potential type exists in [bondedtypes], then its parameters will be used.
All valence bonds should be listed. In our case, do not forget about bond between the nitrogen and the lone pair.
11. Define angles section.
[angles] 1 2 3 1 ...
Meaning of the numbers is the same as in [bonds] sections. It is customary to define all possible angles in the molecule. It is excessive in terms of internal coordinates but it is common.
12. Define dihedrals section.
[dihedrals] 1 2 3 3 ...
Dihedral section is often decorated in the same way as angle section. All atom quads that form dihedral angles are typically listed. However, there is another way when only one atom quad defining the dihedral is listed. The potential should be designed for one or another variant.
Trick: It is possible to use zero for HNN* dihedrals and non-zero for XNNX dihedral (X = lone pair).
13. Define pairs section.
[pairs] 1 2 1 ...
Pairs are decorated in the same way as bonds. It is possible to form energy profile of dihedral rotation using pairs section. When we turn on the pair for 1-4 neighbors, the VdW and coulomb potential between these atoms are turned on. In our case, you may consider to turn on X-X pair — it lead to repulsion (both vdw and coulomb) between lone pairs.
14. Make topology file and include the toy force field.
#include "hzn.itp" [system] Hydrazine [molecules] HZN 1
The arbitrary name can be defined in [system] section. Use residue name in [molecules] section the same as in [moleculetype] section.
15. Let’s perform geometry optimization with this topology.
integrator = steep nsteps = -1 emtol = 1 pbc = no cutoff-scheme = group rvdw = 1.0 rcoulomb = 1.0 rlist = 1.0 nstxout = 1 nstcalcenergy = 1 nstenergy = 1 nstlog = 10
The construction of MDP file is out of our problem’s tasks. Here is geometry optimization file that make GROMACS write every frame of the optimization into trajectory file (TRR) and write the energy values into the binary file (EDR). emtol may be changed to vary the deepness of optimization algorithm.
16. View the result.
gmx grompp -c file.gro -p file.top -f file.mdp -o emin.tpr gmx mdrun -s emin.tpr -deffnm EM0 -v
Open initial gro after running these commands with PyMOL and load the trajectory:
load_traj EM0.trr, file
file — is PyMOL object name. Use «interval=10″ to skip most frames in the trajectory if it is too big.
Compare your structure with ab initio optimized hydrazin.