src/gmxlib/nrama.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,2013, 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 <math.h>
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "string2.h"
  46 #include "typedefs.h"
  47 #include "random.h"
  48 #include "bondf.h"
  49 #include "futil.h"
  50 #include "gmx_fatal.h"
  51 #include "nrama.h"
  52 #include "rmpbc.h"
  53
  54 static const char *pp_pat[] = { "C", "N", "CA", "C", "N" };
  55 #define NPP (sizeof(pp_pat)/sizeof(pp_pat[0]))
  56
  57 static int d_comp(const void *a,const void *b)
  58 {
  59   t_dih *da,*db;
  60
  61   da=(t_dih *)a;
  62   db=(t_dih *)b;
  63
  64   if (da->ai[1] < db->ai[1])
  65     return -1;
  66   else if (da->ai[1] == db->ai[1])
  67     return (da->ai[2] - db->ai[2]);
  68   else
  69     return 1;
  70 }
  71
  72
  73 static void calc_dihs(t_xrama *xr)
  74 {
  75   int    i,t1,t2,t3;
  76   rvec   r_ij,r_kj,r_kl,m,n;
  77   real   sign;
  78   t_dih  *dd;
  79   gmx_rmpbc_t  gpbc=NULL;
  80
  81   gpbc = gmx_rmpbc_init(xr->idef,xr->ePBC,xr->natoms,xr->box);
  82   gmx_rmpbc(gpbc,xr->natoms,xr->box,xr->x);
  83   gmx_rmpbc_done(gpbc);
  84
  85   for(i=0; (i<xr->ndih); i++) {
  86     dd=&(xr->dih[i]);
  87     dd->ang=dih_angle(xr->x[dd->ai[0]],xr->x[dd->ai[1]],
  88                       xr->x[dd->ai[2]],xr->x[dd->ai[3]],
  89                       NULL,
  90                       r_ij,r_kj,r_kl,m,n,&sign,&t1,&t2,&t3);
  91   }
  92 }
  93
  94 gmx_bool new_data(t_xrama *xr)
  95 {
  96   if (!read_next_x(xr->oenv,xr->traj,&xr->t,xr->natoms,xr->x,xr->box))
  97     return FALSE;
  98
  99   calc_dihs(xr);
 100
 101   return TRUE;
 102 }
 103
 104 static int find_atom(const char *find,char ***names,int start,int nr)
 105 {
 106   int i;
 107
 108   for(i=start; (i<nr); i++)
 109     if (strcmp(find,*names[i]) == 0)
 110       return i;
 111   return -1;
 112 }
 113
 114 static void add_xr(t_xrama *xr,int ff[5],t_atoms *atoms)
 115 {
 116   char buf[12];
 117   int i;
 118
 119   srenew(xr->dih,xr->ndih+2);
 120   for(i=0; (i<4); i++)
 121     xr->dih[xr->ndih].ai[i]=ff[i];
 122   for(i=0; (i<4); i++)
 123     xr->dih[xr->ndih+1].ai[i]=ff[i+1];
 124   xr->ndih+=2;
 125
 126   srenew(xr->pp,xr->npp+1);
 127   xr->pp[xr->npp].iphi=xr->ndih-2;
 128   xr->pp[xr->npp].ipsi=xr->ndih-1;
 129   xr->pp[xr->npp].bShow=FALSE;
 130   sprintf(buf,"%s-%d",*atoms->resinfo[atoms->atom[ff[1]].resind].name,
 131           atoms->resinfo[atoms->atom[ff[1]].resind].nr);
 132   xr->pp[xr->npp].label=strdup(buf);
 133   xr->npp++;
 134 }
 135
 136 static void get_dih(t_xrama *xr,t_atoms *atoms)
 137 {
 138   int found,ff[NPP];
 139   int i;
 140   size_t j;
 141
 142   for(i=0; (i<atoms->nr); ) {
 143     found=i;
 144     for(j=0; (j<NPP); j++) {
 145       if ((ff[j]=find_atom(pp_pat[j],atoms->atomname,found,atoms->nr)) == -1)
 146         break;
 147       found=ff[j]+1;
 148     }
 149     if (j != NPP)
 150       break;
 151     add_xr(xr,ff,atoms);
 152     i=ff[0]+1;
 153   }
 154   fprintf(stderr,"Found %d phi-psi combinations\n",xr->npp);
 155 }
 156
 157 static int search_ff(int thisff[NPP],int ndih,int **ff)
 158 {
 159   int  j,k;
 160   gmx_bool bFound=FALSE;
 161
 162   for(j=0; (j<ndih); j++) {
 163     bFound=TRUE;
 164     for(k=1; (k<=3); k++)
 165       bFound = bFound && (thisff[k]==ff[j][k]);
 166     if (bFound) {
 167       if (thisff[0] == -1)
 168         ff[j][4]=thisff[4];
 169       else
 170         ff[j][0]=thisff[0];
 171       break;
 172     }
 173   }
 174   if (!bFound) {
 175     for(k=0; (k<5); k++)
 176       ff[ndih][k]=thisff[k];
 177     ndih++;
 178   }
 179   return ndih;
 180 }
 181
 182 static void min_max(t_xrama *xr)
 183 {
 184   int ai,i,j;
 185
 186   xr->amin=xr->natoms;
 187   xr->amax=0;
 188   for(i=0; (i<xr->ndih); i++)
 189     for(j=0; (j<4); j++) {
 190       ai=xr->dih[i].ai[j];
 191       if (ai < xr->amin)
 192         xr->amin = ai;
 193       else if (ai > xr->amax)
 194         xr->amax = ai;
 195     }
 196 }
 197
 198 static void get_dih_props(t_xrama *xr,t_idef *idef,int mult)
 199 {
 200   int     i,ft,ftype,nra;
 201   t_iatom *ia;
 202   t_dih   *dd,key;
 203
 204   ia=idef->il[F_PDIHS].iatoms;
 205   for (i=0; (i<idef->il[F_PDIHS].nr); ) {
 206     ft=ia[0];
 207     ftype=idef->functype[ft];
 208     nra=interaction_function[ftype].nratoms;
 209
 210     if (ftype != F_PDIHS)
 211       gmx_incons("ftype is not a dihedral");
 212
 213     key.ai[1]=ia[2];
 214     key.ai[2]=ia[3];
 215     if ((dd = (t_dih *)bsearch(&key,xr->dih,xr->ndih,(size_t)sizeof(key),d_comp))
 216         != NULL) {
 217       dd->mult=idef->iparams[ft].pdihs.mult;
 218       dd->phi0=idef->iparams[ft].pdihs.phiA;
 219     }
 220
 221     i+=nra+1;
 222     ia+=nra+1;
 223   }
 224   /* Fill in defaults for values not in the topology */
 225   for(i=0; (i<xr->ndih); i++) {
 226     if (xr->dih[i].mult == 0) {
 227       fprintf(stderr,
 228               "Dihedral around %d,%d not found in topology. Using mult=%d\n",
 229               xr->dih[i].ai[1],xr->dih[i].ai[2],mult);
 230       xr->dih[i].mult=mult;
 231       xr->dih[i].phi0=180;
 232     }
 233   }
 234 }
 235
 236
 237
 238 t_topology *init_rama(const output_env_t oenv,const char *infile,
 239                       const char *topfile, t_xrama *xr,int mult)
 240 {
 241   t_topology *top;
 242   int    ePBC;
 243   real   t;
 244
 245   top=read_top(topfile,&xr->ePBC);
 246
 247   /*get_dih2(xr,top->idef.functype,&(top->idef.bondeds),&(top->atoms));*/
 248   get_dih(xr,&(top->atoms));
 249   get_dih_props(xr,&(top->idef),mult);
 250   xr->natoms=read_first_x(oenv,&xr->traj,infile,&t,&(xr->x),xr->box);
 251   xr->idef=&(top->idef);
 252   xr->oenv=oenv;
 253
 254   min_max(xr);
 255   calc_dihs(xr);
 256
 257   return top;
 258 }
 259