src/gmxlib/nrama.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "sysstuff.h"
  41 #include "smalloc.h"
  42 #include "string2.h"
  43 #include "typedefs.h"
  44 #include "random.h"
  45 #include "bondf.h"
  46 #include "futil.h"
  47 #include "gmx_fatal.h"
  48 #include "nrama.h"
  49 #include "rmpbc.h"
  50
  51 static const char *pp_pat[] = { "C", "N", "CA", "C", "N" };
  52 #define NPP (sizeof(pp_pat)/sizeof(pp_pat[0]))
  53
  54 static int d_comp(const void *a,const void *b)
  55 {
  56   t_dih *da,*db;
  57
  58   da=(t_dih *)a;
  59   db=(t_dih *)b;
  60
  61   if (da->ai[1] < db->ai[1])
  62     return -1;
  63   else if (da->ai[1] == db->ai[1])
  64     return (da->ai[2] - db->ai[2]);
  65   else
  66     return 1;
  67 }
  68
  69
  70 static void calc_dihs(t_xrama *xr)
  71 {
  72   int    i,t1,t2,t3;
  73   rvec   r_ij,r_kj,r_kl,m,n;
  74   real   sign;
  75   t_dih  *dd;
  76
  77   rm_pbc(xr->idef,xr->ePBC,xr->natoms,xr->box,xr->x,xr->x);
  78
  79   for(i=0; (i<xr->ndih); i++) {
  80     dd=&(xr->dih[i]);
  81     dd->ang=dih_angle(xr->x[dd->ai[0]],xr->x[dd->ai[1]],
  82                       xr->x[dd->ai[2]],xr->x[dd->ai[3]],
  83                       NULL,
  84                       r_ij,r_kj,r_kl,m,n,&sign,&t1,&t2,&t3);
  85   }
  86 }
  87
  88 bool new_data(t_xrama *xr)
  89 {
  90   if (!read_next_x(xr->oenv,xr->traj,&xr->t,xr->natoms,xr->x,xr->box))
  91     return FALSE;
  92
  93   calc_dihs(xr);
  94
  95   return TRUE;
  96 }
  97
  98 static int find_atom(const char *find,char ***names,int start,int nr)
  99 {
 100   int i;
 101
 102   for(i=start; (i<nr); i++)
 103     if (strcmp(find,*names[i]) == 0)
 104       return i;
 105   return -1;
 106 }
 107
 108 static void add_xr(t_xrama *xr,int ff[5],t_atoms *atoms)
 109 {
 110   char buf[12];
 111   int i;
 112
 113   srenew(xr->dih,xr->ndih+2);
 114   for(i=0; (i<4); i++)
 115     xr->dih[xr->ndih].ai[i]=ff[i];
 116   for(i=0; (i<4); i++)
 117     xr->dih[xr->ndih+1].ai[i]=ff[i+1];
 118   xr->ndih+=2;
 119
 120   srenew(xr->pp,xr->npp+1);
 121   xr->pp[xr->npp].iphi=xr->ndih-2;
 122   xr->pp[xr->npp].ipsi=xr->ndih-1;
 123   xr->pp[xr->npp].bShow=FALSE;
 124   sprintf(buf,"%s-%d",*atoms->resinfo[atoms->atom[ff[1]].resind].name,
 125           atoms->resinfo[atoms->atom[ff[1]].resind].nr);
 126   xr->pp[xr->npp].label=strdup(buf);
 127   xr->npp++;
 128 }
 129
 130 static void get_dih(t_xrama *xr,t_atoms *atoms)
 131 {
 132   int found,ff[NPP];
 133   int i;
 134   size_t j;
 135
 136   for(i=0; (i<atoms->nr); ) {
 137     found=i;
 138     for(j=0; (j<NPP); j++) {
 139       if ((ff[j]=find_atom(pp_pat[j],atoms->atomname,found,atoms->nr)) == -1)
 140         break;
 141       found=ff[j]+1;
 142     }
 143     if (j != NPP)
 144       break;
 145     add_xr(xr,ff,atoms);
 146     i=ff[0]+1;
 147   }
 148   fprintf(stderr,"Found %d phi-psi combinations\n",xr->npp);
 149 }
 150
 151 static int search_ff(int thisff[NPP],int ndih,int **ff)
 152 {
 153   int  j,k;
 154   bool bFound=FALSE;
 155
 156   for(j=0; (j<ndih); j++) {
 157     bFound=TRUE;
 158     for(k=1; (k<=3); k++)
 159       bFound = bFound && (thisff[k]==ff[j][k]);
 160     if (bFound) {
 161       if (thisff[0] == -1)
 162         ff[j][4]=thisff[4];
 163       else
 164         ff[j][0]=thisff[0];
 165       break;
 166     }
 167   }
 168   if (!bFound) {
 169     for(k=0; (k<5); k++)
 170       ff[ndih][k]=thisff[k];
 171     ndih++;
 172   }
 173   return ndih;
 174 }
 175
 176 static void get_dih2(t_xrama *xr,t_functype functype[],
 177                      t_ilist *bondeds,t_atoms *atoms)
 178 {
 179   int     found,**ff,thisff[NPP];
 180   int     i,type,ftype,nat,ai,aj,ak,al;
 181   int     j,k;
 182   char    *cai,*caj,*cak,*cal;
 183   int     ndih,maxdih;
 184   t_iatom *iatoms;
 185
 186   ndih=0;
 187   maxdih=1+(bondeds->nr/5);
 188   snew(ff,maxdih);
 189   for(j=0; (j<maxdih); j++) {
 190     snew(ff[j],NPP);
 191   }
 192   fprintf(stderr,"There may be as many as %d dihedrals...\n",maxdih);
 193
 194   iatoms=bondeds->iatoms;
 195   for(i=0; (i<bondeds->nr); ) {
 196     type=iatoms[0];
 197     ftype=functype[type];
 198     nat=interaction_function[ftype].nratoms+1;
 199
 200     if (ftype == F_PDIHS) {
 201       ai=iatoms[1]; cai=*atoms->atomname[ai];
 202       aj=iatoms[2]; caj=*atoms->atomname[aj];
 203       ak=iatoms[3]; cak=*atoms->atomname[ak];
 204       al=iatoms[4]; cal=*atoms->atomname[al];
 205
 206       for(j=0; (j<NPP); j++)
 207         thisff[j]=-1;
 208       if (strcasecmp(cai,"C") == 0) {
 209         /* May be a Phi angle */
 210         if ((strcasecmp(caj,"N") == 0) &&
 211             (strcasecmp(cak,"CA") == 0) &&
 212             (strcasecmp(cal,"C") == 0))
 213           thisff[0]=ai,thisff[1]=aj,thisff[2]=ak,thisff[3]=al;
 214       }
 215       else if (strcasecmp(cai,"N") == 0) {
 216         /* May be a Psi angle */
 217         if ((strcasecmp(caj,"CA") == 0) &&
 218             (strcasecmp(cak,"C") == 0) &&
 219             (strcasecmp(cal,"N") == 0))
 220           thisff[1]=ai,thisff[2]=aj,thisff[3]=ak,thisff[4]=al;
 221       }
 222       if (thisff[1] != -1) {
 223         ndih=search_ff(thisff,ndih,ff);
 224         if (ndih > maxdih)
 225           gmx_fatal(FARGS,"More than %n dihedrals found. SEGV?",maxdih);
 226       }
 227       else {
 228         fprintf(stderr,"No Phi or Psi, atoms: %s %s %s %s\n",
 229                 cai,caj,cak,cal);
 230       }
 231     }
 232     i+=nat;
 233     iatoms+=nat;
 234   }
 235   for(j=0; (j<ndih); j++) {
 236     if ((ff[j][0] != -1) && (ff[j][4] != -1))
 237       add_xr(xr,ff[j],atoms);
 238     else {
 239       fprintf(stderr,"Incomplete dihedral(%d) atoms: ",j);
 240       for(k=0; (k<NPP); k++)
 241         fprintf(stderr,"%2s ",
 242                 ff[j][k] == -1 ? "-1" : *atoms->atomname[ff[j][k]]);
 243       fprintf(stderr,"\n");
 244     }
 245   }
 246   fprintf(stderr,"Found %d phi-psi combinations\n",xr->npp);
 247 }
 248
 249 static void min_max(t_xrama *xr)
 250 {
 251   int ai,i,j;
 252
 253   xr->amin=xr->natoms;
 254   xr->amax=0;
 255   for(i=0; (i<xr->ndih); i++)
 256     for(j=0; (j<4); j++) {
 257       ai=xr->dih[i].ai[j];
 258       if (ai < xr->amin)
 259         xr->amin = ai;
 260       else if (ai > xr->amax)
 261         xr->amax = ai;
 262     }
 263 }
 264
 265 static void get_dih_props(t_xrama *xr,t_idef *idef,int mult)
 266 {
 267   int     i,ft,ftype,nra;
 268   t_iatom *ia;
 269   t_dih   *dd,key;
 270
 271   ia=idef->il[F_PDIHS].iatoms;
 272   for (i=0; (i<idef->il[F_PDIHS].nr); ) {
 273     ft=ia[0];
 274     ftype=idef->functype[ft];
 275     nra=interaction_function[ftype].nratoms;
 276
 277     if (ftype != F_PDIHS)
 278       gmx_incons("ftype is not a dihedral");
 279
 280     key.ai[1]=ia[2];
 281     key.ai[2]=ia[3];
 282     if ((dd = (t_dih *)bsearch(&key,xr->dih,xr->ndih,(size_t)sizeof(key),d_comp))
 283         != NULL) {
 284       dd->mult=idef->iparams[ft].pdihs.mult;
 285       dd->phi0=idef->iparams[ft].pdihs.phiA;
 286     }
 287
 288     i+=nra+1;
 289     ia+=nra+1;
 290   }
 291   /* Fill in defaults for values not in the topology */
 292   for(i=0; (i<xr->ndih); i++) {
 293     if (xr->dih[i].mult == 0) {
 294       fprintf(stderr,
 295               "Dihedral around %d,%d not found in topology. Using mult=%d\n",
 296               xr->dih[i].ai[1],xr->dih[i].ai[2],mult);
 297       xr->dih[i].mult=mult;
 298       xr->dih[i].phi0=180;
 299     }
 300   }
 301 }
 302
 303
 304
 305 t_topology *init_rama(const output_env_t oenv,const char *infile,
 306                       const char *topfile, t_xrama *xr,int mult)
 307 {
 308   t_topology *top;
 309   int    ePBC;
 310   real   t;
 311
 312   top=read_top(topfile,&xr->ePBC);
 313
 314   /*get_dih2(xr,top->idef.functype,&(top->idef.bondeds),&(top->atoms));*/
 315   get_dih(xr,&(top->atoms));
 316   get_dih_props(xr,&(top->idef),mult);
 317   xr->natoms=read_first_x(oenv,&xr->traj,infile,&t,&(xr->x),xr->box);
 318   xr->idef=&(top->idef);
 319   xr->oenv=oenv;
 320
 321   min_max(xr);
 322   calc_dihs(xr);
 323
 324   return top;
 325 }
 326