src/tools/orise.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_orise_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "maths.h"
  33 #include "string2.h"
  34 #include "hxprops.h"
  35 #include "fatal.h"
  36 #include "futil.h"
  37 #include "smalloc.h"
  38 #include "readev.h"
  39 #include "macros.h"
  40 #include "confio.h"
  41 #include "copyrite.h"
  42 #include "statutil.h"
  43 #include "mcprop.h"
  44 #include "orise.h"
  45
  46 void optimize(FILE *fp,int nx,t_propfunc *f,t_pinp *p)
  47 {
  48   real *fact,sf;
  49   int  i;
  50
  51   snew(fact,nx);
  52   sf=0.0;
  53   for(i=0; (i<nx); i++) {
  54     fact[i]=1.0/(i+2.0);
  55     sf+=fact[i]*fact[i];
  56   }
  57   for(i=0; (i<nx); i++)
  58     fact[i]/=sqrt(sf);
  59
  60   do_mc(fp,nx,fact,p->step,p->v0,p->tol,p->nsteps,f);
  61
  62   fprintf(fp,"MC Done\n");
  63   fprintf(fp,"Components are:\n");
  64   sf=0.0;
  65   for(i=0; (i<nx); i++) {
  66     fprintf(fp,"EV %5d:  Fac: %8.3f\n",i,fact[i]);
  67     sf+=fact[i]*fact[i];
  68   }
  69   fprintf(fp,"Norm of vector: %8.3f\n",sqrt(sf));
  70 }
  71
  72 static rvec    *xav;
  73 static rvec    **EV;
  74 static real    **evprj;
  75 static int     nca,natoms,nframes;
  76 static atom_id *ca_index;
  77
  78 real risefunc(int nx,real x[])
  79 {
  80   static rvec *xxx=NULL;
  81   real   cx,cy,cz;
  82   double rrr,rav2,rav;
  83   int    i,j,m,n,ai;
  84
  85   if (xxx == NULL)
  86     snew(xxx,natoms);
  87
  88   rav2=0,rav=0;
  89
  90   for(j=0; (j<nframes); j++) {
  91     /* Make a structure, we only have to do Z */
  92     for(i=0; (i<nca); i++) {
  93       ai=ca_index[i];
  94       /*cx=xav[ai][XX];
  95       cy=xav[ai][YY];*/
  96       cz=xav[ai][ZZ];
  97       for(n=0; (n<nx); n++) {
  98         /*cx+=EV[n][ai][XX]*x[n]*evprj[n][j];
  99         cy+=EV[n][ai][YY]*x[n]*evprj[n][j];*/
 100         cz+=EV[n][ai][ZZ]*x[n]*evprj[n][j];
 101       }
 102       /*xxx[ai][XX]=cx;
 103       xxx[ai][YY]=cy;*/
 104       xxx[ai][ZZ]=cz;
 105     }
 106     rrr   = rise(nca,ca_index,xxx);
 107     rav  += rrr;
 108     rav2 += rrr*rrr;
 109   }
 110   rav/=nframes;
 111   rrr=sqrt(rav2/nframes-rav*rav);
 112
 113   return -rrr;
 114 }
 115
 116 real radfunc(int nx,real x[])
 117 {
 118   static rvec *xxx=NULL;
 119   real   cx,cy,cz;
 120   double rrr,rav2,rav;
 121   int    i,j,m,n,ai;
 122
 123   if (xxx == NULL)
 124     snew(xxx,natoms);
 125
 126   rav2=0,rav=0;
 127
 128   for(j=0; (j<nframes); j++) {
 129     /* Make a structure, we only have to do X & Y */
 130     for(i=0; (i<nca); i++) {
 131       ai=ca_index[i];
 132       cx=xav[ai][XX];
 133       cy=xav[ai][YY];
 134       cz=xav[ai][ZZ];
 135       for(n=0; (n<nx); n++) {
 136         cx+=EV[n][ai][XX]*x[n]*evprj[n][j];
 137         cy+=EV[n][ai][YY]*x[n]*evprj[n][j];
 138         cz+=EV[n][ai][ZZ]*x[n]*evprj[n][j];
 139       }
 140       xxx[ai][XX]=cx;
 141       xxx[ai][YY]=cy;
 142       xxx[ai][ZZ]=cz;
 143     }
 144     rrr   = radius(NULL,nca,ca_index,xxx);
 145     rav  += rrr;
 146     rav2 += rrr*rrr;
 147   }
 148   rav/=nframes;
 149   rrr=sqrt(rav2/nframes-rav*rav);
 150
 151   return -rrr;
 152 }
 153
 154 void init_optim(int nx,rvec *xxav,rvec **EEV,
 155                 real **eevprj,int nnatoms,
 156                 int nnca,atom_id *cca_index,
 157                 t_pinp *p)
 158 {
 159   xav      = xxav;
 160   EV       = EEV;
 161   evprj    = eevprj;
 162   natoms   = nnatoms;
 163   nframes  = p->nframes;
 164   nca      = nnca;
 165   ca_index = cca_index;
 166 }
 167
 168 void optim_rise(int nx,rvec *xxav,rvec **EEV,
 169                 real **eevprj,int nnatoms,
 170                 int nnca,atom_id *cca_index,
 171                 t_pinp *p)
 172 {
 173   FILE *fp;
 174
 175   fp       = ffopen("rise.log","w");
 176   init_optim(nx,xxav,EEV,eevprj,nnatoms,nnca,cca_index,p);
 177   optimize(fp,nx,risefunc,p);
 178   fclose(fp);
 179 }
 180
 181 void optim_radius(int nx,rvec *xxav,rvec **EEV,
 182                   real **eevprj,int nnatoms,
 183                   int nnca,atom_id *cca_index,
 184                   t_pinp *p)
 185 {
 186   FILE *fp;
 187
 188   fp       = ffopen("radius.log","w");
 189
 190   init_optim(nx,xxav,EEV,eevprj,nnatoms,nnca,cca_index,p);
 191   optimize(fp,nx,radfunc,p);
 192   fclose(fp);
 193 }
 194