src/contrib/ehdata.c

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