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