2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, 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.
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.
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.
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.
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.
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.
47 #include "gmx_random.h"
57 #include "mtop_util.h"
61 real photo,coh,incoh,incoh_abs;
64 /* THIS TABLE HAS ADDED A 12 keV COLUMN TO HYDROGEN, CARBON, */
65 /* OXYGEN, NITROGEN AND SULPHUR BY FITTING A QUADRATIC TO THE */
66 /* POINTS 8keV, 10keV and 12keV - now contains 6, 8, 10, 12, */
68 /* Units are barn. They are converted to nm^2 by multiplying */
69 /* by 1e-10, which is done in Imax (ionize.c) */
70 /* Update: contains energy 2 KeV and B, Na, Li, Al, Mg */
71 /* Update: data taken from NIST XCOM */
73 static const t_cross cross_sec_h[] = {
74 { 5.21E-01, 3.39E-01, 3.21E-01, -1 },
75 { 2.63e-2, 1.01e-1, 5.49e-1, 7.12e-3 },
76 { 9.79e-3, 6.18e-2, 5.83e-1, 9.60e-3 },
77 { 4.55e-3, 4.16e-2, 5.99e-1, 1.19e-2 },
78 { 1.52e-3, 2.79e-2, 6.08e-1, 1.41e-2 },
79 { 1.12e-3, 1.96e-2, 6.09e-1, 1.73e-2 },
80 { 4.16e-4, 1.13e-2, 6.07e-1, 2.23e-2 }
82 static const t_cross cross_sec_c[] = {
83 { 3.10E+3, 1.42E+1, 1.03E+0, -1 },
84 { 1.99e+2, 5.88e+0, 2.29e+0, 3.06e-2 },
85 { 8.01e+1, 4.22e+0, 2.56e+0, 4.38e-2 },
86 { 3.92e+1, 3.26e+0, 2.74e+0, 5.72e-2 },
87 { 2.19e+1, 2.55e+0, 2.88e+0, 7.07e-2 },
88 { 1.06e+1, 1.97e+0, 3.04e+0, 9.15e-2 },
89 { 4.15e+0, 1.30e+0, 3.20e+0, 1.24e-1 }
91 static const t_cross cross_sec_n[] = {
92 { 5.78E+3, 2.13E+1, 1.11E+0, -1 },
93 { 3.91e+2, 8.99e+0, 2.49e+0, 3.43e-2 },
94 { 1.59e+2, 6.29e+0, 2.86e+0, 5.01e-2 },
95 { 7.88e+1, 4.76e+0, 3.10e+0, 6.57e-2 },
96 { 4.42e+1, 3.66e+0, 3.28e+0, 8.13e-2 },
97 { 2.16e+1, 2.82e+0, 3.46e+0, 1.05e-1 },
98 { 8.52e+0, 1.88e+0, 3.65e+0, 1.43e-1 }
100 static const t_cross cross_sec_o[] = {
101 { 9.74E+3, 3.00E+1, 1.06E+0, -1 },
102 { 6.90e+2, 1.33e+1, 2.66e+0, 3.75e-2 },
103 { 2.84e+2, 9.21e+0, 3.14e+0, 5.62e-2 },
104 { 1.42e+2, 6.85e+0, 3.44e+0, 7.43e-2 },
105 { 8.01e+1, 5.18e+0, 3.66e+0, 9.20e-2 },
106 { 3.95e+1, 3.97e+0, 3.87e+0, 1.18e-1 },
107 { 1.57e+1, 2.64e+0, 4.10e+0, 1.61e-1 }
109 static const t_cross cross_sec_s[] = {
110 { 1.07E+5, 1.15E+2, 2.03E+0, -1 },
111 { 1.10e+4, 5.54e+1, 3.98e+0, 5.42e-2 },
112 { 4.91e+3, 4.29e+1, 4.71e+0, 8.38e-2 },
113 { 2.58e+3, 3.36e+1, 5.32e+0, 1.16e-1 },
114 { 1.52e+3, 2.64e+1, 5.81e+0, 1.48e-1 },
115 { 7.82e+2, 1.97e+1, 6.36e+0, 2.00e-1 },
116 { 3.29e+2, 1.29e+1, 6.94e+0, 2.80e-1 }
118 static const t_cross cross_sec_mg[] = {
119 { 7.79E+4, 7.19E+1, 1.34E+0, -1 },
120 { 3.75E+3, 3.75E+1, 3.18E+0, -1 },
121 { 1.61E+3, 2.75E+1, 3.91E+0, -1 },
122 { 8.25E+2, 2.06E+1, 4.47E+0, -1 },
123 { 4.75E+2, 1.61E+1, 4.88E+0, -1 },
124 { 2.40E+2, 1.16E+1, 5.32E+0, -1 },
125 { 9.82E+1, 7.59E+0, 5.74E+0, -1 }
127 static const t_cross cross_sec_al[] = {
128 { 1.01E+5, 8.24E+1, 1.51E+0, -1 },
129 { 5.12E+3, 4.32E+1, 3.45E+0, -1 },
130 { 2.22E+3, 3.24E+1, 4.16E+0, -1 },
131 { 1.14E+3, 2.47E+1, 4.74E+0, -1 },
132 { 6.63E+2, 1.93E+1, 5.19E+0, -1 },
133 { 3.37E+2, 1.41E+1, 5.67E+0, -1 },
134 { 1.39E+2, 9.17E+0, 6.14E+0, -1 }
136 static const t_cross cross_sec_b[] = {
137 { 2.86E+3, 1.05E+1, 8.20E-1, -1 },
138 { 9.38E+1, 3.76E+0, 1.92E+0, -1 },
139 { 3.72E+1, 2.81E+0, 2.15E+0, -1 },
140 { 1.80E+1, 2.20E+0, 2.31E+0, -1 },
141 { 9.92E+0, 1.77E+0, 2.44E+0, -1 },
142 { 4.77E+0, 1.32E+0, 2.58E+0, -1 },
143 { 1.84E+0, 8.56E-1, 2.71E+0, -1 }
145 static const t_cross cross_sec_na[] = {
146 { 5.80E+4, 6.27E+1, 1.01E+0, -1 },
147 { 2.65E+3, 3.17E+1, 2.95E+0, -1 },
148 { 1.13E+3, 2.26E+1, 3.67E+0, -1 },
149 { 5.74E+2, 1.68E+1, 4.20E+0, -1 },
150 { 3.28E+2, 1.30E+1, 4.58E+0, -1 },
151 { 1.65E+2, 9.43E+0, 4.96E+0, -1 },
152 { 6.71E+1, 6.16E+0, 5.34E+0, -1 }
154 static const t_cross cross_sec_li[] = {
155 { 3.08E+2, 3.37E+0, 6.38E-1, -1 },
156 { 8.60E+0, 1.60E+0, 1.18E+0, -1 },
157 { 3.31E+0, 1.16E+0, 1.36E+0, -1 },
158 { 1.57E+0, 8.63E-1, 1.48E+0, -1 },
159 { 8.50E-1, 6.59E-1, 1.57E+0, -1 },
160 { 4.01E-1, 4.63E-1, 1.64E+0, -1 },
161 { 1.52E-1, 2.85E-1, 1.70E+0, -1 }
167 const t_cross *cross;
170 static const t_element element[] = {
171 { "H", 1, cross_sec_h },
172 { "C", 6, cross_sec_c },
173 { "N", 7, cross_sec_n },
174 { "O", 8, cross_sec_o },
175 { "S", 16, cross_sec_s },
176 { "LI", 3, cross_sec_li },
177 { "B", 5, cross_sec_b },
178 { "NA", 11, cross_sec_na },
179 { "MG", 12, cross_sec_mg },
180 { "AL", 13, cross_sec_al },
181 { "AR", 20, cross_sec_s }, /* This is not correct! */
182 { "EL", 1, cross_sec_h } /* This is not correct! */
184 #define NELEM asize(element)
187 * In the first column the binding energy of the K-electrons;
188 * THIS IS IN eV, which matches the photon energies.
189 * In the second column the binding energy of the outer shell electrons
190 * The third column describes the photoelectric cross sections,
191 * where this now gives the fraction of photoelectric events
192 * which correspond to K-shell events, I called f_j in my
194 * The final column (a new column) now gives the values for the lifetimes
198 real E_K,E_L,Prob_K,tau;
201 const t_recoil recoil[] = {
203 { 0.0136, 0.0, 0.0, 0},
204 { 0.0246, 0.0, 0.0, 0},
205 { 0.055, 0.005, 0.960, 0.012},
206 { 0.117, 0.009, 0.956, 0.012},
207 { 0.192, 0.008, 0.952, 0.012},
208 { 0.284, 0.011, 0.950, 0.0113},
209 { 0.402, 0.015, 0.950, 0.0083},
210 { 0.532, 0.014, 0.936, 0.0066},
211 { 0.687, 0.017, 0.928, 0.0045},
212 { 0.874, 0.031, 0.922, 0.0033},
213 { 1.072, 0.041, 0.933, 0.0028},
214 { 1.305, 0.054, 0.927, 0.0022},
215 { 1.560, 0.077, 0.922, 0.0019},
216 { 1.839, 0.105, 0.918, 0.00165},
217 { 2.146, 0.133, 0.921, 0.00145},
218 { 2.472, 0.166, 0.908, 0.00130},
219 { 2.822, 0.212, 0.902, 0.0012},
220 { 3.203, 0.247, 0.902, 0.0010},
221 { 3.607, 0.298, 0.894, 0.00095},
222 { 4.038, 0.348, 0.890, 0.00085},
223 { 4.490, 0.404, 0.886, 0.00078},
224 { 4.966, 0.458, 0.882, 0.00073},
225 { 5.465, 0.516, 0.885, 0.00062},
226 { 5.989, 0.578, 0.883, 0.00055},
227 { 6.539, 0.645, 0.880, 0.00049},
228 { 7.112, 0.713, 0.877, 0.00044}
231 #define PREFIX "IONIZE: "
233 enum { eionCYL, eionSURF, eionGAUSS, eionNR };
235 enum { ecollPHOTO, ecollINELASTIC, ecollNR };
239 real fj,sigPh,sigIn,vAuger;
242 /* BEGIN GLOBAL VARIABLES */
246 The 2 in this list doesn't really mean 2, but 2.5 keV as
247 it's checked inside the code and added 0.5 when needed.
250 static int Energies[] = { 2, 6, 8, 10, 12, 15, 20 };
251 static int ionize_seed = 1993;
252 #define NENER asize(Energies)
253 /* END GLOBAL VARIABLES */
255 void dump_ca(FILE *fp,t_cross_atom *ca,int i,const char *file,int line)
257 fprintf(fp,PREFIX"(line %d) atom %d, z = %d, n = %d, k = %d\n",
258 line,i,ca->z,ca->n,ca->k);
261 t_cross_atom *mk_cross_atom(FILE *log,t_mdatoms *md,
262 gmx_mtop_t *mtop,int Eindex)
264 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 };
270 fprintf(log,PREFIX"Filling data structure for ionization\n");
271 fprintf(log,PREFIX"Warning: all fj values set to 0.95 for now\n");
273 snew(elemcnt,NELEM+1);
274 for(i=0; (i<md->nr); i++) {
277 /* This code does not work for domain decomposition */
278 gmx_mtop_atominfo_global(mtop,i,&cc,&resnr,&resname);
279 for(j=0; (j<NELEM); j++)
280 if (strncmp(cc,element[j].name,strlen(element[j].name)) == 0) {
281 ca[i].z = element[j].nel;
285 gmx_fatal(FARGS,PREFIX"Don't know number of electrons for %s",cc);
289 ca[i].sigPh = element[elem_index[ca[i].z]].cross[Eindex].photo;
290 ca[i].sigIn = element[elem_index[ca[i].z]].cross[Eindex].incoh;
291 ca[i].fj = recoil[ca[i].z].Prob_K;
294 ca[i].vAuger = 0.904;
297 ca[i].vAuger = 0.920;
300 ca[i].vAuger = 0.929;
302 case 3: /* probably not correct! */
305 case 5: /* probably not correct! */
308 case 11: /* probably not correct! */
309 case 12: /* probably not correct! */
310 case 13: /* probably not correct! */
320 fprintf(log,PREFIX"You have the following elements in your system (%d atoms):\n"PREFIX,md->nr);
321 for(j=0; (j<NELEM); j++)
323 fprintf(log," %s: %d",element[j].name,elemcnt[j]);
324 fprintf(log," atoms\n");
331 int number_K(t_cross_atom *ca)
339 int number_L(t_cross_atom *ca)
341 return ca->k-2+ca->z-ca->n;
344 real xray_cross_section(int eColl,t_cross_atom *ca)
358 c = (nK*0.5*ca->fj + nL/(ca->z-2)*(1-ca->fj))*ca->sigPh;
361 c = (ca->z-ca->n)*ca->sigIn/ca->z;
364 gmx_fatal(FARGS,"No such collision type %d\n",eColl);
369 real prob_K(int eColl,t_cross_atom *ca)
373 if ((ca->z <= 2) || (ca->z == ca->n))
378 Pl = (ca->k-2+ca->z-ca->n)*(1-ca->fj)/(ca->z-2);
379 Pk = (2-ca->k)*ca->fj*0.5;
383 P = (2-ca->k)/(ca->z-ca->n);
386 gmx_fatal(FARGS,"No such collision type %d\n",eColl);
391 double myexp(double x)
399 real ptheta_incoh(int Eindex,real theta)
400 /* theta should be in degrees */
402 /* These numbers generated by fitting 5 gaussians to the real function
403 * that describes the probability for theta.
404 * We use symmetry in the gaussian (see 180-angle) therefore there
405 * are fewer parameters (only 8 per energylevel).
407 static double ppp[NENER][8] = {
408 { -0.00295169, 10.4847, 0.0341099, /*-*/43.1963,
409 -0.0164054, 30.2452, 71.0311, 2.50282 },
410 { -0.00370852, 9.02037, 0.100559, /*-*/42.9962,
411 -0.0537891, 35.5077, 71.4305, 1.05515 },
412 { -0.00427039, 7.86831, 0.118042, /*-*/45.9846,
413 -0.0634505, 38.6134, 70.3857, 0.240082 },
414 { -0.004514, 7.0728, 0.13464, /*-*/48.213,
415 -0.0723, 41.06, 69.38, -0.02 },
416 { -0.00488796, 5.87988, 0.159574, /*-*/51.5556,
417 -0.0855767, 44.7307, 69.0251, -0.414604 },
418 { -0.00504604, 4.56299, 0.201064, /*-*/54.8599,
419 -0.107153, 48.7016, 68.8212, -0.487699 }
421 double g1,g2,g3,g4,g5,ptheta;
423 g1 = myexp(-0.5*sqr((theta-ppp[Eindex][7])/ppp[Eindex][1]));
424 g2 = myexp(-0.5*sqr((theta-180+ppp[Eindex][7])/ppp[Eindex][1]));
425 g3 = myexp(-0.5*sqr((theta-90)/ppp[Eindex][3]));
426 g4 = myexp(-0.5*sqr((theta-ppp[Eindex][6])/ppp[Eindex][5]));
427 g5 = myexp(-0.5*sqr((theta-180+ppp[Eindex][6])/ppp[Eindex][5]));
429 ptheta = ppp[Eindex][0]*(g1+g2) + ppp[Eindex][2]*g3 + ppp[Eindex][4]*(g4+g5);
434 real rand_theta_incoh(int Eindex,int *seed)
438 static gmx_bool bFirst = TRUE;
440 static int i,j,cur=1;
444 dx = 90.0/(real)NINTP;
446 /* Compute cumulative integrals of all probability distributions */
448 for(i=0; (i<NENER); i++) {
449 snew(intp[i],NINTP+1);
450 y[prev] = ptheta_incoh(i,0.0);
452 for(j=1; (j<=NINTP); j++) {
453 y[cur] = ptheta_incoh(i,j*dx);
455 intp[i][j] = intp[i][j-1] + (y[cur]+y[prev])*dx;
460 fprintf(debug,"Integrated probability functions for theta incoherent\n");
461 for(j=0; (j<NINTP); j++) {
462 fprintf(debug,"%10f",dx*j);
463 for(i=0; (i<NENER); i++)
464 fprintf(debug," %10f",intp[i][j]);
472 for(j=0; (j<NINTP) && (rrr > intp[Eindex][j]); j++)
475 return (j-1+(rrr-intp[Eindex][j-1])/(intp[Eindex][j]-intp[Eindex][j-1]))*dx;
478 static void polar2cart(real phi,real theta,rvec v)
480 v[XX] = cos(phi)*sin(theta);
481 v[YY] = sin(phi)*sin(theta);
485 void rand_vector(rvec v,int *seed)
489 theta = 180.0*rando(seed);
490 phi = 360.0*rando(seed);
491 polar2cart(phi,theta,v);
494 gmx_bool khole_decay(FILE *fp,t_cross_atom *ca,rvec x[],rvec v[],int ion,
500 if ((ca->vAuger < 0) || (recoil[ca->z].tau == 0)) {
501 dump_ca(stderr,ca,ion,__FILE__,__LINE__);
504 if ((rando(seed) < dt/recoil[ca->z].tau) && (number_L(ca)>1)) {
506 fprintf(debug,"DECAY: Going to decay a k hole\n");
509 /* Generate random vector */
510 rand_vector(dv,seed);
514 fprintf(debug,"DECAY: factor=%10g, dv = (%8.3f, %8.3f, %8.3f)\n",
515 factor,dv[XX],dv[YY],dv[ZZ]);
525 void ionize(FILE *fp,const output_env_t oenv,t_mdatoms *md,gmx_mtop_t *mtop,
526 real t,t_inputrec *ir, rvec x[],rvec v[],int start,int end,
527 matrix box,t_commrec *cr)
529 static FILE *xvg,*ion;
530 static const char *const_leg[] = { "Probability", "Primary Ionization", "Integral over PI", "KHole-Decay", "Integral over KD" };
531 static gmx_bool bFirst = TRUE;
532 static real t0,imax,width,rho,nphot;
533 static real interval;
534 static int dq_tot,nkd_tot,mode,ephot;
535 static t_cross_atom *ca;
536 static int Eindex=-1;
537 static gmx_rng_t gaussrand=NULL;
539 real factor,E_lost=0;
540 real pt,ptot,pphot,pcoll[ecollNR],tmax;
541 real hboxx,hboxy,rho2;
543 gmx_bool bIonize=FALSE,bKHole,bL,bDOIT;
544 int i,k,kk,m,nK,nL,dq,nkh,nkdecay;
545 int *nionize,*nkhole,*ndecay,nbuf[2];
550 /* Get parameters for gaussian photon pulse from inputrec */
551 t0 = ir->userreal1; /* Peak of the gaussian pulse */
552 nphot = ir->userreal2; /* Intensity */
553 width = ir->userreal3; /* Width of the peak (in time) */
554 rho = ir->userreal4; /* Diameter of the focal spot (nm) */
555 ionize_seed= ir->userint1; /* Random seed for stochastic ionization */
556 ephot = ir->userint2; /* Energy of the photons */
557 mode = ir->userint3; /* Mode of ionizing */
558 interval = 0.001*ir->userint4; /* Interval between pulses (ps) */
559 gaussrand=gmx_rng_init(ionize_seed);
561 snew(leg,asize(const_leg));
562 for(i=0;i<asize(const_leg);i++)
563 leg[i]=strdup(const_leg[i]);
565 if ((width <= 0) || (nphot <= 0))
566 gmx_fatal(FARGS,"Your parameters for ionization are not set properly\n"
567 "width (userreal3) = %f, nphot (userreal2) = %f",
570 if ((mode < 0) || (mode >= eionNR))
571 gmx_fatal(FARGS,"Ionization mode (userint3)"
572 " should be in the range 0 .. %d",eionNR-1);
576 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
579 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
582 if (ionize_seed == 0)
583 ionize_seed = make_seed();
585 for(i=0; (i<cr->nodeid); i++)
586 ionize_seed = INT_MAX*rando(&ionize_seed);
587 fprintf(fp,PREFIX"Modifying seed on parallel processor to %d\n",
591 for(Eindex=0; (Eindex < NENER) && (Energies[Eindex] != ephot); Eindex++)
594 gmx_fatal(FARGS,PREFIX"Energy level of %d keV not supported",ephot);
596 /* Initiate cross section data etc. */
597 ca = mk_cross_atom(fp,md,mtop,Eindex);
602 xvg = xvgropen("ionize.xvg","Ionization Events","Time (ps)","()",oenv);
603 xvgr_legend(xvg,asize(leg),(const char**)leg,oenv);
604 ion = gmx_fio_fopen("ionize.log","w");
606 fprintf(fp,PREFIX"Parameters for ionization events:\n");
607 fprintf(fp,PREFIX"Imax = %g, t0 = %g, width = %g, seed = %d\n"
608 PREFIX"# Photons = %g, rho = %g, ephot = %d (keV)\n",
609 imax,t0,width,ionize_seed,nphot,rho,ephot);
610 fprintf(fp,PREFIX"Electron_mass: %10.3e(keV) Atomic_mass: %10.3e(keV)\n"
611 PREFIX"Speed_of_light: %10.3e(nm/ps)\n",
612 ELECTRONMASS_keV,ATOMICMASS_keV,SPEED_OF_LIGHT);
613 fprintf(fp,PREFIX"Interval between shots: %g ps\n",interval);
614 fprintf(fp,PREFIX"Eindex = %d\n",Eindex);
615 fprintf(fp,PREFIX"Doing ionizations for atoms %d - %d\n",start,end);
622 /******************************************************
624 * H E R E S T A R T S I O N I Z A T I O N
626 ******************************************************/
628 /* Calculate probability */
631 while (t > (tmax+interval*0.5))
633 /* End when t <= t0 + (N+0.5) interval */
635 pt = imax*ir->delta_t*exp(-0.5*sqr((t-tmax)/width));
639 hboxx = 0.5*box[XX][XX];
640 hboxy = 0.5*box[YY][YY];
643 /* Arrays for ionization statistics */
644 snew(nionize,md->nr);
648 /* Loop over atoms */
649 for(i=start; (i<end); i++) {
650 /* Loop over collision types */
653 for(k=0; (k<ecollNR); k++)
654 /* Determine cross section for this collision type */
655 pcoll[k]= pt*xray_cross_section(k,&(ca[i]));
657 /* Total probability of ionisation */
658 ptot = 1 - (1-pcoll[ecollPHOTO])*(1-pcoll[ecollINELASTIC]);
660 fprintf(debug,PREFIX"Ptot = %g, t = %g\n",ptot,t);
662 /* Check whether to ionize this guy */
666 bDOIT = (((rando(&ionize_seed) < ptot) && (ca[i].n < ca[i].z)) &&
667 ((sqr(x[i][XX] - hboxx) + sqr(x[i][YY] - hboxy)) < rho2));
673 gmx_fatal(FARGS,"Unknown ionization mode %d (%s, line %d)",mode,
680 /* The relative probability for a photoellastic event is given by: */
681 pphot = pcoll[ecollPHOTO]/(pcoll[ecollPHOTO]+pcoll[ecollINELASTIC]);
683 if (rando(&ionize_seed) < pphot)
688 /* If a random number is smaller than the probability for
689 * an L ionization than do that. Note that the probability
690 * may be zero (H, He), but the < instead of <= covers that.
692 nK = number_K(&ca[i]);
693 nL = number_L(&ca[i]);
694 bL = (nK == 0) || ( (nL > 0) && (rando(&ionize_seed) > prob_K(k,&(ca[i]))));
698 /* Select which one to take by yet another random numer */
701 /* Get parameters for photoelestic effect */
702 /* Note that in the article this is called 2 theta */
703 theta = DEG2RAD*gmx_rng_gaussian_table(gaussrand)*26.0+70.0;
705 phi = 2*M_PI*rando(&ionize_seed);
708 E_lost = ephot-recoil[ca[i].z].E_L*(ca[i].n+1);
709 if (ephot == 2) E_lost += 0.5; /* Real energy should be 2.5 KeV*/
711 E_lost = ephot-recoil[ca[i].z].E_K;
712 if (ephot == 2) E_lost += 0.5; /* Real energy should be 2.5 KeV*/
713 if ((ca[i].z > 2) && (nL > 0))
717 fprintf(debug,"i = %d, nK = %d, nL = %d, bL = %s, bKHole = %s\n",
718 i,nK,nL,EBOOL(bL),EBOOL(bKHole));
725 /* Compute the components of the velocity vector */
726 factor = ((ELECTRONMASS_keV/(ATOMICMASS_keV*md->massT[i]))*
727 (SPEED_OF_LIGHT*sqrt(2*E_lost/ELECTRONMASS_keV)));
729 /* Subtract momentum of recoiling electron */
730 polar2cart(phi,theta,ddv);
731 for(m=0; (m<DIM); m++)
732 dv[m] -= factor*ddv[m];
735 pr_rvec(debug,0,"ELL",dv,DIM,TRUE);
741 case ecollINELASTIC: {
742 real theta,Ebind,Eelec;
745 Ebind = (ca[i].n+1)*recoil[ca[i].z].E_L;
747 Ebind = recoil[ca[i].z].E_K;
748 if ((ca[i].z > 2) && (nL > 0))
751 theta = DEG2RAD*rand_theta_incoh(Eindex,&ionize_seed);
752 Eelec = (sqr(ephot)/512)*(1-cos(2*theta));
753 if (ephot == 2) Eelec = (sqr(ephot+.5)/512)*(1-cos(2*theta)); /* Real energy should be 2.5 KeV*/
754 bIonize = (Ebind <= Eelec);
755 bKHole = bKHole && bIonize;
757 fprintf(debug,PREFIX"Ebind: %g, Eelectron: %g\n",Ebind,Eelec);
759 /* Subtract momentum of recoiling photon */
760 /*phi = 2*M_PI*rando(&ionize_seed);
763 dv[XX] -= factor*cos(phi)*sin(theta);
764 dv[YY] -= factor*sin(phi)*sin(theta);
765 dv[ZZ] -= factor*cos(theta);
768 pr_rvec(debug,0,"INELL",dv,DIM,TRUE);
773 gmx_fatal(FARGS,"Ga direct naar de gevangenis. Ga niet langs start");
776 /* First increase the charge */
777 if (ca[i].n < ca[i].z) {
778 md->chargeA[i] += 1.0;
779 md->chargeB[i] += 1.0;
784 fprintf(debug,"Random-dv[%3d] = %10.3e,%10.3e,%10.3e,"
785 " ephot = %d, Elost=%10.3e\n",
786 i,dv[XX],dv[YY],dv[ZZ],ephot,E_lost);
789 /* Now actually add the impulse to the velocities */
790 for(m=0; (m<DIM); m++)
800 /* Now check old event: Loop over k holes! */
802 for (kk = 0; (kk < nkh); kk++)
803 if (khole_decay(fp,&(ca[i]),x,v,i,&ionize_seed,ir->delta_t)) {
806 md->chargeA[i] += 1.0;
807 md->chargeB[i] += 1.0;
810 if (debug && (ca[i].n > 0))
811 dump_ca(debug,&(ca[i]),i,__FILE__,__LINE__);
814 /* Sum events for statistics if necessary */
816 gmx_sumi(md->nr,nionize,cr);
817 gmx_sumi(md->nr,nkhole,cr);
818 gmx_sumi(md->nr,ndecay,cr);
819 nbuf[0] = dq; nbuf[1] = nkdecay;
821 dq = nbuf[0]; nkdecay = nbuf[1];
823 /* Now sum global events on this timestep to cumulative numbers */
829 /* Print data to the file that holds ionization events per atom */
830 fprintf(ion,"%12.8f",t);
831 for(i=0; (i<md->nr); i++) {
833 fprintf(ion," I:%d",i+1);
835 fprintf(ion," K:%d",i+1);
837 fprintf(ion," D:%d",i+1);
843 /* Print statictics to file */
844 fprintf(xvg,"%10.5f %10.3e %6d %6d %6d %6d",
845 t,pt,dq,dq_tot,nkdecay,nkd_tot);