src/contrib/ehdata.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2008, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdlib.h>
  43 #include <math.h>
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "copyrite.h"
  48 #include "statutil.h"
  49 #include "gmx_fatal.h"
  50 #include "random.h"
  51 #include "strdb.h"
  52 #include "futil.h"
  53 #include "physics.h"
  54 #include "ehdata.h"
  55
  56 typedef struct {
  57   int  nener,nomega,nq;     /* Number of entries in the table      */
  58   real *ener;               /* Energy values                       */
  59   real **omega;             /* Omega values (energy dependent)     */
  60   real ***prob,***q;        /* Probability and Q                   */
  61 } t_pq_inel;
  62
  63 typedef struct {
  64   int  nener,n2Ddata;       /* Number of entries in the table      */
  65   real *ener;               /* Energy values                       */
  66   real **prob,**data;       /* Probability and data (energy loss)  */
  67 } t_p2Ddata;
  68
  69 static int realcomp(const void *a,const void *b)
  70 {
  71   real *ra,*rb;
  72
  73   ra = (real *)a;
  74   rb = (real *)b;
  75   if (ra < rb)
  76     return -1;
  77   else if (ra > rb)
  78     return 1;
  79   else
  80     return 0;
  81 }
  82
  83 static t_p2Ddata *read_p2Ddata(char *fn)
  84 {
  85   FILE    *fp;
  86   t_p2Ddata *p2Ddata;
  87   int     i,j;
  88   double  e,p,o;
  89
  90   fprintf(stdout,"Going to read %s\n",fn);
  91   fp = ffopen(fn,"r");
  92
  93   /* Allocate memory and set constants */
  94   snew(p2Ddata,1);
  95   if (fscanf(fp,"%d%d",&p2Ddata->nener,&p2Ddata->n2Ddata) != 2)
  96     gmx_fatal(FARGS,"I need two integers: nener, n in file %s",fn);
  97
  98   snew(p2Ddata->ener,p2Ddata->nener);
  99   snew(p2Ddata->prob,p2Ddata->nener);
 100   snew(p2Ddata->data,p2Ddata->nener);
 101
 102   /* Double loop to read data */
 103   for(i=0; (i<p2Ddata->nener); i++) {
 104     fprintf(stderr,"\rEnergy %d/%d",i+1,p2Ddata->nener);
 105     snew(p2Ddata->prob[i],p2Ddata->n2Ddata);
 106     snew(p2Ddata->data[i],p2Ddata->n2Ddata);
 107
 108     for(j=0; (j<p2Ddata->n2Ddata); j++) {
 109       fscanf(fp,"%lf%lf%lf",&e,&p,&o);
 110
 111       /* Consistency check */
 112       if (j==0)
 113         p2Ddata->ener[i] = e;
 114       else if (fabs(p2Ddata->ener[i]-e) > 1e-6*e)
 115         gmx_fatal(FARGS,"Inconsistent energy %f i=%d j=%d",e,i,j);
 116       p2Ddata->prob[i][j] = p;
 117       p2Ddata->data[i][j] = o;
 118     }
 119     /* There is some noise on the data, take it away by sorting,
 120      * because otherwise binary search does not work.
 121      * This is equivalent to shifting in the data slightly along the X-axis
 122      * but better than linear search with the "real" data.
 123      */
 124     qsort(p2Ddata->data[i],p2Ddata->n2Ddata,sizeof(p2Ddata->data[0][0]),
 125           realcomp);
 126   }
 127   fprintf(stderr,"\n");
 128
 129   ffclose(fp);
 130
 131   return p2Ddata;
 132 }
 133
 134 static t_pq_inel *read_pq(char *fn)
 135 {
 136   FILE      *fp;
 137   t_pq_inel *pq;
 138   int       i,j,k;
 139   double    e,p,o,t;
 140
 141   fprintf(stdout,"Going to read %s\n",fn);
 142   fp = ffopen(fn,"r");
 143
 144   /* Allocate memory and set constants */
 145   snew(pq,1);
 146   if (fscanf(fp,"%d%d%d",&pq->nener,&pq->nomega,&pq->nq) != 3)
 147     gmx_fatal(FARGS,"I need three integers: nener, nomega, nq in file %s",fn);
 148
 149   snew(pq->ener,pq->nener);
 150   snew(pq->omega,pq->nener);
 151   snew(pq->prob,pq->nener);
 152   snew(pq->q,pq->nener);
 153
 154   /* Triple loop to read data */
 155   for(i=0; (i<pq->nener); i++) {
 156     fprintf(stderr,"\rEnergy %d/%d",i+1,pq->nener);
 157     snew(pq->prob[i],pq->nomega);
 158     snew(pq->q[i],pq->nomega);
 159     snew(pq->omega[i],pq->nomega);
 160
 161     for(j=0; (j<pq->nomega); j++) {
 162       snew(pq->prob[i][j],pq->nq);
 163       snew(pq->q[i][j],pq->nq);
 164
 165       for(k=0; (k<pq->nq); k++) {
 166         fscanf(fp,"%lf%lf%lf%lf",&e,&o,&p,&t);
 167
 168         /* Consistency check */
 169         if ((j == 0) && (k == 0))
 170           pq->ener[i] = e;
 171         else if (fabs(pq->ener[i]-e) > 1e-6*e)
 172           gmx_fatal(FARGS,"Inconsistent energy %f i=%d j=%d k=%d",e,i,j,k);
 173
 174         if (k == 0)
 175           pq->omega[i][j] = o;
 176         else if (fabs(pq->omega[i][j]-o) > 1e-6*o)
 177           gmx_fatal(FARGS,"Inconsistent omega %f i=%d j=%d k=%d",o,i,j,k);
 178
 179         pq->prob[i][j][k] = p;
 180         pq->q[i][j][k] = t;
 181       }
 182     }
 183   }
 184   fprintf(stderr,"\n");
 185
 186   ffclose(fp);
 187
 188   return pq;
 189 }
 190
 191 static int my_bsearch(real val,int ndata,real data[],gmx_bool bLower)
 192 {
 193   int ilo,ihi,imed;
 194
 195   if (val < data[0])
 196     return -1;
 197   ilo = 0;
 198   ihi = ndata-1;
 199   while ((ihi - ilo) > 1) {
 200     imed = (ilo+ihi)/2;
 201     if (data[imed] > val)
 202       ihi = imed;
 203     else
 204       ilo = imed;
 205   }
 206   /* Now val should be in between data[ilo] and data[ihi] */
 207   /* Decide which one is closest */
 208   if (bLower || ((val-data[ilo]) <= (data[ihi]-val)))
 209     return ilo;
 210   else
 211     return ihi;
 212 }
 213
 214 static real interpolate2D(int nx,int ny,real *dx,real **data,
 215                           real x0,real fy,int nx0,int ny0)
 216 {
 217   real fx;
 218
 219   fx  = (x0-dx[nx0])/(dx[nx0+1]-dx[nx0]);
 220
 221   return (fx*fy*data[nx0][ny0] + fx*(1-fy)*data[nx0][ny0+1] +
 222           (1-fx)*fy*data[nx0+1][ny0] + (1-fx)*(1-fy)*data[nx0+1][ny0+1]);
 223 }
 224
 225 real get_omega(real ekin,int *seed,FILE *fp,char *fn)
 226 {
 227   static t_p2Ddata *p2Ddata = NULL;
 228   real r,ome,fx,fy;
 229   int  eindex,oindex;
 230
 231   if (p2Ddata == NULL)
 232     p2Ddata = read_p2Ddata(fn);
 233
 234   /* Get energy index by binary search */
 235   if ((eindex = my_bsearch(ekin,p2Ddata->nener,p2Ddata->ener,TRUE)) >= 0) {
 236 #ifdef DEBUG
 237     if (eindex >= p2Ddata->nener)
 238       gmx_fatal(FARGS,"eindex (%d) out of range (max %d) in get_omega",
 239                   eindex,p2Ddata->nener);
 240 #endif
 241
 242     /* Start with random number */
 243     r = rando(seed);
 244
 245     /* Do binary search in the energy table */
 246     if ((oindex = my_bsearch(r,p2Ddata->n2Ddata,p2Ddata->prob[eindex],TRUE)) >= 0) {
 247 #ifdef DEBUG
 248       if (oindex >= p2Ddata->n2Ddata)
 249         gmx_fatal(FARGS,"oindex (%d) out of range (max %d) in get_omega",
 250                     oindex,p2Ddata->n2Ddata);
 251 #endif
 252
 253       fy = ((r-p2Ddata->prob[eindex][oindex])/
 254             (p2Ddata->prob[eindex][oindex+1]-p2Ddata->prob[eindex][oindex]));
 255       ome = interpolate2D(p2Ddata->nener,p2Ddata->n2Ddata,p2Ddata->ener,
 256                           p2Ddata->data,ekin,fy,
 257                           eindex,oindex);
 258       /* ome = p2Ddata->data[eindex][oindex];*/
 259
 260       if (fp)
 261         fprintf(fp,"%8.3f  %8.3f\n",ome,r);
 262
 263       return ome;
 264     }
 265   }
 266   return 0;
 267 }
 268
 269 real get_theta_el(real ekin,int *seed,FILE *fp,char *fn)
 270 {
 271   static t_p2Ddata *p2Ddata = NULL;
 272   real r,theta;
 273   int  eindex,tindex;
 274
 275   if (p2Ddata == NULL)
 276     p2Ddata = read_p2Ddata(fn);
 277
 278   /* Get energy index by binary search */
 279   if ((eindex = my_bsearch(ekin,p2Ddata->nener,p2Ddata->ener,TRUE)) >= 0) {
 280
 281     /* Start with random number */
 282     r = rando(seed);
 283
 284     /* Do binary search in the energy table */
 285     if ((tindex = my_bsearch(r,p2Ddata->n2Ddata,p2Ddata->prob[eindex],FALSE)) >= 0) {
 286
 287       theta = p2Ddata->data[eindex][tindex];
 288
 289       if (fp)
 290         fprintf(fp,"%8.3f  %8.3f\n",theta,r);
 291
 292       return theta;
 293     }
 294   }
 295   return 0;
 296 }
 297
 298 real get_q_inel(real ekin,real omega,int *seed,FILE *fp,char *fn)
 299 {
 300   static t_pq_inel *pq = NULL;
 301   int    eindex,oindex,tindex;
 302   real   r,theta;
 303
 304   if (pq == NULL)
 305     pq = read_pq(fn);
 306
 307   /* Get energy index by binary search */
 308   if ((eindex = my_bsearch(ekin,pq->nener,pq->ener,TRUE)) >= 0) {
 309 #ifdef DEBUG
 310     if (eindex >= pq->nener)
 311       gmx_fatal(FARGS,"eindex out of range (%d >= %d)",eindex,pq->nener);
 312 #endif
 313
 314     /* Do binary search in the energy table */
 315     if ((oindex = my_bsearch(omega,pq->nomega,pq->omega[eindex],FALSE)) >= 0) {
 316 #ifdef DEBUG
 317       if (oindex >= pq->nomega)
 318         gmx_fatal(FARGS,"oindex out of range (%d >= %d)",oindex,pq->nomega);
 319 #endif
 320
 321       /* Start with random number */
 322       r = rando(seed);
 323
 324       if ((tindex = my_bsearch(r,pq->nq,pq->prob[eindex][oindex],FALSE)) >= 0) {
 325 #ifdef DEBUG
 326         if (tindex >= pq->nq)
 327           gmx_fatal(FARGS,"tindex out of range (%d >= %d)",tindex,pq->nq);
 328 #endif
 329
 330         theta = pq->q[eindex][oindex][tindex];
 331
 332         if (fp)
 333           fprintf(fp,"get_q_inel: %8.3f  %8.3f\n",theta,r);
 334
 335         return theta;
 336       }
 337     }
 338   }
 339   return 0;
 340 }
 341
 342 static int read_cross(char *fn,real **ener,real **cross,real factor)
 343 {
 344   char   **data=NULL;
 345   double e,c;
 346   int    i,j,n;
 347
 348   fprintf(stdout,"Going to read %s\n",fn);
 349   n = get_file(fn,&data);
 350
 351   /* Allocate memory */
 352   snew(*cross,n);
 353   snew(*ener,n);
 354   for(i=j=0; (i<n); i++) {
 355     if (sscanf(data[i],"%lf%lf",&e,&c) == 2) {
 356       (*ener)[j] = e;
 357       (*cross)[j] = c*factor;
 358       j++;
 359     }
 360     sfree(data[i]);
 361   }
 362   sfree(data);
 363   if (j != n)
 364     fprintf(stderr,"There were %d empty lines in file %s\n",n-j,fn);
 365
 366   return j;
 367 }
 368
 369 real cross_inel(real ekin,real rho,char *fn)
 370 {
 371   static real *ener  = NULL;
 372   static real *cross = NULL;
 373   static int  ninel;
 374   int eindex;
 375
 376   /* Read data at first call, convert A^2 to nm^2 */
 377   if (cross == NULL)
 378     ninel = read_cross(fn,&ener,&cross,rho*0.01);
 379
 380   /* Compute index with binary search */
 381   if ((eindex = my_bsearch(ekin,ninel,ener,FALSE)) >= 0) {
 382 #ifdef DEBUG
 383     if (eindex >= ninel)
 384       gmx_fatal(FARGS,"ekin = %f, ener[0] = %f, ener[%d] = %f",ekin,
 385                   ener[0],ener[ninel-1]);
 386 #endif
 387     return cross[eindex];
 388   }
 389
 390   return 0;
 391 }
 392
 393 real cross_el(real ekin,real rho,char *fn)
 394 {
 395   static real *ener  = NULL;
 396   static real *cross = NULL;
 397   static int  nel;
 398   int eindex;
 399
 400   /* Read data at first call, convert A^2 to nm^2  */
 401   if (cross == NULL)
 402     nel = read_cross(fn,&ener,&cross,rho*0.01);
 403
 404   /* Compute index with binary search */
 405   if ((eindex = my_bsearch(ekin,nel,ener,FALSE)) >= 0) {
 406 #ifdef DEBUG
 407     if (eindex >= nel)
 408       gmx_fatal(FARGS,"ekin = %f, ener[0] = %f, ener[%d] = %f",ekin,
 409                   ener[0],ener[nel-1]);
 410 #endif
 411
 412     return cross[eindex];
 413   }
 414   return 0;
 415 }
 416
 417 real band_ener(int *seed,FILE *fp,char *fn)
 418 {
 419   static real *ener  = NULL;
 420   static real *prob  = NULL;
 421   static int  nener;
 422   int  eindex;
 423   real r;
 424
 425   /* Read data at first call, misuse read_cross function */
 426   if (prob == NULL)
 427     nener = read_cross(fn,&ener,&prob,1.0);
 428
 429   r = rando(seed);
 430
 431   if ((eindex = my_bsearch(r,nener,prob,FALSE)) >= 0) {
 432 #ifdef DEBUG
 433     if (eindex >= nener)
 434       gmx_fatal(FARGS,"r = %f, prob[0] = %f, prob[%d] = %f",r,
 435                   prob[0],prob[nener-1]);
 436 #endif
 437
 438     if (fp)
 439       fprintf(fp,"%8.3f  %8.3f\n",ener[eindex],r);
 440
 441     return ener[eindex];
 442   }
 443   return 0;
 444 }
 445
 446 static void test_omega(FILE *fp,int *seed)
 447 {
 448   int i;
 449
 450   fprintf(fp,"Testing the energy loss tables\n");
 451   for(i=0; (i<1000); i++) {
 452     (void) get_omega(500*rando(seed),seed,fp,NULL);
 453   }
 454 }
 455
 456 static void test_q_inel(FILE *fp,int *seed)
 457 {
 458   int i;
 459
 460   fprintf(fp,"Testing the energy/omega dependent inelastic scattering q tables\n");
 461   for(i=0; (i<1000); i++) {
 462     (void) get_q_inel(500*rando(seed),400*rando(seed),seed,fp,NULL);
 463   }
 464 }
 465
 466 static void test_theta_el(FILE *fp,int *seed)
 467 {
 468   int i;
 469
 470   fprintf(fp,"Testing the energy dependent elastic scattering theta tables\n");
 471   for(i=0; (i<1000); i++) {
 472     (void) get_theta_el(500*rando(seed),seed,fp,NULL);
 473   }
 474 }
 475
 476 static void test_sigma_el(FILE *fp,real rho)
 477 {
 478   int  i;
 479
 480   fprintf(fp,"Testing the elastic cross sections table\n");
 481   for(i=0; (i<500); i++) {
 482     fprintf(fp,"%3d  %8.3f\n",i,cross_el(i,rho,NULL));
 483   }
 484 }
 485
 486 static void test_sigma_inel(FILE *fp,real rho)
 487 {
 488   int  i;
 489
 490   fprintf(fp,"Testing the inelastic cross sections table\n");
 491   for(i=0; (i<500); i++) {
 492     fprintf(fp,"%3d  %8.3f\n",i,cross_inel(i,rho,NULL));
 493   }
 494 }
 495
 496 static void test_band_ener(int *seed,FILE *fp)
 497 {
 498   int i;
 499
 500   for(i=0; (i<500); i++) {
 501     band_ener(seed,fp,NULL);
 502   }
 503 }
 504
 505 void init_tables(int nfile,t_filenm fnm[],real rho)
 506 {
 507   int  seed  = 1993;
 508   real ekin  = 20;
 509   real omega = 10;
 510
 511   (void) band_ener(&seed,NULL,opt2fn("-band",nfile,fnm));
 512   (void) cross_el(ekin,rho,opt2fn("-sigel",nfile,fnm));
 513   (void) cross_inel(ekin,rho,opt2fn("-sigin",nfile,fnm));
 514   (void) get_theta_el(ekin,&seed,NULL,opt2fn("-thetael",nfile,fnm));
 515   (void) get_omega(ekin,&seed,NULL,opt2fn("-eloss",nfile,fnm));
 516   (void) get_q_inel(ekin,omega,&seed,NULL,opt2fn("-qtrans",nfile,fnm));
 517 }
 518
 519 void test_tables(int *seed,char *fn,real rho)
 520 {
 521   FILE *fp;
 522
 523   fp = fopen(fn,"w");
 524
 525   test_omega(fp,seed);
 526   test_q_inel(fp,seed);
 527   test_theta_el(fp,seed);
 528   test_sigma_el(fp,rho);
 529   test_sigma_inel(fp,rho);
 530   test_band_ener(seed,fp);
 531
 532   fclose(fp);
 533 }
 534