src/tools/orise.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  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39
  40 #include "typedefs.h"
  41 #include "maths.h"
  42 #include "string2.h"
  43 #include "hxprops.h"
  44 #include "gmx_fatal.h"
  45 #include "futil.h"
  46 #include "smalloc.h"
  47 #include "readev.h"
  48 #include "macros.h"
  49 #include "confio.h"
  50 #include "copyrite.h"
  51 #include "statutil.h"
  52 #include "mcprop.h"
  53 #include "orise.h"
  54
  55 void optimize(FILE *fp,int nx,t_propfunc *f,t_pinp *p)
  56 {
  57   real *fact,sf;
  58   int  i;
  59
  60   snew(fact,nx);
  61   sf=0.0;
  62   for(i=0; (i<nx); i++) {
  63     fact[i]=1.0/(i+2.0);
  64     sf+=fact[i]*fact[i];
  65   }
  66   for(i=0; (i<nx); i++)
  67     fact[i]/=sqrt(sf);
  68
  69   do_mc(fp,nx,fact,p->step,p->v0,p->tol,p->nsteps,f);
  70
  71   fprintf(fp,"MC Done\n");
  72   fprintf(fp,"Components are:\n");
  73   sf=0.0;
  74   for(i=0; (i<nx); i++) {
  75     fprintf(fp,"EV %5d:  Fac: %8.3f\n",i,fact[i]);
  76     sf+=fact[i]*fact[i];
  77   }
  78   fprintf(fp,"Norm of vector: %8.3f\n",sqrt(sf));
  79 }
  80
  81 static rvec    *xav;
  82 static rvec    **EV;
  83 static real    **evprj;
  84 static int     nca,natoms,nframes;
  85 static atom_id *ca_index;
  86
  87 real risefunc(int nx,real x[])
  88 {
  89   static rvec *xxx=NULL;
  90   real   cx,cy,cz;
  91   double rrr,rav2,rav;
  92   int    i,j,m,n,ai;
  93
  94   if (xxx == NULL)
  95     snew(xxx,natoms);
  96
  97   rav2=0,rav=0;
  98
  99   for(j=0; (j<nframes); j++) {
 100     /* Make a structure, we only have to do Z */
 101     for(i=0; (i<nca); i++) {
 102       ai=ca_index[i];
 103       /*cx=xav[ai][XX];
 104       cy=xav[ai][YY];*/
 105       cz=xav[ai][ZZ];
 106       for(n=0; (n<nx); n++) {
 107         /*cx+=EV[n][ai][XX]*x[n]*evprj[n][j];
 108         cy+=EV[n][ai][YY]*x[n]*evprj[n][j];*/
 109         cz+=EV[n][ai][ZZ]*x[n]*evprj[n][j];
 110       }
 111       /*xxx[ai][XX]=cx;
 112       xxx[ai][YY]=cy;*/
 113       xxx[ai][ZZ]=cz;
 114     }
 115     rrr   = rise(nca,ca_index,xxx);
 116     rav  += rrr;
 117     rav2 += rrr*rrr;
 118   }
 119   rav/=nframes;
 120   rrr=sqrt(rav2/nframes-rav*rav);
 121
 122   return -rrr;
 123 }
 124
 125 real radfunc(int nx,real x[])
 126 {
 127   static rvec *xxx=NULL;
 128   real   cx,cy,cz;
 129   double rrr,rav2,rav;
 130   int    i,j,m,n,ai;
 131
 132   if (xxx == NULL)
 133     snew(xxx,natoms);
 134
 135   rav2=0,rav=0;
 136
 137   for(j=0; (j<nframes); j++) {
 138     /* Make a structure, we only have to do X & Y */
 139     for(i=0; (i<nca); i++) {
 140       ai=ca_index[i];
 141       cx=xav[ai][XX];
 142       cy=xav[ai][YY];
 143       cz=xav[ai][ZZ];
 144       for(n=0; (n<nx); n++) {
 145         cx+=EV[n][ai][XX]*x[n]*evprj[n][j];
 146         cy+=EV[n][ai][YY]*x[n]*evprj[n][j];
 147         cz+=EV[n][ai][ZZ]*x[n]*evprj[n][j];
 148       }
 149       xxx[ai][XX]=cx;
 150       xxx[ai][YY]=cy;
 151       xxx[ai][ZZ]=cz;
 152     }
 153     rrr   = radius(NULL,nca,ca_index,xxx);
 154     rav  += rrr;
 155     rav2 += rrr*rrr;
 156   }
 157   rav/=nframes;
 158   rrr=sqrt(rav2/nframes-rav*rav);
 159
 160   return -rrr;
 161 }
 162
 163 void init_optim(int nx,rvec *xxav,rvec **EEV,
 164                 real **eevprj,int nnatoms,
 165                 int nnca,atom_id *cca_index,
 166                 t_pinp *p)
 167 {
 168   xav      = xxav;
 169   EV       = EEV;
 170   evprj    = eevprj;
 171   natoms   = nnatoms;
 172   nframes  = p->nframes;
 173   nca      = nnca;
 174   ca_index = cca_index;
 175 }
 176
 177 void optim_rise(int nx,rvec *xxav,rvec **EEV,
 178                 real **eevprj,int nnatoms,
 179                 int nnca,atom_id *cca_index,
 180                 t_pinp *p)
 181 {
 182   FILE *fp;
 183
 184   fp       = ffopen("rise.log","w");
 185   init_optim(nx,xxav,EEV,eevprj,nnatoms,nnca,cca_index,p);
 186   optimize(fp,nx,risefunc,p);
 187   ffclose(fp);
 188 }
 189
 190 void optim_radius(int nx,rvec *xxav,rvec **EEV,
 191                   real **eevprj,int nnatoms,
 192                   int nnca,atom_id *cca_index,
 193                   t_pinp *p)
 194 {
 195   FILE *fp;
 196
 197   fp       = ffopen("radius.log","w");
 198
 199   init_optim(nx,xxav,EEV,eevprj,nnatoms,nnca,cca_index,p);
 200   optimize(fp,nx,radfunc,p);
 201   ffclose(fp);
 202 }
 203