src/contrib/ehdata.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 "typedefs.h"
  42 #include "gromacs/utility/smalloc.h"
  43 #include "macros.h"
  44 #include "gromacs/utility/fatalerror.h"
  45 #include "random.h"
  46 #include "gromacs/fileio/strdb.h"
  47 #include "gromacs/utility/futil.h"
  48 #include "ehdata.h"
  49
  50 typedef struct {
  51   int  nener,nomega,nq;     /* Number of entries in the table      */
  52   real *ener;               /* Energy values                       */
  53   real **omega;             /* Omega values (energy dependent)     */
  54   real ***prob,***q;        /* Probability and Q                   */
  55 } t_pq_inel;
  56
  57 typedef struct {
  58   int  nener,n2Ddata;       /* Number of entries in the table      */
  59   real *ener;               /* Energy values                       */
  60   real **prob,**data;       /* Probability and data (energy loss)  */
  61 } t_p2Ddata;
  62
  63 static int realcomp(const void *a,const void *b)
  64 {
  65   real *ra,*rb;
  66
  67   ra = (real *)a;
  68   rb = (real *)b;
  69   if (ra < rb)
  70     return -1;
  71   else if (ra > rb)
  72     return 1;
  73   else
  74     return 0;
  75 }
  76
  77 static t_p2Ddata *read_p2Ddata(char *fn)
  78 {
  79   FILE    *fp;
  80   t_p2Ddata *p2Ddata;
  81   int     i,j;
  82   double  e,p,o;
  83
  84   fprintf(stdout,"Going to read %s\n",fn);
  85   fp = gmx_ffopen(fn,"r");
  86
  87   /* Allocate memory and set constants */
  88   snew(p2Ddata,1);
  89   if (fscanf(fp,"%d%d",&p2Ddata->nener,&p2Ddata->n2Ddata) != 2)
  90     gmx_fatal(FARGS,"I need two integers: nener, n in file %s",fn);
  91
  92   snew(p2Ddata->ener,p2Ddata->nener);
  93   snew(p2Ddata->prob,p2Ddata->nener);
  94   snew(p2Ddata->data,p2Ddata->nener);
  95
  96   /* Double loop to read data */
  97   for(i=0; (i<p2Ddata->nener); i++) {
  98     fprintf(stderr,"\rEnergy %d/%d",i+1,p2Ddata->nener);
  99     snew(p2Ddata->prob[i],p2Ddata->n2Ddata);
 100     snew(p2Ddata->data[i],p2Ddata->n2Ddata);
 101
 102     for(j=0; (j<p2Ddata->n2Ddata); j++) {
 103       fscanf(fp,"%lf%lf%lf",&e,&p,&o);
 104
 105       /* Consistency check */
 106       if (j==0)
 107         p2Ddata->ener[i] = e;
 108       else if (fabs(p2Ddata->ener[i]-e) > 1e-6*e)
 109         gmx_fatal(FARGS,"Inconsistent energy %f i=%d j=%d",e,i,j);
 110       p2Ddata->prob[i][j] = p;
 111       p2Ddata->data[i][j] = o;
 112     }
 113     /* There is some noise on the data, take it away by sorting,
 114      * because otherwise binary search does not work.
 115      * This is equivalent to shifting in the data slightly along the X-axis
 116      * but better than linear search with the "real" data.
 117      */
 118     qsort(p2Ddata->data[i],p2Ddata->n2Ddata,sizeof(p2Ddata->data[0][0]),
 119           realcomp);
 120   }
 121   fprintf(stderr,"\n");
 122
 123   gmx_ffclose(fp);
 124
 125   return p2Ddata;
 126 }
 127
 128 static t_pq_inel *read_pq(char *fn)
 129 {
 130   FILE      *fp;
 131   t_pq_inel *pq;
 132   int       i,j,k;
 133   double    e,p,o,t;
 134
 135   fprintf(stdout,"Going to read %s\n",fn);
 136   fp = gmx_ffopen(fn,"r");
 137
 138   /* Allocate memory and set constants */
 139   snew(pq,1);
 140   if (fscanf(fp,"%d%d%d",&pq->nener,&pq->nomega,&pq->nq) != 3)
 141     gmx_fatal(FARGS,"I need three integers: nener, nomega, nq in file %s",fn);
 142
 143   snew(pq->ener,pq->nener);
 144   snew(pq->omega,pq->nener);
 145   snew(pq->prob,pq->nener);
 146   snew(pq->q,pq->nener);
 147
 148   /* Triple loop to read data */
 149   for(i=0; (i<pq->nener); i++) {
 150     fprintf(stderr,"\rEnergy %d/%d",i+1,pq->nener);
 151     snew(pq->prob[i],pq->nomega);
 152     snew(pq->q[i],pq->nomega);
 153     snew(pq->omega[i],pq->nomega);
 154
 155     for(j=0; (j<pq->nomega); j++) {
 156       snew(pq->prob[i][j],pq->nq);
 157       snew(pq->q[i][j],pq->nq);
 158
 159       for(k=0; (k<pq->nq); k++) {
 160         fscanf(fp,"%lf%lf%lf%lf",&e,&o,&p,&t);
 161
 162         /* Consistency check */
 163         if ((j == 0) && (k == 0))
 164           pq->ener[i] = e;
 165         else if (fabs(pq->ener[i]-e) > 1e-6*e)
 166           gmx_fatal(FARGS,"Inconsistent energy %f i=%d j=%d k=%d",e,i,j,k);
 167
 168         if (k == 0)
 169           pq->omega[i][j] = o;
 170         else if (fabs(pq->omega[i][j]-o) > 1e-6*o)
 171           gmx_fatal(FARGS,"Inconsistent omega %f i=%d j=%d k=%d",o,i,j,k);
 172
 173         pq->prob[i][j][k] = p;
 174         pq->q[i][j][k] = t;
 175       }
 176     }
 177   }
 178   fprintf(stderr,"\n");
 179
 180   gmx_ffclose(fp);
 181
 182   return pq;
 183 }
 184
 185 static int my_bsearch(real val,int ndata,real data[],gmx_bool bLower)
 186 {
 187   int ilo,ihi,imed;
 188
 189   if (val < data[0])
 190     return -1;
 191   ilo = 0;
 192   ihi = ndata-1;
 193   while ((ihi - ilo) > 1) {
 194     imed = (ilo+ihi)/2;
 195     if (data[imed] > val)
 196       ihi = imed;
 197     else
 198       ilo = imed;
 199   }
 200   /* Now val should be in between data[ilo] and data[ihi] */
 201   /* Decide which one is closest */
 202   if (bLower || ((val-data[ilo]) <= (data[ihi]-val)))
 203     return ilo;
 204   else
 205     return ihi;
 206 }
 207
 208 static real interpolate2D(int nx,int ny,real *dx,real **data,
 209                           real x0,real fy,int nx0,int ny0)
 210 {
 211   real fx;
 212
 213   fx  = (x0-dx[nx0])/(dx[nx0+1]-dx[nx0]);
 214
 215   return (fx*fy*data[nx0][ny0] + fx*(1-fy)*data[nx0][ny0+1] +
 216           (1-fx)*fy*data[nx0+1][ny0] + (1-fx)*(1-fy)*data[nx0+1][ny0+1]);
 217 }
 218
 219 real get_omega(real ekin,int *seed,FILE *fp,char *fn)
 220 {
 221   static t_p2Ddata *p2Ddata = NULL;
 222   real r,ome,fx,fy;
 223   int  eindex,oindex;
 224
 225   if (p2Ddata == NULL)
 226     p2Ddata = read_p2Ddata(fn);
 227
 228   /* Get energy index by binary search */
 229   if ((eindex = my_bsearch(ekin,p2Ddata->nener,p2Ddata->ener,TRUE)) >= 0) {
 230 #ifdef DEBUG
 231     if (eindex >= p2Ddata->nener)
 232       gmx_fatal(FARGS,"eindex (%d) out of range (max %d) in get_omega",
 233                   eindex,p2Ddata->nener);
 234 #endif
 235
 236     /* Start with random number */
 237     r = rando(seed);
 238
 239     /* Do binary search in the energy table */
 240     if ((oindex = my_bsearch(r,p2Ddata->n2Ddata,p2Ddata->prob[eindex],TRUE)) >= 0) {
 241 #ifdef DEBUG
 242       if (oindex >= p2Ddata->n2Ddata)
 243         gmx_fatal(FARGS,"oindex (%d) out of range (max %d) in get_omega",
 244                     oindex,p2Ddata->n2Ddata);
 245 #endif
 246
 247       fy = ((r-p2Ddata->prob[eindex][oindex])/
 248             (p2Ddata->prob[eindex][oindex+1]-p2Ddata->prob[eindex][oindex]));
 249       ome = interpolate2D(p2Ddata->nener,p2Ddata->n2Ddata,p2Ddata->ener,
 250                           p2Ddata->data,ekin,fy,
 251                           eindex,oindex);
 252       /* ome = p2Ddata->data[eindex][oindex];*/
 253
 254       if (fp)
 255         fprintf(fp,"%8.3f  %8.3f\n",ome,r);
 256
 257       return ome;
 258     }
 259   }
 260   return 0;
 261 }
 262
 263 real get_theta_el(real ekin,int *seed,FILE *fp,char *fn)
 264 {
 265   static t_p2Ddata *p2Ddata = NULL;
 266   real r,theta;
 267   int  eindex,tindex;
 268
 269   if (p2Ddata == NULL)
 270     p2Ddata = read_p2Ddata(fn);
 271
 272   /* Get energy index by binary search */
 273   if ((eindex = my_bsearch(ekin,p2Ddata->nener,p2Ddata->ener,TRUE)) >= 0) {
 274
 275     /* Start with random number */
 276     r = rando(seed);
 277
 278     /* Do binary search in the energy table */
 279     if ((tindex = my_bsearch(r,p2Ddata->n2Ddata,p2Ddata->prob[eindex],FALSE)) >= 0) {
 280
 281       theta = p2Ddata->data[eindex][tindex];
 282
 283       if (fp)
 284         fprintf(fp,"%8.3f  %8.3f\n",theta,r);
 285
 286       return theta;
 287     }
 288   }
 289   return 0;
 290 }
 291
 292 real get_q_inel(real ekin,real omega,int *seed,FILE *fp,char *fn)
 293 {
 294   static t_pq_inel *pq = NULL;
 295   int    eindex,oindex,tindex;
 296   real   r,theta;
 297
 298   if (pq == NULL)
 299     pq = read_pq(fn);
 300
 301   /* Get energy index by binary search */
 302   if ((eindex = my_bsearch(ekin,pq->nener,pq->ener,TRUE)) >= 0) {
 303 #ifdef DEBUG
 304     if (eindex >= pq->nener)
 305       gmx_fatal(FARGS,"eindex out of range (%d >= %d)",eindex,pq->nener);
 306 #endif
 307
 308     /* Do binary search in the energy table */
 309     if ((oindex = my_bsearch(omega,pq->nomega,pq->omega[eindex],FALSE)) >= 0) {
 310 #ifdef DEBUG
 311       if (oindex >= pq->nomega)
 312         gmx_fatal(FARGS,"oindex out of range (%d >= %d)",oindex,pq->nomega);
 313 #endif
 314
 315       /* Start with random number */
 316       r = rando(seed);
 317
 318       if ((tindex = my_bsearch(r,pq->nq,pq->prob[eindex][oindex],FALSE)) >= 0) {
 319 #ifdef DEBUG
 320         if (tindex >= pq->nq)
 321           gmx_fatal(FARGS,"tindex out of range (%d >= %d)",tindex,pq->nq);
 322 #endif
 323
 324         theta = pq->q[eindex][oindex][tindex];
 325
 326         if (fp)
 327           fprintf(fp,"get_q_inel: %8.3f  %8.3f\n",theta,r);
 328
 329         return theta;
 330       }
 331     }
 332   }
 333   return 0;
 334 }
 335
 336 static int read_cross(char *fn,real **ener,real **cross,real factor)
 337 {
 338   char   **data=NULL;
 339   double e,c;
 340   int    i,j,n;
 341
 342   fprintf(stdout,"Going to read %s\n",fn);
 343   n = get_file(fn,&data);
 344
 345   /* Allocate memory */
 346   snew(*cross,n);
 347   snew(*ener,n);
 348   for(i=j=0; (i<n); i++) {
 349     if (sscanf(data[i],"%lf%lf",&e,&c) == 2) {
 350       (*ener)[j] = e;
 351       (*cross)[j] = c*factor;
 352       j++;
 353     }
 354     sfree(data[i]);
 355   }
 356   sfree(data);
 357   if (j != n)
 358     fprintf(stderr,"There were %d empty lines in file %s\n",n-j,fn);
 359
 360   return j;
 361 }
 362
 363 real cross_inel(real ekin,real rho,char *fn)
 364 {
 365   static real *ener  = NULL;
 366   static real *cross = NULL;
 367   static int  ninel;
 368   int eindex;
 369
 370   /* Read data at first call, convert A^2 to nm^2 */
 371   if (cross == NULL)
 372     ninel = read_cross(fn,&ener,&cross,rho*0.01);
 373
 374   /* Compute index with binary search */
 375   if ((eindex = my_bsearch(ekin,ninel,ener,FALSE)) >= 0) {
 376 #ifdef DEBUG
 377     if (eindex >= ninel)
 378       gmx_fatal(FARGS,"ekin = %f, ener[0] = %f, ener[%d] = %f",ekin,
 379                   ener[0],ener[ninel-1]);
 380 #endif
 381     return cross[eindex];
 382   }
 383
 384   return 0;
 385 }
 386
 387 real cross_el(real ekin,real rho,char *fn)
 388 {
 389   static real *ener  = NULL;
 390   static real *cross = NULL;
 391   static int  nel;
 392   int eindex;
 393
 394   /* Read data at first call, convert A^2 to nm^2  */
 395   if (cross == NULL)
 396     nel = read_cross(fn,&ener,&cross,rho*0.01);
 397
 398   /* Compute index with binary search */
 399   if ((eindex = my_bsearch(ekin,nel,ener,FALSE)) >= 0) {
 400 #ifdef DEBUG
 401     if (eindex >= nel)
 402       gmx_fatal(FARGS,"ekin = %f, ener[0] = %f, ener[%d] = %f",ekin,
 403                   ener[0],ener[nel-1]);
 404 #endif
 405
 406     return cross[eindex];
 407   }
 408   return 0;
 409 }
 410
 411 real band_ener(int *seed,FILE *fp,char *fn)
 412 {
 413   static real *ener  = NULL;
 414   static real *prob  = NULL;
 415   static int  nener;
 416   int  eindex;
 417   real r;
 418
 419   /* Read data at first call, misuse read_cross function */
 420   if (prob == NULL)
 421     nener = read_cross(fn,&ener,&prob,1.0);
 422
 423   r = rando(seed);
 424
 425   if ((eindex = my_bsearch(r,nener,prob,FALSE)) >= 0) {
 426 #ifdef DEBUG
 427     if (eindex >= nener)
 428       gmx_fatal(FARGS,"r = %f, prob[0] = %f, prob[%d] = %f",r,
 429                   prob[0],prob[nener-1]);
 430 #endif
 431
 432     if (fp)
 433       fprintf(fp,"%8.3f  %8.3f\n",ener[eindex],r);
 434
 435     return ener[eindex];
 436   }
 437   return 0;
 438 }
 439
 440 static void test_omega(FILE *fp,int *seed)
 441 {
 442   int i;
 443
 444   fprintf(fp,"Testing the energy loss tables\n");
 445   for(i=0; (i<1000); i++) {
 446     (void) get_omega(500*rando(seed),seed,fp,NULL);
 447   }
 448 }
 449
 450 static void test_q_inel(FILE *fp,int *seed)
 451 {
 452   int i;
 453
 454   fprintf(fp,"Testing the energy/omega dependent inelastic scattering q tables\n");
 455   for(i=0; (i<1000); i++) {
 456     (void) get_q_inel(500*rando(seed),400*rando(seed),seed,fp,NULL);
 457   }
 458 }
 459
 460 static void test_theta_el(FILE *fp,int *seed)
 461 {
 462   int i;
 463
 464   fprintf(fp,"Testing the energy dependent elastic scattering theta tables\n");
 465   for(i=0; (i<1000); i++) {
 466     (void) get_theta_el(500*rando(seed),seed,fp,NULL);
 467   }
 468 }
 469
 470 static void test_sigma_el(FILE *fp,real rho)
 471 {
 472   int  i;
 473
 474   fprintf(fp,"Testing the elastic cross sections table\n");
 475   for(i=0; (i<500); i++) {
 476     fprintf(fp,"%3d  %8.3f\n",i,cross_el(i,rho,NULL));
 477   }
 478 }
 479
 480 static void test_sigma_inel(FILE *fp,real rho)
 481 {
 482   int  i;
 483
 484   fprintf(fp,"Testing the inelastic cross sections table\n");
 485   for(i=0; (i<500); i++) {
 486     fprintf(fp,"%3d  %8.3f\n",i,cross_inel(i,rho,NULL));
 487   }
 488 }
 489
 490 static void test_band_ener(int *seed,FILE *fp)
 491 {
 492   int i;
 493
 494   for(i=0; (i<500); i++) {
 495     band_ener(seed,fp,NULL);
 496   }
 497 }
 498
 499 void init_tables(int nfile,t_filenm fnm[],real rho)
 500 {
 501   int  seed  = 1993;
 502   real ekin  = 20;
 503   real omega = 10;
 504
 505   (void) band_ener(&seed,NULL,opt2fn("-band",nfile,fnm));
 506   (void) cross_el(ekin,rho,opt2fn("-sigel",nfile,fnm));
 507   (void) cross_inel(ekin,rho,opt2fn("-sigin",nfile,fnm));
 508   (void) get_theta_el(ekin,&seed,NULL,opt2fn("-thetael",nfile,fnm));
 509   (void) get_omega(ekin,&seed,NULL,opt2fn("-eloss",nfile,fnm));
 510   (void) get_q_inel(ekin,omega,&seed,NULL,opt2fn("-qtrans",nfile,fnm));
 511 }
 512
 513 void test_tables(int *seed,char *fn,real rho)
 514 {
 515   FILE *fp;
 516
 517   fp = fopen(fn,"w");
 518
 519   test_omega(fp,seed);
 520   test_q_inel(fp,seed);
 521   test_theta_el(fp,seed);
 522   test_sigma_el(fp,rho);
 523   test_sigma_inel(fp,rho);
 524   test_band_ener(seed,fp);
 525
 526   fclose(fp);
 527 }
 528