?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' of git://git.gromacs.org/gromacs
[gromacs/adressmacs.git] / src / kernel / ionize.c
blob25dc3e5898e38bc2e2340765bb127e39cd6ab8c7
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.2.0
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-2004, 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 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <string.h>
40 #include "smalloc.h"
41 #include "typedefs.h"
42 #include "macros.h"
43 #include "random.h"
44 #include "gmx_random.h"
45 #include "physics.h"
46 #include "xvgr.h"
47 #include "vec.h"
48 #include "pbc.h"
49 #include "txtdump.h"
50 #include "ionize.h"
51 #include "names.h"
52 #include "futil.h"
53 #include "network.h"
54 #include "mtop_util.h"
55 #include "gmxfio.h"
56
57 typedef struct {
58 real photo,coh,incoh,incoh_abs;
59 } t_cross;
60
61 /* THIS TABLE HAS ADDED A 12 keV COLUMN TO HYDROGEN, CARBON, */
62 /* OXYGEN, NITROGEN AND SULPHUR BY FITTING A QUADRATIC TO THE */
63 /* POINTS 8keV, 10keV and 12keV - now contains 6, 8, 10, 12, */
64 /* 15 and 20 keV */
65 /* Units are barn. They are converted to nm^2 by multiplying */
66 /* by 1e-10, which is done in Imax (ionize.c) */
67 /* Update: contains energy 2 KeV and B, Na, Li, Al, Mg */
68 /* Update: data taken from NIST XCOM */
69
70 static const t_cross cross_sec_h[] = {
71 { 5.21E-01, 3.39E-01, 3.21E-01, -1 },
72 { 2.63e-2, 1.01e-1, 5.49e-1, 7.12e-3 },
73 { 9.79e-3, 6.18e-2, 5.83e-1, 9.60e-3 },
74 { 4.55e-3, 4.16e-2, 5.99e-1, 1.19e-2 },
75 { 1.52e-3, 2.79e-2, 6.08e-1, 1.41e-2 },
76 { 1.12e-3, 1.96e-2, 6.09e-1, 1.73e-2 },
77 { 4.16e-4, 1.13e-2, 6.07e-1, 2.23e-2 }
78 };
79 static const t_cross cross_sec_c[] = {
80 { 3.10E+3, 1.42E+1, 1.03E+0, -1 },
81 { 1.99e+2, 5.88e+0, 2.29e+0, 3.06e-2 },
82 { 8.01e+1, 4.22e+0, 2.56e+0, 4.38e-2 },
83 { 3.92e+1, 3.26e+0, 2.74e+0, 5.72e-2 },
84 { 2.19e+1, 2.55e+0, 2.88e+0, 7.07e-2 },
85 { 1.06e+1, 1.97e+0, 3.04e+0, 9.15e-2 },
86 { 4.15e+0, 1.30e+0, 3.20e+0, 1.24e-1 }
87 };
88 static const t_cross cross_sec_n[] = {
89 { 5.78E+3, 2.13E+1, 1.11E+0, -1 },
90 { 3.91e+2, 8.99e+0, 2.49e+0, 3.43e-2 },
91 { 1.59e+2, 6.29e+0, 2.86e+0, 5.01e-2 },
92 { 7.88e+1, 4.76e+0, 3.10e+0, 6.57e-2 },
93 { 4.42e+1, 3.66e+0, 3.28e+0, 8.13e-2 },
94 { 2.16e+1, 2.82e+0, 3.46e+0, 1.05e-1 },
95 { 8.52e+0, 1.88e+0, 3.65e+0, 1.43e-1 }
96 };
97 static const t_cross cross_sec_o[] = {
98 { 9.74E+3, 3.00E+1, 1.06E+0, -1 },
99 { 6.90e+2, 1.33e+1, 2.66e+0, 3.75e-2 },
100 { 2.84e+2, 9.21e+0, 3.14e+0, 5.62e-2 },
101 { 1.42e+2, 6.85e+0, 3.44e+0, 7.43e-2 },
102 { 8.01e+1, 5.18e+0, 3.66e+0, 9.20e-2 },
103 { 3.95e+1, 3.97e+0, 3.87e+0, 1.18e-1 },
104 { 1.57e+1, 2.64e+0, 4.10e+0, 1.61e-1 }
105 };
106 static const t_cross cross_sec_s[] = {
107 { 1.07E+5, 1.15E+2, 2.03E+0, -1 },
108 { 1.10e+4, 5.54e+1, 3.98e+0, 5.42e-2 },
109 { 4.91e+3, 4.29e+1, 4.71e+0, 8.38e-2 },
110 { 2.58e+3, 3.36e+1, 5.32e+0, 1.16e-1 },
111 { 1.52e+3, 2.64e+1, 5.81e+0, 1.48e-1 },
112 { 7.82e+2, 1.97e+1, 6.36e+0, 2.00e-1 },
113 { 3.29e+2, 1.29e+1, 6.94e+0, 2.80e-1 }
114 };
115 static const t_cross cross_sec_mg[] = {
116 { 7.79E+4, 7.19E+1, 1.34E+0, -1 },
117 { 3.75E+3, 3.75E+1, 3.18E+0, -1 },
118 { 1.61E+3, 2.75E+1, 3.91E+0, -1 },
119 { 8.25E+2, 2.06E+1, 4.47E+0, -1 },
120 { 4.75E+2, 1.61E+1, 4.88E+0, -1 },
121 { 2.40E+2, 1.16E+1, 5.32E+0, -1 },
122 { 9.82E+1, 7.59E+0, 5.74E+0, -1 }
123 };
124 static const t_cross cross_sec_al[] = {
125 { 1.01E+5, 8.24E+1, 1.51E+0, -1 },
126 { 5.12E+3, 4.32E+1, 3.45E+0, -1 },
127 { 2.22E+3, 3.24E+1, 4.16E+0, -1 },
128 { 1.14E+3, 2.47E+1, 4.74E+0, -1 },
129 { 6.63E+2, 1.93E+1, 5.19E+0, -1 },
130 { 3.37E+2, 1.41E+1, 5.67E+0, -1 },
131 { 1.39E+2, 9.17E+0, 6.14E+0, -1 }
132 };
133 static const t_cross cross_sec_b[] = {
134 { 2.86E+3, 1.05E+1, 8.20E-1, -1 },
135 { 9.38E+1, 3.76E+0, 1.92E+0, -1 },
136 { 3.72E+1, 2.81E+0, 2.15E+0, -1 },
137 { 1.80E+1, 2.20E+0, 2.31E+0, -1 },
138 { 9.92E+0, 1.77E+0, 2.44E+0, -1 },
139 { 4.77E+0, 1.32E+0, 2.58E+0, -1 },
140 { 1.84E+0, 8.56E-1, 2.71E+0, -1 }
141 };
142 static const t_cross cross_sec_na[] = {
143 { 5.80E+4, 6.27E+1, 1.01E+0, -1 },
144 { 2.65E+3, 3.17E+1, 2.95E+0, -1 },
145 { 1.13E+3, 2.26E+1, 3.67E+0, -1 },
146 { 5.74E+2, 1.68E+1, 4.20E+0, -1 },
147 { 3.28E+2, 1.30E+1, 4.58E+0, -1 },
148 { 1.65E+2, 9.43E+0, 4.96E+0, -1 },
149 { 6.71E+1, 6.16E+0, 5.34E+0, -1 }
150 };
151 static const t_cross cross_sec_li[] = {
152 { 3.08E+2, 3.37E+0, 6.38E-1, -1 },
153 { 8.60E+0, 1.60E+0, 1.18E+0, -1 },
154 { 3.31E+0, 1.16E+0, 1.36E+0, -1 },
155 { 1.57E+0, 8.63E-1, 1.48E+0, -1 },
156 { 8.50E-1, 6.59E-1, 1.57E+0, -1 },
157 { 4.01E-1, 4.63E-1, 1.64E+0, -1 },
158 { 1.52E-1, 2.85E-1, 1.70E+0, -1 }
159 };
160
161 typedef struct {
162 const char *name;
163 int nel;
164 const t_cross *cross;
165 } t_element;
166
167 static const t_element element[] = {
168 { "H", 1, cross_sec_h },
169 { "C", 6, cross_sec_c },
170 { "N", 7, cross_sec_n },
171 { "O", 8, cross_sec_o },
172 { "S", 16, cross_sec_s },
173 { "LI", 3, cross_sec_li },
174 { "B", 5, cross_sec_b },
175 { "NA", 11, cross_sec_na },
176 { "MG", 12, cross_sec_mg },
177 { "AL", 13, cross_sec_al },
178 { "AR", 20, cross_sec_s }, /* This is not correct! */
179 { "EL", 1, cross_sec_h } /* This is not correct! */
180 };
181 #define NELEM asize(element)
182
183 /*
184 * In the first column the binding energy of the K-electrons;
185 * THIS IS IN eV, which matches the photon energies.
186 * In the second column the binding energy of the outer shell electrons
187 * The third column describes the photoelectric cross sections,
188 * where this now gives the fraction of photoelectric events
189 * which correspond to K-shell events, I called f_j in my
190 * notes:
191 * The final column (a new column) now gives the values for the lifetimes
192 * in ps.
193 */
194 typedef struct {
195 real E_K,E_L,Prob_K,tau;
196 } t_recoil;
197
198 const t_recoil recoil[] = {
199 { 0.0, 0.0, 0.0, 0},
200 { 0.0136, 0.0, 0.0, 0},
201 { 0.0246, 0.0, 0.0, 0},
202 { 0.055, 0.005, 0.960, 0.012},
203 { 0.117, 0.009, 0.956, 0.012},
204 { 0.192, 0.008, 0.952, 0.012},
205 { 0.284, 0.011, 0.950, 0.0113},
206 { 0.402, 0.015, 0.950, 0.0083},
207 { 0.532, 0.014, 0.936, 0.0066},
208 { 0.687, 0.017, 0.928, 0.0045},
209 { 0.874, 0.031, 0.922, 0.0033},
210 { 1.072, 0.041, 0.933, 0.0028},
211 { 1.305, 0.054, 0.927, 0.0022},
212 { 1.560, 0.077, 0.922, 0.0019},
213 { 1.839, 0.105, 0.918, 0.00165},
214 { 2.146, 0.133, 0.921, 0.00145},
215 { 2.472, 0.166, 0.908, 0.00130},
216 { 2.822, 0.212, 0.902, 0.0012},
217 { 3.203, 0.247, 0.902, 0.0010},
218 { 3.607, 0.298, 0.894, 0.00095},
219 { 4.038, 0.348, 0.890, 0.00085},
220 { 4.490, 0.404, 0.886, 0.00078},
221 { 4.966, 0.458, 0.882, 0.00073},
222 { 5.465, 0.516, 0.885, 0.00062},
223 { 5.989, 0.578, 0.883, 0.00055},
224 { 6.539, 0.645, 0.880, 0.00049},
225 { 7.112, 0.713, 0.877, 0.00044}
226 };
227
228 #define PREFIX "IONIZE: "
229
230 enum { eionCYL, eionSURF, eionGAUSS, eionNR };
231
232 enum { ecollPHOTO, ecollINELASTIC, ecollNR };
233
234 typedef struct {
235 int z,n,k;
236 real fj,sigPh,sigIn,vAuger;
237 } t_cross_atom;
238
239 /* BEGIN GLOBAL VARIABLES */
240
241 /*
242 UGLY HACK
243 The 2 in this list doesn't really mean 2, but 2.5 keV as
244 it's checked inside the code and added 0.5 when needed.
245 */
246
247 static int Energies[] = { 2, 6, 8, 10, 12, 15, 20 };
248 static int ionize_seed = 1993;
249 #define NENER asize(Energies)
250 /* END GLOBAL VARIABLES */
251
252 void dump_ca(FILE *fp,t_cross_atom *ca,int i,const char *file,int line)
253 {
254 fprintf(fp,PREFIX"(line %d) atom %d, z = %d, n = %d, k = %d\n",
255 line,i,ca->z,ca->n,ca->k);
256 }
257
258 t_cross_atom *mk_cross_atom(FILE *log,t_mdatoms *md,
259 gmx_mtop_t *mtop,int Eindex)
260 {
261 int elem_index[] = { 0, 0, 0, 5, 0, 6, 1, 2, 3, 0, 0, 7, 8, 9, 0, 0, 4, 0, 0, 0, 10, 11 };
262 t_cross_atom *ca;
263 int *elemcnt;
264 char *cc,*resname;
265 int i,j,resnr;
266
267 fprintf(log,PREFIX"Filling data structure for ionization\n");
268 fprintf(log,PREFIX"Warning: all fj values set to 0.95 for now\n");
269 snew(ca,md->nr);
270 snew(elemcnt,NELEM+1);
271 for(i=0; (i<md->nr); i++) {
272 ca[i].n = 0;
273 ca[i].k = 0;
274 /* This code does not work for domain decomposition */
275 gmx_mtop_atominfo_global(mtop,i,&cc,&resnr,&resname);
276 for(j=0; (j<NELEM); j++)
277 if (strncmp(cc,element[j].name,strlen(element[j].name)) == 0) {
278 ca[i].z = element[j].nel;
279 break;
280 }
281 if (j == NELEM)
282 gmx_fatal(FARGS,PREFIX"Don't know number of electrons for %s",cc);
283
284 elemcnt[j]++;
285
286 ca[i].sigPh = element[elem_index[ca[i].z]].cross[Eindex].photo;
287 ca[i].sigIn = element[elem_index[ca[i].z]].cross[Eindex].incoh;
288 ca[i].fj = recoil[ca[i].z].Prob_K;
289 switch (ca[i].z) {
290 case 6:
291 ca[i].vAuger = 0.904;
292 break;
293 case 7:
294 ca[i].vAuger = 0.920;
295 break;
296 case 8:
297 ca[i].vAuger = 0.929;
298 break;
299 case 3: /* probably not correct! */
300 ca[i].vAuger = 0.9;
301 break;
302 case 5: /* probably not correct! */
303 ca[i].vAuger = 0.9;
304 break;
305 case 11: /* probably not correct! */
306 case 12: /* probably not correct! */
307 case 13: /* probably not correct! */
308 case 16:
309 case 20:
310 ca[i].vAuger = 1.0;
311 break;
312 default:
313 ca[i].vAuger= -1;
314 }
315 }
316
317 fprintf(log,PREFIX"You have the following elements in your system (%d atoms):\n"PREFIX,md->nr);
318 for(j=0; (j<NELEM); j++)
319 if (elemcnt[j] > 0)
320 fprintf(log," %s: %d",element[j].name,elemcnt[j]);
321 fprintf(log," atoms\n");
322
323 sfree(elemcnt);
324
325 return ca;
326 }
327
328 int number_K(t_cross_atom *ca)
329 {
330 if (ca->z <= 2)
331 return ca->z-ca->n;
332 else
333 return 2-ca->k;
334 }
335
336 int number_L(t_cross_atom *ca)
337 {
338 return ca->k-2+ca->z-ca->n;
339 }
340
341 real xray_cross_section(int eColl,t_cross_atom *ca)
342 {
343 real c=0;
344 int nK,nL;
345
346 switch (eColl) {
347 case ecollPHOTO:
348 nK = number_K(ca);
349 nL = number_L(ca);
350 if (ca->z == 1)
351 c = ca->sigPh;
352 else if (ca->z == 2)
353 c = ca->sigPh*0.5;
354 else
355 c = (nK*0.5*ca->fj + nL/(ca->z-2)*(1-ca->fj))*ca->sigPh;
356 break;
357 case ecollINELASTIC:
358 c = (ca->z-ca->n)*ca->sigIn/ca->z;
359 break;
360 default:
361 gmx_fatal(FARGS,"No such collision type %d\n",eColl);
362 }
363 return c;
364 }
365
366 real prob_K(int eColl,t_cross_atom *ca)
367 {
368 real Pl,Pk,P=0;
369
370 if ((ca->z <= 2) || (ca->z == ca->n))
371 return 0;
372
373 switch (eColl) {
374 case ecollPHOTO:
375 Pl = (ca->k-2+ca->z-ca->n)*(1-ca->fj)/(ca->z-2);
376 Pk = (2-ca->k)*ca->fj*0.5;
377 P = Pk/(Pl+Pk);
378 break;
379 case ecollINELASTIC:
380 P = (2-ca->k)/(ca->z-ca->n);
381 break;
382 default:
383 gmx_fatal(FARGS,"No such collision type %d\n",eColl);
384 }
385 return P;
386 }
387
388 double myexp(double x)
389 {
390 if (x < -70)
391 return 0.0;
392 else
393 return exp(x);
394 }
395
396 real ptheta_incoh(int Eindex,real theta)
397 /* theta should be in degrees */
398 {
399 /* These numbers generated by fitting 5 gaussians to the real function
400 * that describes the probability for theta.
401 * We use symmetry in the gaussian (see 180-angle) therefore there
402 * are fewer parameters (only 8 per energylevel).
403 */
404 static double ppp[NENER][8] = {
405 { -0.00295169, 10.4847, 0.0341099, /*-*/43.1963,
406 -0.0164054, 30.2452, 71.0311, 2.50282 },
407 { -0.00370852, 9.02037, 0.100559, /*-*/42.9962,
408 -0.0537891, 35.5077, 71.4305, 1.05515 },
409 { -0.00427039, 7.86831, 0.118042, /*-*/45.9846,
410 -0.0634505, 38.6134, 70.3857, 0.240082 },
411 { -0.004514, 7.0728, 0.13464, /*-*/48.213,
412 -0.0723, 41.06, 69.38, -0.02 },
413 { -0.00488796, 5.87988, 0.159574, /*-*/51.5556,
414 -0.0855767, 44.7307, 69.0251, -0.414604 },
415 { -0.00504604, 4.56299, 0.201064, /*-*/54.8599,
416 -0.107153, 48.7016, 68.8212, -0.487699 }
417 };
418 double g1,g2,g3,g4,g5,ptheta;
419
420 g1 = myexp(-0.5*sqr((theta-ppp[Eindex][7])/ppp[Eindex][1]));
421 g2 = myexp(-0.5*sqr((theta-180+ppp[Eindex][7])/ppp[Eindex][1]));
422 g3 = myexp(-0.5*sqr((theta-90)/ppp[Eindex][3]));
423 g4 = myexp(-0.5*sqr((theta-ppp[Eindex][6])/ppp[Eindex][5]));
424 g5 = myexp(-0.5*sqr((theta-180+ppp[Eindex][6])/ppp[Eindex][5]));
425
426 ptheta = ppp[Eindex][0]*(g1+g2) + ppp[Eindex][2]*g3 + ppp[Eindex][4]*(g4+g5);
427
428 return ptheta;
429 }
430
431 real rand_theta_incoh(int Eindex,int *seed)
432 {
433 #define NINTP 450
434 #define prev (1-cur)
435 static gmx_bool bFirst = TRUE;
436 static real **intp;
437 static int i,j,cur=1;
438 real rrr,dx;
439 real y[2];
440
441 dx = 90.0/(real)NINTP;
442 if (bFirst) {
443 /* Compute cumulative integrals of all probability distributions */
444 snew(intp,NENER);
445 for(i=0; (i<NENER); i++) {
446 snew(intp[i],NINTP+1);
447 y[prev] = ptheta_incoh(i,0.0);
448 /*sum = y[prev];*/
449 for(j=1; (j<=NINTP); j++) {
450 y[cur] = ptheta_incoh(i,j*dx);
451 /*sum += y[cur];*/
452 intp[i][j] = intp[i][j-1] + (y[cur]+y[prev])*dx;
453 cur = prev;
454 }
455 }
456 if (debug) {
457 fprintf(debug,"Integrated probability functions for theta incoherent\n");
458 for(j=0; (j<NINTP); j++) {
459 fprintf(debug,"%10f",dx*j);
460 for(i=0; (i<NENER); i++)
461 fprintf(debug," %10f",intp[i][j]);
462 fprintf(debug,"\n");
463 }
464 }
465 bFirst = FALSE;
466 }
467
468 rrr = rando(seed);
469 for(j=0; (j<NINTP) && (rrr > intp[Eindex][j]); j++)
470 ;
471
472 return (j-1+(rrr-intp[Eindex][j-1])/(intp[Eindex][j]-intp[Eindex][j-1]))*dx;
473 }
474
475 static void polar2cart(real phi,real theta,rvec v)
476 {
477 v[XX] = cos(phi)*sin(theta);
478 v[YY] = sin(phi)*sin(theta);
479 v[ZZ] = cos(theta);
480 }
481
482 void rand_vector(rvec v,int *seed)
483 {
484 real theta,phi;
485
486 theta = 180.0*rando(seed);
487 phi = 360.0*rando(seed);
488 polar2cart(phi,theta,v);
489 }
490
491 real electron_cross_section(FILE *fp,rvec v,real mass,int nelec)
492 {
493 /* Compute cross section for electrons */
494 real T,B,U,S,Q,R,N,t,u,lnt,sigma;
495 real a0 = 0.05292; /* nm */
496
497 /* Have to determine T (kinetic energy of electron) */
498 T = 0.5*mass*iprod(v,v);
499
500 /* R is the binding energy of the electron in hydrogen */
501 R = 13.61*ELECTRONVOLT;
502
503 /* Have to determine the binding energy B, differs per orbital of course */
504 B = R;
505
506 /* Have to determine the orbital kinetic energy U */
507 U = R;
508
509 /* Have to know number of electrons */
510 N = nelec;
511
512 /* Magic constant Q */
513 Q = 1;
514
515 /* Some help variables */
516 t = T/B;
517 u = U/B;
518 S = 4*M_PI*sqr(a0)*N*sqr(R/B);
519 lnt = log(t);
520
521 /* Resulting variable */
522 sigma = (S/(t+u+1))*( 0.5*Q*lnt*(1-1/sqr(t)) + (2-Q)*(1-1/t-lnt/(t+1)) );
523
524 return sigma;
525 }
526
527 gmx_bool khole_decay(FILE *fp,t_cross_atom *ca,rvec x[],rvec v[],int ion,
528 int *seed,real dt)
529 {
530 rvec dv;
531 real factor;
532
533 if ((ca->vAuger < 0) || (recoil[ca->z].tau == 0)) {
534 dump_ca(stderr,ca,ion,__FILE__,__LINE__);
535 exit(1);
536 }
537 if ((rando(seed) < dt/recoil[ca->z].tau) && (number_L(ca)>1)) {
538 if (debug)
539 fprintf(debug,"DECAY: Going to decay a k hole\n");
540 ca->n++;
541 ca->k--;
542 /* Generate random vector */
543 rand_vector(dv,seed);
544
545 factor = ca->vAuger;
546 if (debug)
547 fprintf(debug,"DECAY: factor=%10g, dv = (%8.3f, %8.3f, %8.3f)\n",
548 factor,dv[XX],dv[YY],dv[ZZ]);
549 svmul(factor,dv,dv);
550 rvec_inc(v[ion],dv);
551
552 return TRUE;
553 }
554 else
555 return FALSE;
556 }
557
558 void ionize(FILE *fp,const output_env_t oenv,t_mdatoms *md,gmx_mtop_t *mtop,
559 real t,t_inputrec *ir, rvec x[],rvec v[],int start,int end,
560 matrix box,t_commrec *cr)
561 {
562 static FILE *xvg,*ion;
563 static const char *const_leg[] = { "Probability", "Primary Ionization", "Integral over PI", "KHole-Decay", "Integral over KD" };
564 static gmx_bool bFirst = TRUE;
565 static real t0,imax,width,rho,nphot;
566 static real interval;
567 static int dq_tot,nkd_tot,mode,ephot;
568 static t_cross_atom *ca;
569 static int Eindex=-1;
570 static gmx_rng_t gaussrand=NULL;
571
572 real factor,E_lost=0;
573 real pt,ptot,pphot,pcoll[ecollNR],tmax;
574 real hboxx,hboxy,rho2;
575 rvec dv,ddv;
576 gmx_bool bIonize=FALSE,bKHole,bL,bDOIT;
577 int i,k,kk,m,nK,nL,dq,nkh,nkdecay;
578 int *nionize,*nkhole,*ndecay,nbuf[2];
579 char **leg;
580
581
582 if (bFirst) {
583 /* Get parameters for gaussian photon pulse from inputrec */
584 t0 = ir->userreal1; /* Peak of the gaussian pulse */
585 nphot = ir->userreal2; /* Intensity */
586 width = ir->userreal3; /* Width of the peak (in time) */
587 rho = ir->userreal4; /* Diameter of the focal spot (nm) */
588 ionize_seed= ir->userint1; /* Random seed for stochastic ionization */
589 ephot = ir->userint2; /* Energy of the photons */
590 mode = ir->userint3; /* Mode of ionizing */
591 interval = 0.001*ir->userint4; /* Interval between pulses (ps) */
592 gaussrand=gmx_rng_init(ionize_seed);
593
594 snew(leg,asize(const_leg));
595 for(i=0;i<asize(const_leg);i++)
596 leg[i]=strdup(const_leg[i]);
597
598 if ((width <= 0) || (nphot <= 0))
599 gmx_fatal(FARGS,"Your parameters for ionization are not set properly\n"
600 "width (userreal3) = %f, nphot (userreal2) = %f",
601 width,nphot);
602
603 if ((mode < 0) || (mode >= eionNR))
604 gmx_fatal(FARGS,"Ionization mode (userint3)"
605 " should be in the range 0 .. %d",eionNR-1);
606
607 switch (mode) {
608 case eionCYL:
609 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
610 break;
611 case eionSURF:
612 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
613 break;
614 }
615 if (ionize_seed == 0)
616 ionize_seed = make_seed();
617 if (PAR(cr)) {
618 for(i=0; (i<cr->nodeid); i++)
619 ionize_seed = INT_MAX*rando(&ionize_seed);
620 fprintf(fp,PREFIX"Modifying seed on parallel processor to %d\n",
621 ionize_seed);
622 }
623
624 for(Eindex=0; (Eindex < NENER) && (Energies[Eindex] != ephot); Eindex++)
625 ;
626 if (Eindex == NENER)
627 gmx_fatal(FARGS,PREFIX"Energy level of %d keV not supported",ephot);
628
629 /* Initiate cross section data etc. */
630 ca = mk_cross_atom(fp,md,mtop,Eindex);
631
632 dq_tot = 0;
633 nkd_tot = 0;
634
635 xvg = xvgropen("ionize.xvg","Ionization Events","Time (ps)","()",oenv);
636 xvgr_legend(xvg,asize(leg),(const char**)leg,oenv);
637 ion = gmx_fio_fopen("ionize.log","w");
638
639 fprintf(fp,PREFIX"Parameters for ionization events:\n");
640 fprintf(fp,PREFIX"Imax = %g, t0 = %g, width = %g, seed = %d\n"
641 PREFIX"# Photons = %g, rho = %g, ephot = %d (keV)\n",
642 imax,t0,width,ionize_seed,nphot,rho,ephot);
643 fprintf(fp,PREFIX"Electron_mass: %10.3e(keV) Atomic_mass: %10.3e(keV)\n"
644 PREFIX"Speed_of_light: %10.3e(nm/ps)\n",
645 ELECTRONMASS_keV,ATOMICMASS_keV,SPEED_OF_LIGHT);
646 fprintf(fp,PREFIX"Interval between shots: %g ps\n",interval);
647 fprintf(fp,PREFIX"Eindex = %d\n",Eindex);
648 fprintf(fp,PREFIX"Doing ionizations for atoms %d - %d\n",start,end);
649
650 fflush(fp);
651
652 bFirst = FALSE;
653 }
654
655 /******************************************************
656 *
657 * H E R E S T A R T S I O N I Z A T I O N
658 *
659 ******************************************************/
660
661 /* Calculate probability */
662 tmax = t0;
663 if (interval > 0)
664 while (t > (tmax+interval*0.5))
665 tmax += interval;
666 /* End when t <= t0 + (N+0.5) interval */
667
668 pt = imax*ir->delta_t*exp(-0.5*sqr((t-tmax)/width));
669 dq = 0;
670 nkdecay = 0;
671
672 hboxx = 0.5*box[XX][XX];
673 hboxy = 0.5*box[YY][YY];
674 rho2 = sqr(rho);
675
676 /* Arrays for ionization statistics */
677 snew(nionize,md->nr);
678 snew(nkhole,md->nr);
679 snew(ndecay,md->nr);
680
681 /* Loop over atoms */
682 for(i=start; (i<end); i++) {
683 /* Loop over collision types */
684 bKHole = FALSE;
685 bIonize = FALSE;
686 for(k=0; (k<ecollNR); k++)
687 /* Determine cross section for this collision type */
688 pcoll[k]= pt*xray_cross_section(k,&(ca[i]));
689
690 /* Total probability of ionisation */
691 ptot = 1 - (1-pcoll[ecollPHOTO])*(1-pcoll[ecollINELASTIC]);
692 if (debug && (i==0))
693 fprintf(debug,PREFIX"Ptot = %g, t = %g\n",ptot,t);
694
695 /* Check whether to ionize this guy */
696 bDOIT = FALSE;
697 switch (mode) {
698 case eionCYL:
699 bDOIT = (((rando(&ionize_seed) < ptot) && (ca[i].n < ca[i].z)) &&
700 ((sqr(x[i][XX] - hboxx) + sqr(x[i][YY] - hboxy)) < rho2));
701 break;
702 case eionSURF:
703 bDOIT = FALSE;
704 break;
705 default:
706 gmx_fatal(FARGS,"Unknown ionization mode %d (%s, line %d)",mode,
707 __FILE__,__LINE__);
708 }
709
710 if (bDOIT) {
711 clear_rvec(dv);
712
713 /* The relative probability for a photoellastic event is given by: */
714 pphot = pcoll[ecollPHOTO]/(pcoll[ecollPHOTO]+pcoll[ecollINELASTIC]);
715
716 if (rando(&ionize_seed) < pphot)
717 k = ecollPHOTO;
718 else
719 k = ecollINELASTIC;
720
721 /* If a random number is smaller than the probability for
722 * an L ionization than do that. Note that the probability
723 * may be zero (H, He), but the < instead of <= covers that.
724 */
725 nK = number_K(&ca[i]);
726 nL = number_L(&ca[i]);
727 bL = (nK == 0) || ( (nL > 0) && (rando(&ionize_seed) > prob_K(k,&(ca[i]))));
728
729 switch (k) {
730 case ecollPHOTO: {
731 /* Select which one to take by yet another random numer */
732 real theta,phi;
733
734 /* Get parameters for photoelestic effect */
735 /* Note that in the article this is called 2 theta */
736 theta = DEG2RAD*gmx_rng_gaussian_table(gaussrand)*26.0+70.0;
737
738 phi = 2*M_PI*rando(&ionize_seed);
739
740 if (bL)
741 E_lost = ephot-recoil[ca[i].z].E_L*(ca[i].n+1);
742 if (ephot == 2) E_lost += 0.5; /* Real energy should be 2.5 KeV*/
743 else {
744 E_lost = ephot-recoil[ca[i].z].E_K;
745 if (ephot == 2) E_lost += 0.5; /* Real energy should be 2.5 KeV*/
746 if ((ca[i].z > 2) && (nL > 0))
747 bKHole = TRUE;
748 }
749 if (debug)
750 fprintf(debug,"i = %d, nK = %d, nL = %d, bL = %s, bKHole = %s\n",
751 i,nK,nL,BOOL(bL),BOOL(bKHole));
752 if (E_lost < 0) {
753 E_lost = 0.0;
754 bIonize = FALSE;
755 bKHole = FALSE;
756 }
757 else {
758 /* Compute the components of the velocity vector */
759 factor = ((ELECTRONMASS_keV/(ATOMICMASS_keV*md->massT[i]))*
760 (SPEED_OF_LIGHT*sqrt(2*E_lost/ELECTRONMASS_keV)));
761
762 /* Subtract momentum of recoiling electron */
763 polar2cart(phi,theta,ddv);
764 for(m=0; (m<DIM); m++)
765 dv[m] -= factor*ddv[m];
766
767 if (debug)
768 pr_rvec(debug,0,"ELL",dv,DIM,TRUE);
769
770 bIonize = TRUE;
771 }
772 break;
773 }
774 case ecollINELASTIC: {
775 real theta,Ebind,Eelec;
776
777 if (bL)
778 Ebind = (ca[i].n+1)*recoil[ca[i].z].E_L;
779 else {
780 Ebind = recoil[ca[i].z].E_K;
781 if ((ca[i].z > 2) && (nL > 0))
782 bKHole = TRUE;
783 }
784 theta = DEG2RAD*rand_theta_incoh(Eindex,&ionize_seed);
785 Eelec = (sqr(ephot)/512)*(1-cos(2*theta));
786 if (ephot == 2) Eelec = (sqr(ephot+.5)/512)*(1-cos(2*theta)); /* Real energy should be 2.5 KeV*/
787 bIonize = (Ebind <= Eelec);
788 bKHole = bKHole && bIonize;
789 if (debug)
790 fprintf(debug,PREFIX"Ebind: %g, Eelectron: %g\n",Ebind,Eelec);
791 if (!bIonize) {
792 /* Subtract momentum of recoiling photon */
793 /*phi = 2*M_PI*rando(&ionize_seed);
794 bKHole = FALSE;
795 factor = ephot*438;
796 dv[XX] -= factor*cos(phi)*sin(theta);
797 dv[YY] -= factor*sin(phi)*sin(theta);
798 dv[ZZ] -= factor*cos(theta);
799 */
800 if (debug)
801 pr_rvec(debug,0,"INELL",dv,DIM,TRUE);
802 }
803 break;
804 }
805 default:
806 gmx_fatal(FARGS,"Ga direct naar de gevangenis. Ga niet langs start");
807 }
808 if (bIonize) {
809 /* First increase the charge */
810 if (ca[i].n < ca[i].z) {
811 md->chargeA[i] += 1.0;
812 md->chargeB[i] += 1.0;
813 ca[i].n++;
814 dq ++;
815 }
816 if (debug) {
817 fprintf(debug,"Random-dv[%3d] = %10.3e,%10.3e,%10.3e,"
818 " ephot = %d, Elost=%10.3e\n",
819 i,dv[XX],dv[YY],dv[ZZ],ephot,E_lost);
820 }
821 }
822 /* Now actually add the impulse to the velocities */
823 for(m=0; (m<DIM); m++)
824 v[i][m] += dv[m];
825 if (bKHole) {
826 ca[i].k ++;
827 nkhole[i]++;
828 }
829 else if (bIonize)
830 nionize[i]++;
831 }
832
833 /* Now check old event: Loop over k holes! */
834 nkh = ca[i].k;
835 for (kk = 0; (kk < nkh); kk++)
836 if (khole_decay(fp,&(ca[i]),x,v,i,&ionize_seed,ir->delta_t)) {
837 nkdecay ++;
838 ndecay[i]++;
839 md->chargeA[i] += 1.0;
840 md->chargeB[i] += 1.0;
841 }
842
843 if (debug && (ca[i].n > 0))
844 dump_ca(debug,&(ca[i]),i,__FILE__,__LINE__);
845 }
846
847 /* Sum events for statistics if necessary */
848 if (PAR(cr)) {
849 gmx_sumi(md->nr,nionize,cr);
850 gmx_sumi(md->nr,nkhole,cr);
851 gmx_sumi(md->nr,ndecay,cr);
852 nbuf[0] = dq; nbuf[1] = nkdecay;
853 gmx_sumi(2,nbuf,cr);
854 dq = nbuf[0]; nkdecay = nbuf[1];
855 }
856 /* Now sum global events on this timestep to cumulative numbers */
857 dq_tot += dq;
858 nkd_tot += nkdecay;
859
860 /* Printing time */
861 if (MASTER(cr)) {
862 /* Print data to the file that holds ionization events per atom */
863 fprintf(ion,"%12.8f",t);
864 for(i=0; (i<md->nr); i++) {
865 if (nionize[i])
866 fprintf(ion," I:%d",i+1);
867 if (nkhole[i])
868 fprintf(ion," K:%d",i+1);
869 if (ndecay[i])
870 fprintf(ion," D:%d",i+1);
871 }
872 fprintf(ion,"\n");
873 if (debug)
874 fflush(ion);
875
876 /* Print statictics to file */
877 fprintf(xvg,"%10.5f %10.3e %6d %6d %6d %6d",
878 t,pt,dq,dq_tot,nkdecay,nkd_tot);
879 fprintf(xvg,"\n");
880 if (debug)
881 fflush(xvg);
882 }
883 sfree(nionize);
884 sfree(nkhole);
885 sfree(ndecay);
886 }
887
got merge a long time ago