Hello everyone,
I want to just simulate a box full of dimethoxyethane (dme) molecules in gromacs. At first, after wasting a week on just shifting columns and modifying spaces in pdb files (this is so so so unnecessarily cumbersome and complicated in gromacs), i was finally able to create a box containing 64 dme solvent molecules called dme_box.gro (file attached) from the dme.pdb. Then I rescaled this box to get a smaller one called dme_scaled.gro so as to match the density of the dme solvent using the command:
$ gmx_mpi editconf -f dme_box.gro -density 868.3 -o dme_scaled.gro > logfile
The dme_scaled.box has all bond lengths for the molecule corrupted. I do not understand why.
Anyway, I wanted to generate the topology of the dme molecule with a forcefield file dme.itp which I placed in oplsaa.ff directory that i created, and also wrote #include "dme.itp" in the forcefield.itp file. I edited the residuetypes.dat to include DME as Protein and also modified the amoniacids.rtp file to include the forcefield data from the file dme.itp using the tyle that was there. All files are attached. Now, I want to use this goddamned pathological routine called pdb2gmx. When I use:
$ gmx_mpi pdb2gmx -f dme_scaled.gro -ff oplsaa
The kind of output i get is :
===========================================
Using the Oplsaa force field in directory ./oplsaa.ff
Opening force field file /usr/local_rwth/sw/gromacs/gromacs-5.1.2/intel/openmpi/share/gromacs/top//./oplsaa.ff/watermodels.dat
Select the Water Model:
1: TIP4P TIP 4-point, recommended
2: TIP4PEW TIP 4-point with Ewald
3: TIP3P TIP 3-point
4: TIP5P TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
5: TIP5P TIP 5-point improved for Ewald sums
6: SPC simple point charge
7: SPC/E extended simple point charge
8: None
8
Reading dme_scaled.gro...
Read 'Grunge ROck MAChoS', 1024 atoms
Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.
There are 1 chains and 0 blocks of water and 64 residues with 1024 atoms
chain #res #atoms
1 ' ' 64 1024
No occupancies in dme_scaled.gro
Opening force field file ./oplsaa.ff/atomtypes.atp
Atomtype 814
Reading residue database... (oplsaa)
Opening force field file ./oplsaa.ff/aminoacids.rtp
Residue 53
Sorting it all out...
Processing chain 1 (1024 atoms, 64 residues)
Warning: Starting residue DME1 in chain not identified as Protein/RNA/DNA.
Warning: Starting residue DME2 in chain not identified as Protein/RNA/DNA.
Warning: Starting residue DME3 in chain not identified as Protein/RNA/DNA.
Warning: Starting residue DME4 in chain not identified as Protein/RNA/DNA.
Warning: Starting residue DME5 in chain not identified as Protein/RNA/DNA.
More than 5 unidentified residues at start of chain - disabling further warnings.
Problem with chain definition, or missing terminal residues.
This chain does not appear to contain a recognized chain molecule.
If this is incorrect, you can edit residuetypes.dat to modify the behavior.
8 out of 8 lines of specbond.dat converted successfully
Checking for duplicate atoms....
Generating any missing hydrogen atoms and/or adding termini.
Now there are 64 residues with 1024 atoms
Making bonds...
Number of bonds was 960, now 960
Generating angles, dihedrals and pairs...
Before cleaning: 1344 pairs
Before cleaning: 1344 dihedrals
Keeping all generated dihedrals
Making cmap torsions...
There are 1344 dihedrals, 0 impropers, 1664 angles
1344 pairs, 960 bonds and 0 virtual sites
Total mass 5767.859 a.m.u.
Total charge -0.000 e
Writing topology
Writing coordinate file...
--------- PLEASE NOTE ------------
You have successfully generated a topology from: dme_scaled.gro.
The Oplsaa force field is used.
--------- ETON ESAELP ------------
gcq#497: "Research ! A mere excuse for idleness; it has never achieved, and will never achieve any results of the slightest value." (Benjamin Jowett, British theologian, 1817-93)
==============================================
What does atomtype 814 mean here and what does residue 53 mean ? And why can it not get the residue DME, when i have included it in the residuetypes.dat file. I sure that the dme.itp file was not read.
What can i do to get DME recognised as a residue? this is not even an "exotic" molecule, as people who do gromacs often like to always point out. I want to do these simulations for a bunch of systems that do not belong to the gromacs library. i have the forcefield files for these systems and can modify the pdb files to somehow generate the solvent boxes. However, there must be a simpler way of getting the dme.itp file to work directly, isn't it .? Could someone please explain what mistake I am making? This ordeal makes it seem to me like gcq#497 couldn't be closer to truth.
You can simplify this procedure quite a bit.
gmx insert-molecules -f dme.pdb -ci dme.pdb -nmol 63 -o box.pdb -box 4 4 4
(no wasted time editing anything, the coordinate file conforms to the right format)
Don't use editconf. Scaling is useful in only limited cases, probably monoatomic species. If the force field parameters are good, a short simulation will achieve the correct density without totally distorting your molecules.
Don't use pdb2gmx, even though it appears to have worked. You already have dme.itp so you have the topology. The purpose of pdb2gmx is to write a topology; that's work that you don't need to do. There is also no need to define DME as a protein residue.
Your topology can be written simply by hand:
==============
#include "oplsaa.ff/forcefield.itp"
#include "dme.itp"
[ system ]
DME
[ molecules ]
DME 64
==============
Done and done. But you seem to have a functional topol.top anyway, so if you go back to the unscaled coordinates, everything should work fine.
You can simplify this procedure quite a bit.
gmx insert-molecules -f dme.pdb -ci dme.pdb -nmol 63 -o box.pdb -box 4 4 4
(no wasted time editing anything, the coordinate file conforms to the right format)
Don't use editconf. Scaling is useful in only limited cases, probably monoatomic species. If the force field parameters are good, a short simulation will achieve the correct density without totally distorting your molecules.
Don't use pdb2gmx, even though it appears to have worked. You already have dme.itp so you have the topology. The purpose of pdb2gmx is to write a topology; that's work that you don't need to do. There is also no need to define DME as a protein residue.
Your topology can be written simply by hand:
==============
#include "oplsaa.ff/forcefield.itp"
#include "dme.itp"
[ system ]
DME
[ molecules ]
DME 64
==============
Done and done. But you seem to have a functional topol.top anyway, so if you go back to the unscaled coordinates, everything should work fine.
Hi, I have found this discussion because I have a very similar problem now. Have you managed to solve it? If yes, could you share your experience with me? I a box filled with 50 OLED molecules in a random way. For generating, I've tried both packmol and gmx-utility. I also have a special force-field which I have converted from CHARMM-format. The grompp seems to find my force-field which is localized in the calculation directory. But now if I follow the advice of Justin Lemkul and white a topology-file by myselft, I got the following error:
Fatal error:
number of coordinates in coordinate file (box_50.gro, 1900)
does not match topology (topol.top, 1938)
which means double-counting. One molecule contains exactly 38 atoms, and my box is filled with 50 molecules. The topology-file looks like this:
#include "./mCBP.ff/forcefield.itp"
#include "mCBP.itp"
[ system ]
Glass
[molecules]
mCBP 50
What should I change to make a working setup? I would really appreciate help from experts.
Regards,
Olga
- Badges
- Science topic