src/contrib/ehanal.c

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