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