src/tools/dlist.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdlib.h>
  43
  44 #include "string2.h"
  45 #include "smalloc.h"
  46 #include "gstat.h"
  47 #include "gmx_fatal.h"
  48 #include "index.h"
  49
  50 t_dlist *mk_dlist(FILE *log,
  51                   t_atoms *atoms, int *nlist,
  52                   gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bHChi,
  53                   int maxchi, int r0, gmx_residuetype_t rt)
  54 {
  55   int     ires,i,j,k,ii;
  56   t_dihatms atm,prev;
  57   int     nl=0,nc[edMax];
  58   char    *thisres;
  59   t_dlist *dl;
  60
  61   snew(dl,atoms->nres+1);
  62   prev.C = prev.O = -1;
  63   for(i=0; (i<edMax); i++)
  64     nc[i]=0;
  65   ires = -1;
  66   i    =  0;
  67   while (i<atoms->nr) {
  68     ires = atoms->atom[i].resind;
  69
  70     /* Initiate all atom numbers to -1 */
  71     atm.minC=atm.H=atm.N=atm.C=atm.O=atm.minO=-1;
  72     for(j=0; (j<MAXCHI+3); j++)
  73       atm.Cn[j]=-1;
  74
  75     /* Look for atoms in this residue */
  76     /* maybe should allow for chis to hydrogens? */
  77     while ((i<atoms->nr) && (atoms->atom[i].resind == ires)) {
  78       if ((strcmp(*(atoms->atomname[i]),"H") == 0) ||
  79           (strcmp(*(atoms->atomname[i]),"H1") == 0) )
  80         atm.H=i;
  81       else if (strcmp(*(atoms->atomname[i]),"N") == 0)
  82         atm.N=i;
  83       else if (strcmp(*(atoms->atomname[i]),"C") == 0)
  84         atm.C=i;
  85       else if ((strcmp(*(atoms->atomname[i]),"O") == 0) ||
  86                (strcmp(*(atoms->atomname[i]),"O1") == 0))
  87         atm.O=i;
  88       else if (strcmp(*(atoms->atomname[i]),"CA") == 0)
  89         atm.Cn[1]=i;
  90       else if (strcmp(*(atoms->atomname[i]),"CB") == 0)
  91         atm.Cn[2]=i;
  92       else if ((strcmp(*(atoms->atomname[i]),"CG") == 0)  ||
  93                (strcmp(*(atoms->atomname[i]),"CG1") == 0) ||
  94                (strcmp(*(atoms->atomname[i]),"OG") == 0)  ||
  95                (strcmp(*(atoms->atomname[i]),"OG1") == 0) ||
  96                (strcmp(*(atoms->atomname[i]),"SG") == 0))
  97         atm.Cn[3]=i;
  98       else if ((strcmp(*(atoms->atomname[i]),"CD") == 0)  ||
  99                (strcmp(*(atoms->atomname[i]),"CD1") == 0) ||
 100                (strcmp(*(atoms->atomname[i]),"SD") == 0)  ||
 101                (strcmp(*(atoms->atomname[i]),"OD1") == 0) ||
 102                (strcmp(*(atoms->atomname[i]),"ND1") == 0))
 103         atm.Cn[4]=i;
 104       /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
 105       else if (bHChi && ((strcmp(*(atoms->atomname[i]),"HG")  == 0) ||
 106                (strcmp(*(atoms->atomname[i]),"HG1")  == 0)) )
 107         atm.Cn[4]=i;
 108       else if ((strcmp(*(atoms->atomname[i]),"CE") == 0) ||
 109                (strcmp(*(atoms->atomname[i]),"CE1") == 0) ||
 110                (strcmp(*(atoms->atomname[i]),"OE1") == 0) ||
 111                (strcmp(*(atoms->atomname[i]),"NE") == 0))
 112         atm.Cn[5]=i;
 113       else if ((strcmp(*(atoms->atomname[i]),"CZ") == 0) ||
 114                (strcmp(*(atoms->atomname[i]),"NZ") == 0))
 115         atm.Cn[6]=i;
 116       /* HChi flag here too */
 117       else if (bHChi && (strcmp(*(atoms->atomname[i]),"NH1") == 0))
 118         atm.Cn[7]=i;
 119       i++;
 120     }
 121
 122     thisres = *(atoms->resinfo[ires].name);
 123
 124     /* added by grs - special case for aromatics, whose chis above 2 are
 125        not real and produce rubbish output - so set back to -1 */
 126     if (strcmp(thisres,"PHE") == 0 ||
 127         strcmp(thisres,"TYR") == 0 ||
 128         strcmp(thisres,"PTR") == 0 ||
 129         strcmp(thisres,"TRP") == 0 ||
 130         strcmp(thisres,"HIS") == 0 ||
 131         strcmp(thisres,"HISA") == 0 ||
 132         strcmp(thisres,"HISB") == 0 )  {
 133       for (ii=5 ; ii<=7 ; ii++)
 134         atm.Cn[ii]=-1;
 135     }
 136     /* end fixing aromatics */
 137
 138     /* Special case for Pro, has no H */
 139     if (strcmp(thisres,"PRO") == 0)
 140       atm.H=atm.Cn[4];
 141     /* Carbon from previous residue */
 142     if (prev.C != -1)
 143       atm.minC = prev.C;
 144     if (prev.O != -1)
 145       atm.minO = prev.O;
 146     prev = atm;
 147
 148     /* Check how many dihedrals we have */
 149     if ((atm.N != -1) && (atm.Cn[1] != -1) && (atm.C != -1) &&
 150       (atm.O != -1) && ((atm.H != -1) || (atm.minC != -1))) {
 151       dl[nl].resnr     = ires+1;
 152       dl[nl].atm       = atm;
 153       dl[nl].atm.Cn[0] = atm.N;
 154       if ((atm.Cn[3] != -1) && (atm.Cn[2] != -1) && (atm.Cn[1] != -1)) {
 155         nc[0]++;
 156         if (atm.Cn[4] != -1) {
 157           nc[1]++;
 158           if (atm.Cn[5] != -1) {
 159             nc[2]++;
 160             if (atm.Cn[6] != -1) {
 161               nc[3]++;
 162               if (atm.Cn[7] != -1) {
 163                 nc[4]++;
 164                 if (atm.Cn[8] != -1) {
 165                   nc[5]++;
 166                 }
 167               }
 168             }
 169           }
 170         }
 171       }
 172       if ((atm.minC != -1) && (atm.minO != -1))
 173         nc[6]++;
 174       dl[nl].index=gmx_residuetype_get_index(rt,thisres);
 175
 176       sprintf(dl[nl].name,"%s%d",thisres,ires+r0);
 177       nl++;
 178     }
 179     else if (debug)
 180       fprintf(debug,"Could not find N atom but could find other atoms"
 181               " in residue %s%d\n",thisres,ires+r0);
 182   }
 183   fprintf(stderr,"\n");
 184   fprintf(log,"\n");
 185   fprintf(log,"There are %d residues with dihedrals\n",nl);
 186   j=0;
 187   if (bPhi) j+=nl;
 188   if (bPsi) j+=nl;
 189   if (bChi)
 190     for(i=0; (i<maxchi); i++)
 191       j+=nc[i];
 192   fprintf(log,"There are %d dihedrals\n",j);
 193   fprintf(log,"Dihedral: ");
 194   if (bPhi)
 195     fprintf(log," Phi  ");
 196   if (bPsi)
 197     fprintf(log," Psi  ");
 198   if (bChi)
 199     for(i=0; (i<maxchi); i++)
 200       fprintf(log,"Chi%d  ",i+1);
 201   fprintf(log,"\nNumber:   ");
 202   if (bPhi)
 203     fprintf(log,"%4d  ",nl);
 204   if (bPsi)
 205     fprintf(log,"%4d  ",nl);
 206   if (bChi)
 207     for(i=0; (i<maxchi); i++)
 208       fprintf(log,"%4d  ",nc[i]);
 209   fprintf(log,"\n");
 210
 211   *nlist=nl;
 212
 213   return dl;
 214 }
 215
 216 gmx_bool has_dihedral(int Dih,t_dlist *dl)
 217 {
 218   gmx_bool b = FALSE;
 219   int  ddd;
 220
 221   switch (Dih) {
 222   case edPhi:
 223     b = ((dl->atm.H!=-1) && (dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1) && (dl->atm.C!=-1));
 224     break;
 225   case edPsi:
 226     b = ((dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1) && (dl->atm.C!=-1) && (dl->atm.O!=-1));
 227     break;
 228   case edOmega:
 229     b = ((dl->atm.minO!=-1) && (dl->atm.minC!=-1) && (dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1));
 230     break;
 231   case edChi1:
 232   case edChi2:
 233   case edChi3:
 234   case edChi4:
 235   case edChi5:
 236   case edChi6:
 237     ddd = Dih - edChi1;
 238     b   = ((dl->atm.Cn[ddd]!=-1) &&  (dl->atm.Cn[ddd+1]!=-1)&&
 239            (dl->atm.Cn[ddd+2]!=-1) && (dl->atm.Cn[ddd+3]!=-1));
 240     break;
 241   default:
 242     pr_dlist(stdout,1,dl,1,0,TRUE,TRUE,TRUE,TRUE,MAXCHI);
 243     gmx_fatal(FARGS,"Non existant dihedral %d in file %s, line %d",
 244                 Dih,__FILE__,__LINE__);
 245   }
 246   return b;
 247 }
 248
 249 static void pr_one_ro(FILE *fp,t_dlist *dl,int nDih,real dt)
 250 {
 251   int k ;
 252   for(k=0;k<NROT;k++)
 253     fprintf(fp,"  %6.2f",dl->rot_occ[nDih][k]);
 254   fprintf(fp,"\n");
 255 }
 256
 257 static void pr_ntr_s2(FILE *fp,t_dlist *dl,int nDih,real dt)
 258 {
 259   fprintf(fp,"  %6.2f  %6.2f\n",(dt == 0) ? 0 : dl->ntr[nDih]/dt,dl->S2[nDih]);
 260 }
 261
 262 void pr_dlist(FILE *fp,int nl,t_dlist dl[],real dt, int printtype,
 263 gmx_bool bPhi, gmx_bool bPsi,gmx_bool bChi,gmx_bool bOmega, int maxchi)
 264 {
 265   int i, Xi;
 266
 267   void  (*pr_props)(FILE *, t_dlist *, int, real);
 268
 269   /* Analysis of dihedral transitions etc */
 270
 271   if (printtype == edPrintST){
 272     pr_props=pr_ntr_s2 ;
 273     fprintf(stderr,"Now printing out transitions and OPs...\n");
 274   }else{
 275     pr_props=pr_one_ro ;
 276     fprintf(stderr,"Now printing out rotamer occupancies...\n");
 277     fprintf(fp,"\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n\n");
 278   }
 279
 280   /* change atom numbers from 0 based to 1 based */
 281   for(i=0; (i<nl); i++) {
 282     fprintf(fp,"Residue %s\n",dl[i].name);
 283     if (printtype == edPrintST){
 284       fprintf(fp," Angle [   AI,   AJ,   AK,   AL]  #tr/ns  S^2D  \n"
 285                  "--------------------------------------------\n");
 286     } else {
 287       fprintf(fp," Angle [   AI,   AJ,   AK,   AL]  rotamers  0  g(-)  t  g(+)\n"
 288                  "--------------------------------------------\n");
 289     }
 290     if (bPhi) {
 291       fprintf(fp,"   Phi [%5d,%5d,%5d,%5d]",
 292             (dl[i].atm.H == -1) ? 1+dl[i].atm.minC : 1+dl[i].atm.H,
 293             1+dl[i].atm.N, 1+dl[i].atm.Cn[1], 1+dl[i].atm.C);
 294       pr_props(fp,&dl[i],edPhi,dt);
 295     }
 296     if (bPsi) {
 297       fprintf(fp,"   Psi [%5d,%5d,%5d,%5d]",1+dl[i].atm.N, 1+dl[i].atm.Cn[1],
 298             1+dl[i].atm.C, 1+dl[i].atm.O);
 299       pr_props(fp,&dl[i],edPsi,dt);
 300     }
 301     if (bOmega && has_dihedral(edOmega,&(dl[i]))) {
 302       fprintf(fp," Omega [%5d,%5d,%5d,%5d]",1+dl[i].atm.minO, 1+dl[i].atm.minC,
 303               1+dl[i].atm.N, 1+dl[i].atm.Cn[1]);
 304       pr_props(fp,&dl[i],edOmega,dt);
 305     }
 306     for(Xi=0; Xi<MAXCHI; Xi++)
 307       if (bChi && (Xi < maxchi) && (dl[i].atm.Cn[Xi+3] != -1) ) {
 308         fprintf(fp,"   Chi%d[%5d,%5d,%5d,%5d]",Xi+1, 1+dl[i].atm.Cn[Xi],
 309                 1+dl[i].atm.Cn[Xi+1], 1+dl[i].atm.Cn[Xi+2],
 310                 1+dl[i].atm.Cn[Xi+3]);
 311         pr_props(fp,&dl[i],Xi+edChi1,dt); /* Xi+2 was wrong here */
 312       }
 313     fprintf(fp,"\n");
 314   }
 315 }
 316
 317
 318
 319 int pr_trans(FILE *fp,int nl,t_dlist dl[],real dt,int Xi)
 320 {
 321   /* never called at the moment */
 322
 323   int  i,nn,nz;
 324
 325   nz=0;
 326   fprintf(fp,"\\begin{table}[h]\n");
 327   fprintf(fp,"\\caption{Number of dihedral transitions per nanosecond}\n");
 328   fprintf(fp,"\\begin{tabular}{|l|l|}\n");
 329   fprintf(fp,"\\hline\n");
 330   fprintf(fp,"Residue\t&$\\chi_%d$\t\\\\\n",Xi+1);
 331   for(i=0; (i<nl); i++) {
 332     nn=dl[i].ntr[Xi]/dt;
 333
 334     if (nn == 0) {
 335       fprintf(fp,"%s\t&\\HL{%d}\t\\\\\n",dl[i].name,nn);
 336       nz++;
 337     }
 338     else if (nn > 0)
 339       fprintf(fp,"%s\t&\\%d\t\\\\\n",dl[i].name,nn);
 340   }
 341   fprintf(fp,"\\hline\n");
 342   fprintf(fp,"\\end{tabular}\n");
 343   fprintf(fp,"\\end{table}\n\n");
 344
 345   return nz;
 346 }
 347