src/contrib/ehanal.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.3.3
  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-2008, 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 Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdlib.h>
  40 #include <math.h>
  41 #include <string.h>
  42 #include "typedefs.h"
  43 #include "gromacs/utility/smalloc.h"
  44 #include "macros.h"
  45 #include "gromacs/utility/fatalerror.h"
  46 #include "random.h"
  47 #include "gromacs/fileio/pdbio.h"
  48 #include "gromacs/utility/futil.h"
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/fileio/xvgr.h"
  51 #include "gromacs/math/vec.h"
  52 #include "names.h"
  53 #include "ehdata.h"
  54 #include "gromacs/fileio/pdbio.h"
  55
  56 t_histo *init_histo(int np,real minx,real maxx)
  57 {
  58   t_histo *h;
  59
  60   snew(h,1);
  61   snew(h->y,np+1);
  62   snew(h->nh,np+1);
  63   h->np   = np;
  64   if (maxx <= minx)
  65     gmx_fatal(FARGS,"minx (%f) should be less than maxx (%f) in init_histo",minx,maxx);
  66   h->minx = minx;
  67   h->maxx = maxx;
  68   h->dx_1 = np/(maxx-minx);
  69   h->dx   = (maxx-minx)/np;
  70
  71   return h;
  72 }
  73
  74 void done_histo(t_histo *h)
  75 {
  76   sfree(h->y);
  77   sfree(h->nh);
  78   h->np = 0;
  79 }
  80
  81 void add_histo(t_histo *h,real x,real y)
  82 {
  83   int n;
  84
  85   n = (x-h->minx)*h->dx_1;
  86   if ((n < 0) || (n > h->np))
  87     gmx_fatal(FARGS,"Invalid x (%f) in add_histo. Should be in %f - %f",x,h->minx,h->maxx);
  88   h->y[n] += y;
  89   h->nh[n]++;
  90 }
  91
  92 void dump_histo(t_histo *h,char *fn,char *title,char *xaxis,char *yaxis,
  93                 int enorm,real norm_fac)
  94 {
  95   FILE *fp;
  96   int  i,nn;
  97
  98   for(nn=h->np; (nn > 0); nn--)
  99     if (h->nh[nn] != 0)
 100       break;
 101   for(i=0; (i<nn); i++)
 102     if (h->nh[i] > 0)
 103       break;
 104   fp = xvgropen(fn,title,xaxis,yaxis);
 105   for(  ; (i<nn); i++) {
 106     switch (enorm) {
 107     case enormNO:
 108       fprintf(fp,"%12f  %12f  %d\n",h->minx+h->dx*i,h->y[i],h->nh[i]);
 109       break;
 110     case enormFAC:
 111       fprintf(fp,"%12f  %12f  %d\n",h->minx+h->dx*i,h->y[i]*norm_fac,h->nh[i]);
 112       break;
 113     case enormNP:
 114       if (h->nh[i] > 0)
 115         fprintf(fp,"%12f  %12f  %d\n",
 116                 h->minx+h->dx*i,h->y[i]*norm_fac/h->nh[i],h->nh[i]);
 117       break;
 118     default:
 119       gmx_fatal(FARGS,"Wrong value for enorm (%d)",enorm);
 120     }
 121   }
 122   gmx_ffclose(fp);
 123 }
 124
 125 /*******************************************************************
 126  *
 127  * Functions to analyse and monitor scattering
 128  *
 129  *******************************************************************/
 130
 131 void add_scatter_event(t_ana_scat *scatter,rvec pos,gmx_bool bInel,
 132                        real t,real ekin)
 133 {
 134   int np = scatter->np;
 135
 136   if (np == scatter->maxp) {
 137     scatter->maxp += 32;
 138     srenew(scatter->time,scatter->maxp);
 139     srenew(scatter->ekin,scatter->maxp);
 140     srenew(scatter->bInel,scatter->maxp);
 141     srenew(scatter->pos,scatter->maxp);
 142   }
 143   scatter->time[np]  = t;
 144   scatter->bInel[np] = np;
 145   scatter->ekin[np]  = ekin;
 146   copy_rvec(pos,scatter->pos[np]);
 147   scatter->np++;
 148 }
 149
 150 void reset_ana_scat(t_ana_scat *scatter)
 151 {
 152   scatter->np = 0;
 153 }
 154
 155 void done_scatter(t_ana_scat *scatter)
 156 {
 157   scatter->np = 0;
 158   sfree(scatter->time);
 159   sfree(scatter->ekin);
 160   sfree(scatter->bInel);
 161   sfree(scatter->pos);
 162 }
 163
 164 void analyse_scatter(t_ana_scat *scatter,t_histo *hmfp)
 165 {
 166   int   i,n;
 167   rvec  dx;
 168
 169   if (scatter->np > 1) {
 170     for(i=1; (i<scatter->np); i++) {
 171       rvec_sub(scatter->pos[i],scatter->pos[i-1],dx);
 172       add_histo(hmfp,scatter->ekin[i],norm(dx));
 173     }
 174   }
 175 }
 176
 177 /*******************************************************************
 178  *
 179  * Functions to analyse structural changes
 180  *
 181  *******************************************************************/
 182
 183 t_ana_struct *init_ana_struct(int nstep,int nsave,real timestep,
 184                               int maxparticle)
 185 {
 186   t_ana_struct *anal;
 187
 188   snew(anal,1);
 189   anal->nanal = 1.2*((nstep / nsave)+1);
 190   anal->index = 0;
 191   anal->dt    = nsave*timestep;
 192   snew(anal->t,anal->nanal);
 193   snew(anal->maxdist,anal->nanal);
 194   snew(anal->d2_com,anal->nanal);
 195   snew(anal->d2_origin,anal->nanal);
 196   snew(anal->nion,anal->nanal);
 197   anal->nstruct   = 1;
 198   anal->nparticle = 1;
 199   anal->maxparticle = maxparticle;
 200   snew(anal->q,1);
 201   snew(anal->x,1);
 202   snew(anal->x[0],maxparticle);
 203
 204   return anal;
 205 }
 206
 207 void done_ana_struct(t_ana_struct *anal)
 208 {
 209   int i;
 210
 211   sfree(anal->t);
 212   sfree(anal->maxdist);
 213   sfree(anal->d2_com);
 214   sfree(anal->d2_origin);
 215   sfree(anal->nion);
 216   sfree(anal->q);
 217   for(i=0; (i<anal->nstruct); i++)
 218     sfree(anal->x[i]);
 219   sfree(anal->x);
 220 }
 221
 222 void reset_ana_struct(t_ana_struct *anal)
 223 {
 224   int i;
 225
 226   for(i=0; (i<anal->nanal); i++) {
 227     anal->t[i] = 0;
 228     anal->maxdist[i] = 0;
 229     clear_rvec(anal->d2_com[i]);
 230     clear_rvec(anal->d2_origin[i]);
 231     anal->nion[i] = 0;
 232   }
 233   anal->index = 0;
 234 }
 235
 236 void add_ana_struct(t_ana_struct *total,t_ana_struct *add)
 237 {
 238   int i,m;
 239
 240   if (total->index == 0)
 241     total->index = add->index;
 242   else if (total->index != add->index)
 243     gmx_fatal(FARGS,"Analysis incompatible (total: %d, add: %d) %s, %d",
 244                 total->index,add->index,__FILE__,__LINE__);
 245   for(i=0; (i<total->index); i++) {
 246     if (total->t[i] == 0)
 247       total->t[i] = add->t[i];
 248     else if (total->t[i] != add->t[i])
 249       gmx_fatal(FARGS,"Inconsistent times in analysis (%f-%f) %s, %d",
 250                   total->t[i],add->t[i],__FILE__,__LINE__);
 251     if (add->maxdist[i] > total->maxdist[i])
 252       total->maxdist[i]  = add->maxdist[i];
 253     for(m=0; (m<DIM); m++) {
 254       total->d2_com[i][m]    += add->d2_com[i][m]/add->nion[i];
 255       total->d2_origin[i][m] += add->d2_origin[i][m]/add->nion[i];
 256     }
 257     total->nion[i]     += add->nion[i];
 258   }
 259 }
 260
 261 static void do_add_struct(t_ana_struct *anal,int nparticle,rvec x[])
 262 {
 263   int i,j;
 264
 265   if (nparticle > anal->nparticle) {
 266     for(i=0; (i<anal->nstruct); i++) {
 267       for(j=anal->nparticle; (j<nparticle); j++)
 268         copy_rvec(x[j],anal->x[i][j]);
 269     }
 270   }
 271   anal->nparticle=nparticle;
 272   srenew(anal->x,anal->nstruct+1);
 273   snew(anal->x[anal->nstruct],anal->maxparticle);
 274   for(j=0; (j<nparticle); j++)
 275     copy_rvec(x[j],anal->x[anal->nstruct][j]);
 276   anal->nstruct++;
 277 }
 278
 279 void analyse_structure(t_ana_struct *anal,real t,rvec center,
 280                        rvec x[],int nparticle,real charge[])
 281 {
 282   int  i,j,m,nel,n=0;
 283   rvec dx,com;
 284   real dx2,dx1;
 285
 286   j = anal->index;
 287   if (j >= anal->nanal)
 288     gmx_fatal(FARGS,"Too many points in analyse_structure");
 289   anal->t[j]       = t;
 290   anal->maxdist[j] = 0;
 291   clear_rvec(com);
 292   nel = 0;
 293   for(i=0; (i<nparticle); i++) {
 294     if (charge[i] < 0) {
 295       rvec_inc(com,x[i]);
 296       nel++;
 297     }
 298   }
 299   if (nel > 0)
 300     for(m=0; (m<3); m++)
 301       com[m] /= nel;
 302   for(i=0; (i<nparticle); i++) {
 303     if (charge[i] < 0) {
 304       rvec_sub(x[i],com,dx);
 305       for(m=0; (m<DIM); m++) {
 306         anal->d2_com[j][m]    += sqr(dx[m]);
 307         anal->d2_origin[j][m] += sqr(x[i][m]);
 308       }
 309       dx2 = iprod(x[i],x[i]);
 310       dx1 = sqrt(dx2);
 311       if (dx1 > anal->maxdist[j])
 312         anal->maxdist[j] = dx1;
 313       n++;
 314     }
 315   }
 316   do_add_struct(anal,nparticle,x);
 317   anal->nion[j] = n;
 318   anal->index++;
 319 }
 320
 321 void dump_ana_struct(char *rmax,char *nion,char *gyr_com,char *gyr_origin,
 322                      t_ana_struct *anal,int nsim)
 323 {
 324   FILE *fp,*gp,*hp,*kp;
 325   int  i,j;
 326   real t,d2;
 327   char *legend[] = { "Rg", "RgX", "RgY", "RgZ" };
 328
 329   fp = xvgropen(rmax,"rmax","Time (fs)","r (nm)");
 330   gp = xvgropen(nion,"N ion","Time (fs)","N ions");
 331   hp = xvgropen(gyr_com,"Radius of gyration wrt. C.O.M.",
 332                 "Time (fs)","Rg (nm)");
 333   xvgr_legend(hp,asize(legend),legend);
 334   kp = xvgropen(gyr_origin,"Radius of gyration wrt. Origin",
 335                 "Time (fs)","Rg (nm)");
 336   xvgr_legend(kp,asize(legend),legend);
 337   for(i=0; (i<anal->index); i++) {
 338     t = 1000*anal->t[i];
 339     fprintf(fp,"%12g  %10.3f\n",t,anal->maxdist[i]);
 340     fprintf(gp,"%12g  %10.3f\n",t,(1.0*anal->nion[i])/nsim-1);
 341     d2 = anal->d2_com[i][XX] + anal->d2_com[i][YY] + anal->d2_com[i][ZZ];
 342     fprintf(hp,"%12g  %10.3f  %10.3f  %10.3f  %10.3f\n",
 343             t,sqrt(d2/nsim),
 344             sqrt(anal->d2_com[i][XX]/nsim),
 345             sqrt(anal->d2_com[i][YY]/nsim),
 346             sqrt(anal->d2_com[i][ZZ]/nsim));
 347     d2 = anal->d2_origin[i][XX] + anal->d2_origin[i][YY] + anal->d2_origin[i][ZZ];
 348     fprintf(kp,"%12g  %10.3f  %10.3f  %10.3f  %10.3f\n",
 349             t,sqrt(d2/nsim),
 350             sqrt(anal->d2_origin[i][XX]/nsim),
 351             sqrt(anal->d2_origin[i][YY]/nsim),
 352             sqrt(anal->d2_origin[i][ZZ]/nsim));
 353   }
 354   gmx_ffclose(hp);
 355   gmx_ffclose(gp);
 356   gmx_ffclose(fp);
 357   gmx_ffclose(kp);
 358 }
 359
 360 void dump_as_pdb(char *pdb,t_ana_struct *anal)
 361 {
 362   FILE *kp;
 363   int  i,j;
 364   real t;
 365
 366   kp = gmx_ffopen(pdb,"w");
 367   for(i=0; (i<anal->nstruct); i++) {
 368     t = 1000*anal->t[i];
 369     fprintf(kp,"MODEL  %d  time %g fs\n",i+1,t);
 370     for(j=0; (j<anal->nparticle); j++) {
 371       fprintf(kp,get_pdbformat(),"ATOM",i+1,(j < anal->nion[i]) ? "O" : "N",
 372               "PLS",' ',1,
 373               anal->x[i][j][XX]/100,
 374               anal->x[i][j][YY]/100,
 375               anal->x[i][j][ZZ]/100);
 376       fprintf(kp,"\n");
 377     }
 378     fprintf(kp,"ENDMDL\n");
 379   }
 380   gmx_ffclose(kp);
 381 }
 382
 383 char *enms[eNR] = {
 384   "Coulomb", "Repulsion", "Potential",
 385   "EkHole",  "EkElectron", "EkLattice", "Kinetic",
 386   "Total"
 387 };
 388
 389 void add_ana_ener(t_ana_ener *ae,int nn,real e[])
 390 {
 391   int i;
 392
 393   /* First time around we are constantly increasing the array size */
 394   if (nn >= ae->nx) {
 395     if (ae->nx == ae->maxx) {
 396       ae->maxx += 1024;
 397       srenew(ae->e,ae->maxx);
 398     }
 399     for(i=0; (i<eNR); i++)
 400       ae->e[ae->nx][i] = e[i];
 401     ae->nx++;
 402   }
 403   else {
 404     for(i=0; (i<eNR); i++)
 405       ae->e[nn][i] += e[i];
 406   }
 407 }
 408
 409 void dump_ana_ener(t_ana_ener *ae,int nsim,real dt,char *edump,
 410                    t_ana_struct *total)
 411 {
 412   FILE *fp;
 413   int  i,j;
 414   real fac;
 415
 416   fac = 1.0/(nsim*ELECTRONVOLT);
 417   fp=xvgropen(edump,"Energies","Time (fs)","E (eV)");
 418   xvgr_legend(fp,eNR,enms);
 419   fprintf(fp,"@ s%d legend \"Ek/Nelec\"\n",eNR);
 420   fprintf(fp,"@ type nxy\n");
 421   for(i=0; (i<ae->nx); i++) {
 422     fprintf(fp,"%10f",1000.0*dt*i);
 423     for(j=0; (j<eNR); j++)
 424       fprintf(fp,"  %8.3f",ae->e[i][j]*fac);
 425     fprintf(fp,"  %8.3f\n",ae->e[i][eELECTRON]/(ELECTRONVOLT*total->nion[i]));
 426   }
 427   fprintf(fp,"&\n");
 428   gmx_ffclose(fp);
 429 }
 430