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