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 get_dih2(t_xrama *xr,t_functype functype[],
 180                      t_ilist *bondeds,t_atoms *atoms)
 181 {
 182   int     found,**ff,thisff[NPP];
 183   int     i,type,ftype,nat,ai,aj,ak,al;
 184   int     j,k;
 185   char    *cai,*caj,*cak,*cal;
 186   int     ndih,maxdih;
 187   t_iatom *iatoms;
 188
 189   ndih=0;
 190   maxdih=1+(bondeds->nr/5);
 191   snew(ff,maxdih);
 192   for(j=0; (j<maxdih); j++) {
 193     snew(ff[j],NPP);
 194   }
 195   fprintf(stderr,"There may be as many as %d dihedrals...\n",maxdih);
 196
 197   iatoms=bondeds->iatoms;
 198   for(i=0; (i<bondeds->nr); ) {
 199     type=iatoms[0];
 200     ftype=functype[type];
 201     nat=interaction_function[ftype].nratoms+1;
 202
 203     if (ftype == F_PDIHS) {
 204       ai=iatoms[1]; cai=*atoms->atomname[ai];
 205       aj=iatoms[2]; caj=*atoms->atomname[aj];
 206       ak=iatoms[3]; cak=*atoms->atomname[ak];
 207       al=iatoms[4]; cal=*atoms->atomname[al];
 208
 209       for(j=0; (j<NPP); j++)
 210         thisff[j]=-1;
 211       if (gmx_strcasecmp(cai,"C") == 0) {
 212         /* May be a Phi angle */
 213         if ((gmx_strcasecmp(caj,"N") == 0) &&
 214             (gmx_strcasecmp(cak,"CA") == 0) &&
 215             (gmx_strcasecmp(cal,"C") == 0))
 216           thisff[0]=ai,thisff[1]=aj,thisff[2]=ak,thisff[3]=al;
 217       }
 218       else if (gmx_strcasecmp(cai,"N") == 0) {
 219         /* May be a Psi angle */
 220         if ((gmx_strcasecmp(caj,"CA") == 0) &&
 221             (gmx_strcasecmp(cak,"C") == 0) &&
 222             (gmx_strcasecmp(cal,"N") == 0))
 223           thisff[1]=ai,thisff[2]=aj,thisff[3]=ak,thisff[4]=al;
 224       }
 225       if (thisff[1] != -1) {
 226         ndih=search_ff(thisff,ndih,ff);
 227         if (ndih > maxdih)
 228           gmx_fatal(FARGS,"More than %n dihedrals found. SEGV?",maxdih);
 229       }
 230       else {
 231         fprintf(stderr,"No Phi or Psi, atoms: %s %s %s %s\n",
 232                 cai,caj,cak,cal);
 233       }
 234     }
 235     i+=nat;
 236     iatoms+=nat;
 237   }
 238   for(j=0; (j<ndih); j++) {
 239     if ((ff[j][0] != -1) && (ff[j][4] != -1))
 240       add_xr(xr,ff[j],atoms);
 241     else {
 242       fprintf(stderr,"Incomplete dihedral(%d) atoms: ",j);
 243       for(k=0; (k<NPP); k++)
 244         fprintf(stderr,"%2s ",
 245                 ff[j][k] == -1 ? "-1" : *atoms->atomname[ff[j][k]]);
 246       fprintf(stderr,"\n");
 247     }
 248   }
 249   fprintf(stderr,"Found %d phi-psi combinations\n",xr->npp);
 250 }
 251
 252 static void min_max(t_xrama *xr)
 253 {
 254   int ai,i,j;
 255
 256   xr->amin=xr->natoms;
 257   xr->amax=0;
 258   for(i=0; (i<xr->ndih); i++)
 259     for(j=0; (j<4); j++) {
 260       ai=xr->dih[i].ai[j];
 261       if (ai < xr->amin)
 262         xr->amin = ai;
 263       else if (ai > xr->amax)
 264         xr->amax = ai;
 265     }
 266 }
 267
 268 static void get_dih_props(t_xrama *xr,t_idef *idef,int mult)
 269 {
 270   int     i,ft,ftype,nra;
 271   t_iatom *ia;
 272   t_dih   *dd,key;
 273
 274   ia=idef->il[F_PDIHS].iatoms;
 275   for (i=0; (i<idef->il[F_PDIHS].nr); ) {
 276     ft=ia[0];
 277     ftype=idef->functype[ft];
 278     nra=interaction_function[ftype].nratoms;
 279
 280     if (ftype != F_PDIHS)
 281       gmx_incons("ftype is not a dihedral");
 282
 283     key.ai[1]=ia[2];
 284     key.ai[2]=ia[3];
 285     if ((dd = (t_dih *)bsearch(&key,xr->dih,xr->ndih,(size_t)sizeof(key),d_comp))
 286         != NULL) {
 287       dd->mult=idef->iparams[ft].pdihs.mult;
 288       dd->phi0=idef->iparams[ft].pdihs.phiA;
 289     }
 290
 291     i+=nra+1;
 292     ia+=nra+1;
 293   }
 294   /* Fill in defaults for values not in the topology */
 295   for(i=0; (i<xr->ndih); i++) {
 296     if (xr->dih[i].mult == 0) {
 297       fprintf(stderr,
 298               "Dihedral around %d,%d not found in topology. Using mult=%d\n",
 299               xr->dih[i].ai[1],xr->dih[i].ai[2],mult);
 300       xr->dih[i].mult=mult;
 301       xr->dih[i].phi0=180;
 302     }
 303   }
 304 }
 305
 306
 307
 308 t_topology *init_rama(const output_env_t oenv,const char *infile,
 309                       const char *topfile, t_xrama *xr,int mult)
 310 {
 311   t_topology *top;
 312   int    ePBC;
 313   real   t;
 314
 315   top=read_top(topfile,&xr->ePBC);
 316
 317   /*get_dih2(xr,top->idef.functype,&(top->idef.bondeds),&(top->atoms));*/
 318   get_dih(xr,&(top->atoms));
 319   get_dih_props(xr,&(top->idef),mult);
 320   xr->natoms=read_first_x(oenv,&xr->traj,infile,&t,&(xr->x),xr->box);
 321   xr->idef=&(top->idef);
 322   xr->oenv=oenv;
 323
 324   min_max(xr);
 325   calc_dihs(xr);
 326
 327   return top;
 328 }
 329