?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Redefine the default boolean type to gmx_bool.
[gromacs.git] / src / contrib / ehdata.c
blob1eb6f1a67050e9efbbaae5032d7285a4c9c2adc1
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[],gmx_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
The ultimate molecular dynamics simulation package