GROMACS version: gmx, version 2020.1-Ubuntu-2020.1-1
GROMACS modification: Yes/No
Here post your question I have trying to generate my itp and top file using the pdb2gmx command. But as my structure does not have any parameters in the .rtp file of the forcefield I am getting an error.
So I want to know how can I add a residue into the gromacs forcefield, so that I can generate the .itp and .top.
I know the values of all the parameters required, like the force constants of each of each of the bonds, or the angles or dihedral parameters etc.
I looked into gromacs documentation, but it is not quite clear. If someone could show an example on how to do this, then it would be very helpful.
1 Like
Hi,
The information may be a little scattered in GROMACS manual and example may be useful.
I guess you have looked at Adding a Residue to a Force Field — GROMACS 2021.3 documentation
Also maybe the tool x2top or other online/offline tools may be useful to build the step (connectivity for the new residue) that you have to plug in the rtp file. You can also look at the third party category in the forum ( Third party tools and files - GROMACS forums) if some useful tool is there.
\Alessandra
I tried x2top with a simple methane molecule with the values of all the parameters as given in the instructions present in the link gmx x2top — GROMACS 5.1.1 documentation
Here is my pdb file for methane:
COMPND UNNAMED
AUTHOR GENERATED BY OPEN BABEL 2.3.90
HETATM 1 C UNL 1 -2.432 0.404 -0.000 1.00 0.00 C
HETATM 2 H UNL 1 -1.323 0.404 0.000 1.00 0.00 H
HETATM 3 H UNL 1 -2.802 0.634 -1.020 1.00 0.00 H
HETATM 4 H UNL 1 -2.802 1.173 0.709 1.00 0.00 H
HETATM 5 H UNL 1 -2.802 -0.594 0.312 1.00 0.00 H
CONECT 1 2 3 4 5
CONECT 2 1
CONECT 3 1
CONECT 4 1
CONECT 5 1
MASTER 0 0 0 0 0 0 0 0 5 0 5 0
END
and this is what I got after I gave the x2top command
gmx x2top -f methane.pdb -o out.top -r out.rtp -v -nexcl 3 -alldih no -remdih yes -name UNL -param yes -pdbq no -round yes -kb 10000 -kt 400 -kp 5
WARNING: all CONECT records are ignored
Opening force field file /usr/share/gromacs/top/oplsaa.ff/atomname2type.n2t
There are 23 name to type translations in file oplsaa.ff
Generating bonds from distances…
atom 5
Can not find forcefield for atom C-1 with 4 bonds
Program: gmx x2top, version 2020.1-Ubuntu-2020.1-1
Source file: src/gromacs/gmxpreprocess/x2top.cpp (line 188)
Fatal error:
Could only find a forcefield type for 4 out of 5 atoms
For more information and tips for troubleshooting, please check the GROMACS
website at Errors - Gromacs
Hi,
Is the molecule you want to describe is a methane or similar? Maybe it is easy to build directly an itp file and include the parameters there. Maybe someone else has already a template for methane.
\Alessandra
No…mine is a complicated molecule. I was just trying with a simple example… Just wanted to about the procedure. One thing I want to ask is… in the rtp files for gromos, there are terms like gb_42 or ga_14 in the bonds and angles section. What does that mean? Here I am sharing that portion…
[ ALA ]
[ atoms ]
N N -0.31000 0
H H 0.31000 0
CA CH1 0.00000 1
CB CH3 0.00000 1
C C 0.450 2
O O -0.450 2
[ bonds ]
N H gb_2
N CA gb_21
CA CB gb_27
CA C gb_27
C O gb_5
C +N gb_10
[ angles ]
; ai aj ak gromos type
-C N H ga_32
-C N CA ga_31
H N CA ga_18
N CA CB ga_13
N CA C ga_13
CB CA C ga_13
CA C O ga_30
CA C +N ga_19
O C +N ga_33
[ impropers ]
; ai aj ak al gromos type
N -C CA H gi_1
CA N C CB gi_2
C CA +N O gi_1
[ dihedrals ]
; ai aj ak al gromos type
-CA -C N CA gd_14
-C N CA C gd_39
N CA C +N gd_40
as shown above…gb_27, gd_14… if someone could answer me even that then also I think I can try generating the topology
Hi,
the labels gb_xx ga_xx gd_xx corresponds to the bond, angle and dihedral parameters defined in the ffbonded.itp. You find the file ffbonded.itp in the share/gromacs/top/{forcefield_name}.ff/ffbonded.itp where {forcefield_name} is the force field you are looking at.
ffbonded.itp and ffnonbonded.itp are included in forcefield.itp
I hope it helps
Alessandra
I tried adding the parameters in the .rtp file as well as the non bonded parameters in the corresponding ffnonbonbed.itp file but still I am facing an error. I am working with a silica slab, so my system has large number of atoms. Now, as far as I know the parameters for a unit cell should work for the entire periodic lattice, but when I am running the pdb2gmx command, only four atoms are presented in the final gro and top file, all other atoms are removed. Moreover, the bonding between these files are also getting removed. When viewed in VMD I am only getting four isolated atoms, 1 Si, 2 O atoms and 1 H atom
To create a topology for a surface or other crystalline material that has non-linear connectivity, you should not be using pdb2gmx, you should be using x2top.
Ok… Is the procedure same for both? I am stating what I want to achieve here. My purpose is to generate a silica slab (which has a periodicity in it’s structure). Hence, I want to keep the terminal atoms(present at the edge of the simulation box) to have their connectivities open, as they would be generating a continuous sequence of atoms when pbc conditions are applied.
Now, as my system has a large number of atoms I simply cannot use the external servers because they have a limitation on the number of atoms that their servers can process. Hence, I am struck with generating the topology manually.
I have added the atoms, bonds, angles and dihedrals in the .rtp file and the LJ parameters in the nonbonded.itp file. What should I do next?
Explicit inclusion of the parameters in the .rtp file is not required. All bonded and nonbonded terms can be read from ffbonded.itp and ffnonbonded.itp.
Ok… So I can just add the parameters in the ffbonded.itp and ffnonbonded.itp?
I tried adding the parameters in the corresponding files and ran the pdb2gmx command but only 4 atoms were present in the final gro file and all other atoms (from the original silica slab) were removed. What was more annoying was that those 4 atoms were not connected by any bonds. All I got was 4 atoms isolated in space.
Here I am attaching the files I added in the forcefield files:
1)Added in the .rtp file:
[ SIL ]
[ atoms ]
; name type charge chargegroup
SI SI 1.020 0
OW OW -0.510 1
H H 0.255 0
O O -0.993 1
[ bonds ]
; atom1 atom2 c0 c1
SI OW 0.1720 2.7042e+06
O H 0.0938 2.9551e+07
SI O 0.1690 1.9257e+06
[ angles ]
; atom1 atom2 atom3 th0 cth
SI O SI 166.00 64885.36
SI OW H 119.00 2211.40
O SI O 180.00 500.00
OW SI O 109.50 450.00
[ dihedrals ]
; atom1 atom2 atom3 atom4 phi0 cp mult
H OW SI O 180.00 1.00 3
O SI O SI 0.00 1.05 3
-
Added in the ffnonbonded.itp file:
SI 14 0.000 0.000 A 5.803853e-03 5.35358e-06
O 8 0.000 0.000 A 2.24192e-03 1.36137e-06
OW 8 0.000 0.000 A 2.24192e-03 1.36137e-06 -
the final gro file that I got:
3D Atomistic
4
1SIL SI 1 2.169 7.014 -0.182
1SIL OW 2 6.337 -0.752 -0.456
1SIL H 3 0.494 5.754 -0.455
1SIL O 4 2.952 -1.104 -0.091
5.84310 8.11840 0.36460
I am attaching all these for your reference. The pdb file I am attaching with this message (CONVERTED TO .dat format)topol.top (2.4 KB)
silica_slab_edited_OW.dat (718.4 KB)
Correct.
Like I said before, this is not the correct approach; it will never work. Use x2top with a suitable .n2t file.
Ok. So basically the pdb2gmx command is applicable for non-periodic systems. If I were to use a simple ligand, then this would have worked, right?
And one more thing, I tried looking for the format of n2t file…This is what came up:
H H 0.408 1.008 1 O 0.1
O OA -0.674 15.9994 2 C 0.14 H 0.1
C CH3 0.000 15.035 1 C 0.15
C CH0 0.266 12.011 4 C 0.15 C 0.15 C 0.15 O 0.14
I didn’t quite understand what column 6 and onwards mean? And my research paper mentions about specific parameters for the LJ potential. How can I add these into my forcefield? I mean when I run x2top, does this take only the .n2t files or does it also takes the parameters given in the .rtp and ffbonded.itp and ffnonbonded.itp files as well?
I tried the x2top command. Here is what I wrote in the atomname2types.n2t file:
H H 0.2550 1.008 1 OW 0.0938
SI SI 1.0200 28.08 4 O 0.1690 OW 0.1720 O 0.1690 O 0.1690
OW OW -0.510 15.9994 2 SI 0.1720 H 0.0938
O O -0.993 15.9994 2 SI 0.1720 SI 0.1720
This is the error that I got… I have modified my pdb file this time. Please have a look (again I have uploaded it in slab.dat )slab.dat (684.2 KB)
Program: gmx x2top, version 2021.3
Source file: src/gromacs/gmxpreprocess/x2top.cpp (line 186)
Fatal error:
Could only find a forcefield type for 2560 out of 4608 atoms
For more information and tips for troubleshooting, please check the GROMACS
website at Errors - Gromacs
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
See https://manual.gromacs.org/current/reference-manual/file-formats.html#n2t
Pair-specific interactions are listed in [nonbond_params] in ffnonbonded.itp.
This means your .n2t file does not cover some of the bonds that need to be made.
One question I have to ask here. My system is a periodic one right. So do I have to cap the terminal Si atoms with any H or any ghost atom or sort? Because the terminal atoms still have their valencies open (as they are periodic cell and they will get connected when the pbc condition is applied).
x2top applies periodicity, so if that’s the way you want to model them in the simulation, that’s how the topology has to be created, not with any capping groups.
I tried the x2top command but I got an error. Actually I had my terminal silicon atoms uncapped, which resulted in an error, but when I capped the terminal atoms with H, the command ran successfully and I got the corresponding the .top and .rtp file. So now I am quite unsure about how the command works
There’s not a lot I can tell based on the information provided except that x2top reads through each atom in the coordinate file and attempts to assign atom type, charge, and bonded connectivity based on the information you have provided in the .n2t file. Distance tolerances are ±10% as is standard for all GROMACS distance searches. If it does not find a match, you get an error. Likely your uncapped structure had interatomic distances that were beyond the threshold.
So let me get this straight, all I have to do is add my atom name and the corresponding atom types, followed by their charge, masses and connectivities (in terms of bond lengths) and I am good to go. Is that right?
And if I have certain atoms which are not connected to other atoms the way I mentioned it, then it will throw an error. Is that right?
Then can it be possible that the terminal atoms which have their connectivities open (as they are connected to the atoms of the next box) might be leading to the problem that I am facing?