src/gmxlib/nrama.c

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