[gmx-users] Re: grompp: Atoms in the .top are not numbered consecutively
D. Ensign
densign at stanford.edu
Thu May 4 01:00:13 CEST 2006
> > > I'm having a lovely time with my old pal grompp, except for the grumpiness after
> > > replacing eight water molecules with ions in my system and discovering -- to my
> > > astonishment -- the following message:
> > >
> > > "Atoms in the .top are not numbered consecutively from 1"
> > >
> > > Well, this is interesting, seeing as how it's preventing me from advancing, so is
> > > garnering all my attention. So I have a look at my topol.top, and I discover the
> > > following:
> > >
> > > ; BEGIN
> > > ; Topology for human aldose reductase with NAD+ and inhibitor
> > >
> > > ; Include force field
> > > #include "ffamber99.itp"
> > >
> > > ; Include protein topology
> > > #include "haldr.itp"
> > >
> > > ; Include inhibitor topology
> > > #include "inh1.itp"
> > >
> > > ; Include NAD+ topology
> > > #include "nad.itp"
> > >
> > > ; Include water topology
> > > #include "ffamber_tip3p.itp"
> > >
> > > #include "my_ions.itp"
> > >
> > > [ system ]
> > > hAldR with inhibitor in water
> > >
> > > [ molecules ]
> > > ;molecule name number
> > > AldoseReductase 1
> > > Inhibitor1 1
> > > NAD 1
> > > SOL 27048
> > > ;SOL 27056
> > > Ions 1
> > > ;END
> > >
> > > Now, the odd thing is, when I comment out references to the Ions molecule and
> > > my_ions.itp, and remove the ions from the .gro file, everything works just fine!
> > > (Well, it works fine enough.) That's odd because I haven't treated the ion atoms
> > > in my_ions.itp any differently than I treated the atoms in Inhibitor1 or in NAD,
> > > other
> > > than they're after the water.
> > >
> > >
> > paste your my_ions.itp here?
> > Yang Ye
>
> As requested, kind person, here is my ions topology include file:
>
> ; BEGIN
> [ moleculetype ]
> ; Name nrexcl
> Ions 3
>
> [ atoms ]
> ; nr type resnr residue atom cgnr charge mass typeB chargeB
> massB
> 1 amber99_31 1 Na Na 1 1 22.99 ; qtot 1
> 2 amber99_31 2 Na Na 2 1 22.99 ; qtot 2
> 3 amber99_31 3 Na Na 3 1 22.99 ; qtot 1
> 4 amber99_31 4 Na Na 4 1 22.99 ; qtot 2
> 5 amber99_31 5 Na Na 5 1 22.99 ; qtot 2
> 6 amber99_31 6 Na Na 6 1 22.99 ; qtot 2
> 7 amber99_30 7 Cl Cl 7 -1 35.45 ; qtot 1
> 8 amber99_30 8 Cl Cl 8 -1 35.45 ; qtot -2
> ; END
>
> For your comparison enjoyment, here is the .itp that I used as a guide for the template
> above. This .itp was generated after I ran genion on a previous, less interesting
> system,
> then ran pdb2gmx over again which built the .itp below. Unfortunately, running pdb2gmx
> in
> the present case would be a bit of a pain, since I'm using GAFF parameters for my
> inhibitor and the cofactor, and I want to mess with the ffamber99 files as little as
> possible.
>
> ; BEGIN -- template itp file with ions that works, whereas the one above doesn't
> [ moleculetype ]
> ; Name nrexcl
> Protein_B 3
>
> [ atoms ]
> ; nr type resnr residue atom cgnr charge mass typeB chargeB
> massB
> 1 amber99_31 1 Na Na 1 1 22.99 ; qtot 1
> 2 amber99_31 2 Na Na 2 1 22.99 ; qtot 2
> 3 amber99_30 3 Cl Cl 3 -1 35.45 ; qtot 1
> 4 amber99_30 4 Cl Cl 4 -1 35.45 ; qtot 0
> 5 amber99_30 5 Cl Cl 5 -1 35.45 ; qtot -1
> 6 amber99_30 6 Cl Cl 6 -1 35.45 ; qtot -2
>
> ; Include Position restraint file
> #ifdef POSRES
> #include "posre_B.itp"
> #endif
> ; END -- template itp file that works.
>
> I'm really wondering if having the ions after the solvent is what's messing this up --
> in
> the original case, the ions come after the solvent, whereas in the second case (the
> template case) the ions come before the solvent in the gro and in the top. Which means
> I'll be spending a portion of the afternoon writing a Perl script! Woohoo!
I know everyone will be relieved I have solved my problem. For the possibility of aid to
those who come after who might run into this problem, I provide a solution.
Recall, the problem was an error while running grompp that, "Atoms in the .top are not
numbered consecutively from 1" after adding ions with genion. Before adding ions, (ie,
after adding the solvent) everything worked just fine.
The problem seemed that grompp doesn't like to find atoms after the solvent. So I put
'#include "my_ions.itp' just before the '#include water_topology.itp' in topol.top, and
ran the Perl script below (which could probably be made more efficient, but it worked
really fast on my 90,000-atom system) to place the ions before the water in the .gro
file.
Hope this helps future GROMACS users!
Dan Ensign
Stanford University
#!/usr/bin/perl -w
# Dan Ensign, 5/3/06
# Script to take a .gro file which was generated by genion, identify the
# solvent and the ion atoms (ions listed in $ions below) and place the ions
# before the solvent in the list, and renumber them accordingly.
# It is important to identify the ions by "residue name," rather than by atom
# name, because there might be organic Cl!
use strict;
######################## MODIFY ##########################
# set this variable to the target .gro file #
my $gro="foo.gro"; #
#my $gro="sys_solv_ions.gro"; #
# add to this list of ions as necessary #
my $ions="Cl Na"; #
# what the hell, one for solvent, too #
my $solvent="SOL"; #
##########################################################
unless (open GRO, $gro) { print "Sorry -- can't find gro file $gro.\n"; exit; }
my @gro=<GRO>; close GRO;
my @newgro;
#########
# step 1. Loop through the GRO file, placing residues that are neither solvent
# nor ions into @newgro.
foreach $line (@gro) {
if ($line eq $gro[0] || $line eq $gro[1]) { # header and number of atoms
push @newgro, $line;
} else {
my $residue=substr($line,5,4);
$residue=~s/\s+//; #remove whitespace
# check to make sure this isn't solvent
unless ($solvent=~m/$residue/ || $ions=~m/$residue/) {
push @newgro, $line;
}
}
}
#########
# step 2. Identify the last residue number, and the last atom number, of the
# non-solvent, non-ion parts of the gro file
my $boxvectors=pop @newgro; # take off the last line, which should be box vectors
my $lastline=$newgro[-1];
my $lastResNo=substr($lastline,0,5);
my $lastAtomNo=substr($lastline,15,5);
# print "The last residue is '$lastResNo' and the last atom is '$lastAtomNo'\n";
#########
# step 3. Loop through the @gro again, and when the ions are found add them to
# @newgro, being sure to replace their residue numbers and atom numbers
my $currentResNo=$lastResNo+1;
my $currentAtomNo=$lastAtomNo+1;
foreach $line (@gro) {
unless ($line eq $gro[0] || $line eq $gro[1]) { # weird why I have to do this
my $residue=substr($line,5,4);
$residue=~s/\s+//; # remove whitespace
# is this residue one of the ions?
if ($ions=~m/$residue/) {
my $resNo = &make5charactersLong($currentResNo);
my $atomNo = &make5charactersLong($currentAtomNo);
substr($line,0,5)=$resNo;
substr($line,15,5)=$atomNo;
push @newgro, $line;
$currentResNo++;
$currentAtomNo++;
}
}
}
#########
# step 4. Loop through the @gro a final time, and add solvent atoms as they
# are found. Increment $currentResNo after every three atoms
my $waterAtomCount;
foreach $line (@gro) {
unless ($line eq $gro[0] || $line eq $gro[1]) { # still weird
my $residue=substr($line,5,4);
$residue=~s/\s+//; #remove whitespace
# is it SOL?
if ($solvent=~m/$residue/) {
my $resNo = &make5charactersLong($currentResNo);
my $atomNo = &make5charactersLong($currentAtomNo);
substr($line,0,5)=$resNo;
substr($line,15,5)=$atomNo;
push @newgro, $line;
$currentAtomNo++;
$waterAtomCount++;
if ($waterAtomCount==3) {
$waterAtomCount=0;
$currentResNo++;
}
}
}
}
#########
# step 5. Just add back in the box vectors!
push @newgro, $boxvectors;
unlink "rs-out";
open OUT, ">> rs-out";
foreach (@newgro) { print OUT; }
print "Printed output file called 'rs-out'.";
# all computer programs should be friendly and grateful
print " Thanks, and have a nice day!\n";
##### subroutine to make residue and atom numbers take up three spaces
sub make5charactersLong {
my $string=$_[0];
while (length($string)<5) { $string=" ".$string }
return $string;
}
More information about the gromacs.org_gmx-users
mailing list