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