Improved pdb2gmx error message

Increased the extent to which pdb2gmx makes use of the encapsulation
available with the t_atoms struct. Now the function that writes the
error message has an easier time being more descriptive of what
is wrong when the atoms in the input coordinate file do not match
those of the .rtp entry.

Change-Id: Ie39dea18efa358f90decad6491ddc04e3756ac39

diff --git a/src/kernel/gen_ad.c b/src/kernel/gen_ad.c
index 4b6bdb700cbee24adb56846d65ad6bcc444be322..a3234252f95c08574211410e51e6521fe37a383a 100644 (file)
--- a/src/kernel/gen_ad.c
+++ b/src/kernel/gen_ad.c
@@ -452,7 +452,7 @@ static int get_impropers(t_atoms *atoms,t_hackblock hb[],t_param **idih,
        bStop=FALSE;
        for(k=0; (k<4) && !bStop; k++) {
          ai[k] = search_atom(idihs->b[j].a[k],start,
-                             atoms->nr,atoms->atom,atoms->atomname,
+                          atoms,
                              "improper",bAllowMissing);
          if (ai[k] == NO_ATID)
            bStop = TRUE;
@@ -535,10 +535,10 @@ static void gen_excls(t_atoms *atoms, t_excls *excls, t_hackblock hb[],
       
       for(e=0; e<hbexcl->nb; e++) {
        anm = hbexcl->b[e].a[0];
-       i1 = search_atom(anm,astart,atoms->nr,atoms->atom,atoms->atomname,
+       i1 = search_atom(anm,astart,atoms,
                         "exclusion",bAllowMissing);
        anm = hbexcl->b[e].a[1];
-       i2 = search_atom(anm,astart,atoms->nr,atoms->atom,atoms->atomname,
+       i2 = search_atom(anm,astart,atoms,
                         "exclusion",bAllowMissing);
        if (i1!=NO_ATID && i2!=NO_ATID) {
          if (i1 > i2) {

diff --git a/src/kernel/genhydro.c b/src/kernel/genhydro.c
index 94997bfcb9db550475691108241fdd9cfb6545f6..a09693dd03f3e01d3ed3ac0ecfe0e5233b150011 100644 (file)
--- a/src/kernel/genhydro.c
+++ b/src/kernel/genhydro.c
@@ -74,7 +74,7 @@ static atom_id pdbasearch_atom(const char *name,int resind,t_atoms *pdba,
   for(i=0; (i<pdba->nr) && (pdba->atom[i].resind != resind); i++)
     ;
     
-  return search_atom(name,i,pdba->nr,pdba->atom,pdba->atomname,
+  return search_atom(name,i,pdba,
                     searchtype,bAllowMissing);
 }
 

diff --git a/src/kernel/pdb2top.c b/src/kernel/pdb2top.c
index 535c2e53f93e42820949a264e3b162ce10c0510c..98b5d13347601f7353edcff672f91b28219ce2f5 100644 (file)
--- a/src/kernel/pdb2top.c
+++ b/src/kernel/pdb2top.c
@@ -703,20 +703,23 @@ void write_top(FILE *out, char *pr,char *molname,
 }
 
 static atom_id search_res_atom(const char *type,int resind,
-                              int natom,t_atom at[],
-                              char ** const *aname,
+                   t_atoms *atoms,
                               const char *bondtype,gmx_bool bAllowMissing)
 {
   int i;
 
-  for(i=0; (i<natom); i++)
-    if (at[i].resind == resind)
-      return search_atom(type,i,natom,at,aname,bondtype,bAllowMissing);
+  for(i=0; (i<atoms->nr); i++)
+  {
+    if (atoms->atom[i].resind == resind)
+    {
+      return search_atom(type,i,atoms,bondtype,bAllowMissing);
+    }
+  }
   
   return NO_ATID;
 }
 
-static void do_ssbonds(t_params *ps,int natoms,t_atom atom[],char **aname[],
+static void do_ssbonds(t_params *ps,t_atoms *atoms,
                       int nssbonds,t_ssbond *ssbonds,gmx_bool bAllowMissing)
 {
   int     i,ri,rj;
@@ -725,9 +728,9 @@ static void do_ssbonds(t_params *ps,int natoms,t_atom atom[],char **aname[],
   for(i=0; (i<nssbonds); i++) {
     ri = ssbonds[i].res1;
     rj = ssbonds[i].res2;
-    ai = search_res_atom(ssbonds[i].a1,ri,natoms,atom,aname,
+    ai = search_res_atom(ssbonds[i].a1,ri,atoms,
                         "special bond",bAllowMissing);
-    aj = search_res_atom(ssbonds[i].a2,rj,natoms,atom,aname,
+    aj = search_res_atom(ssbonds[i].a2,rj,atoms,
                         "special bond",bAllowMissing);
     if ((ai == NO_ATID) || (aj == NO_ATID))
       gmx_fatal(FARGS,"Trying to make impossible special bond (%s-%s)!",
@@ -743,8 +746,8 @@ static gmx_bool inter_res_bond(const t_rbonded *b)
 }
 
 static void at2bonds(t_params *psb, t_hackblock *hb,
-                     int natoms, t_atom atom[], char **aname[], 
-                     int nres, rvec x[], 
+                     t_atoms *atoms,
+                     rvec x[],
                      real long_bond_dist, real short_bond_dist,
                      gmx_bool bAllowMissing)
 {
@@ -763,16 +766,16 @@ static void at2bonds(t_params *psb, t_hackblock *hb,
 
   fprintf(stderr,"Making bonds...\n");
   i=0;
-  for(resind=0; (resind < nres) && (i<natoms); resind++) {
+  for(resind=0; (resind < atoms->nres) && (i<atoms->nr); resind++) {
     /* add bonds from list of bonded interactions */
     for(j=0; j < hb[resind].rb[ebtsBONDS].nb; j++) {
       /* Unfortunately we can not issue errors or warnings
        * for missing atoms in bonds, as the hydrogens and terminal atoms
        * have not been added yet.
        */
-      ai=search_atom(hb[resind].rb[ebtsBONDS].b[j].AI,i,natoms,atom,aname,
+      ai=search_atom(hb[resind].rb[ebtsBONDS].b[j].AI,i,atoms,
                     ptr,TRUE);
-      aj=search_atom(hb[resind].rb[ebtsBONDS].b[j].AJ,i,natoms,atom,aname,
+      aj=search_atom(hb[resind].rb[ebtsBONDS].b[j].AJ,i,atoms,
                     ptr,TRUE);
       if (ai != NO_ATID && aj != NO_ATID) {
           dist2 = distance2(x[ai],x[aj]);
@@ -790,11 +793,11 @@ static void at2bonds(t_params *psb, t_hackblock *hb,
       }
     }
     /* add bonds from list of hacks (each added atom gets a bond) */
-    while( (i<natoms) && (atom[i].resind == resind) ) {
+    while( (i<atoms->nr) && (atoms->atom[i].resind == resind) ) {
       for(j=0; j < hb[resind].nhack; j++)
        if ( ( hb[resind].hack[j].tp > 0 ||
               hb[resind].hack[j].oname==NULL ) &&
-            strcmp(hb[resind].hack[j].AI,*(aname[i])) == 0 ) {
+            strcmp(hb[resind].hack[j].AI,*(atoms->atomname[i])) == 0 ) {
          switch(hb[resind].hack[j].tp) {
          case 9:          /* COOH terminus */
            add_param(psb,i,i+1,NULL,NULL);     /* C-O  */
@@ -1350,9 +1353,6 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
 {
     int residx,i,j,k;
     const char *ptr;
-    int natoms = atoms->nr;
-    t_atom *atom = atoms->atom;
-    char ***aname = atoms->atomname;
     t_resinfo *resinfo = atoms->resinfo;
     int nres = atoms->nres;
     gmx_bool bAddCMAP;
@@ -1380,7 +1380,7 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
             for(k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
             {
                 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
-                                             i,natoms,atom,aname,ptr,TRUE);
+                                             i,atoms,ptr,TRUE);
                 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
                 if (!bAddCMAP)
                 {
@@ -1389,7 +1389,7 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
                      * into the atom array. */
                     break;
                 }
-                this_residue_index = atom[cmap_atomid[k]].resind;
+                this_residue_index = atoms->atom[cmap_atomid[k]].resind;
                 if (0 == k)
                 {
                     cmap_chainnum = resinfo[this_residue_index].chainnum;
@@ -1413,7 +1413,7 @@ static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuet
                
                if(residx<nres-1)
                {
-                       while(atom[i].resind<residx+1)
+                       while(atoms->atom[i].resind<residx+1)
                        {
                                i++;
                        }
@@ -1472,12 +1472,12 @@ void pdb2top(FILE *top_file, char *posre_fn, char *molname,
   
   /* Make bonds */
   at2bonds(&(plist[F_BONDS]), hb, 
-           atoms->nr, atoms->atom, atoms->atomname, atoms->nres, *x, 
+           atoms, *x,
            long_bond_dist, short_bond_dist, bAllowMissing);
   
   /* specbonds: disulphide bonds & heme-his */
   do_ssbonds(&(plist[F_BONDS]),
-            atoms->nr, atoms->atom, atoms->atomname, nssbonds, ssbonds,
+            atoms, nssbonds, ssbonds,
             bAllowMissing);
   
   nmissat = name2type(atoms, &cgnr, atype, restp, rt);

diff --git a/src/kernel/pgutil.c b/src/kernel/pgutil.c
index 7a3e7d86818beb272b4a08cd946016726a58bdb6..e948dea471b19053e92e5ba9d5435c0d2d5ae47a 100644 (file)
--- a/src/kernel/pgutil.c
+++ b/src/kernel/pgutil.c
@@ -45,6 +45,7 @@
 #define BUFSIZE 1024
 static void atom_not_found(int fatal_errno,const char *file,int line,
                           const char *atomname,int resind,
+                           const char *resname,
                           const char *bondtype,gmx_bool bAllowMissing)
 {
     char message_buffer[BUFSIZE];
@@ -53,17 +54,21 @@ static void atom_not_found(int fatal_errno,const char *file,int line,
         if (0 != strcmp(bondtype, "atom"))
         {
             snprintf(message_buffer, 1024,
-                     "Atom %s is used in an interaction of type %s in the topology\n"
-                     "database, but an atom of that name was not found in residue\n"
-                     "number %d.\n",
-                     atomname,bondtype,resind+1);
+                     "Residue %d named %s of a molecule in the input file was mapped\n"
+                     "to an entry in the topology database, but the atom %s used in\n"
+                     "an interaction of type %s in that entry is not found in the\n"
+                     "input file. Perhaps your atom and/or residue naming needs to be\n"
+                     "fixed.\n",
+                     resind+1, resname, atomname, bondtype);
         }
         else
         {
             snprintf(message_buffer, 1024,
-                     "Atom %s is used in the topology database, but an atom of that\n"
-                     "name was not found in residue number %d.\n",
-                     atomname,resind+1);
+                     "Residue %d named %s of a molecule in the input file was mapped\n"
+                     "to an entry in the topology database, but the atom %s used in\n"
+                     "that entry is not found in the input file. Perhaps your atom\n"
+                     "and/or residue naming needs to be fixed.\n",
+                     resind+1, resname, atomname);
         }
         if (bAllowMissing)
         {
@@ -76,12 +81,15 @@ static void atom_not_found(int fatal_errno,const char *file,int line,
     }
 }
        
-atom_id search_atom(const char *type,int start,int natoms,t_atom at[],
-                   char ** const * anm,
+atom_id search_atom(const char *type,int start,
+                    t_atoms *atoms,
                    const char *bondtype,gmx_bool bAllowMissing)
 {
   int     i,resind=-1;
   gmx_bool    bPrevious,bNext;
+  int natoms = atoms->nr;
+  t_atom *at = atoms->atom;
+  char ** const * anm = atoms->atomname;
 
   bPrevious = (strchr(type,'-') != NULL);
   bNext     = (strchr(type,'+') != NULL);
@@ -102,7 +110,9 @@ atom_id search_atom(const char *type,int start,int natoms,t_atom at[],
        return (atom_id) i;
     }
     if (!(bNext && at[start].resind==at[natoms-1].resind))
-      atom_not_found(FARGS,type,at[start].resind,bondtype,bAllowMissing);
+    {
+        atom_not_found(FARGS,type,at[start].resind,*atoms->resinfo[resind].name,bondtype,bAllowMissing);
+    }
   }
   else {
     /* The previous residue */
@@ -113,7 +123,9 @@ atom_id search_atom(const char *type,int start,int natoms,t_atom at[],
       if (gmx_strcasecmp(type,*(anm[i]))==0)
        return (atom_id) i;
     if (start > 0)
-      atom_not_found(FARGS,type,at[start].resind,bondtype,bAllowMissing);
+    {
+        atom_not_found(FARGS,type,at[start].resind,*atoms->resinfo[resind].name,bondtype,bAllowMissing);
+    }
   }
   return NO_ATID;
 }

diff --git a/src/kernel/pgutil.h b/src/kernel/pgutil.h
index b4ac091894b803475b2cecd7e5b82aa8f65967df..040c310f165a2250a55d2524e2966571a247b11a 100644 (file)
--- a/src/kernel/pgutil.h
+++ b/src/kernel/pgutil.h
@@ -39,8 +39,7 @@
 #include "typedefs.h"
 
 extern atom_id search_atom(const char *type,int start,
-                          int natoms,t_atom at[],
-                          char ** const * anm,
+                          t_atoms *atoms,
                           const char *bondtype,gmx_bool bAllowMissing);
 /* Search an atom in array of pointers to strings, starting from start
  * if type starts with '-' then searches backwards from start.