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