?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
removed group non-boneded call with verlet scheme
[gromacs/AngularHB.git] / src / kernel / ionize.c
blob41d44f1514214b2802538ce3b58ad608400adf13
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 gmx_bool khole_decay(FILE *fp,t_cross_atom *ca,rvec x[],rvec v[],int ion,
492 int *seed,real dt)
493 {
494 rvec dv;
495 real factor;
496
497 if ((ca->vAuger < 0) || (recoil[ca->z].tau == 0)) {
498 dump_ca(stderr,ca,ion,__FILE__,__LINE__);
499 exit(1);
500 }
501 if ((rando(seed) < dt/recoil[ca->z].tau) && (number_L(ca)>1)) {
502 if (debug)
503 fprintf(debug,"DECAY: Going to decay a k hole\n");
504 ca->n++;
505 ca->k--;
506 /* Generate random vector */
507 rand_vector(dv,seed);
508
509 factor = ca->vAuger;
510 if (debug)
511 fprintf(debug,"DECAY: factor=%10g, dv = (%8.3f, %8.3f, %8.3f)\n",
512 factor,dv[XX],dv[YY],dv[ZZ]);
513 svmul(factor,dv,dv);
514 rvec_inc(v[ion],dv);
515
516 return TRUE;
517 }
518 else
519 return FALSE;
520 }
521
522 void ionize(FILE *fp,const output_env_t oenv,t_mdatoms *md,gmx_mtop_t *mtop,
523 real t,t_inputrec *ir, rvec x[],rvec v[],int start,int end,
524 matrix box,t_commrec *cr)
525 {
526 static FILE *xvg,*ion;
527 static const char *const_leg[] = { "Probability", "Primary Ionization", "Integral over PI", "KHole-Decay", "Integral over KD" };
528 static gmx_bool bFirst = TRUE;
529 static real t0,imax,width,rho,nphot;
530 static real interval;
531 static int dq_tot,nkd_tot,mode,ephot;
532 static t_cross_atom *ca;
533 static int Eindex=-1;
534 static gmx_rng_t gaussrand=NULL;
535
536 real factor,E_lost=0;
537 real pt,ptot,pphot,pcoll[ecollNR],tmax;
538 real hboxx,hboxy,rho2;
539 rvec dv,ddv;
540 gmx_bool bIonize=FALSE,bKHole,bL,bDOIT;
541 int i,k,kk,m,nK,nL,dq,nkh,nkdecay;
542 int *nionize,*nkhole,*ndecay,nbuf[2];
543 char **leg;
544
545
546 if (bFirst) {
547 /* Get parameters for gaussian photon pulse from inputrec */
548 t0 = ir->userreal1; /* Peak of the gaussian pulse */
549 nphot = ir->userreal2; /* Intensity */
550 width = ir->userreal3; /* Width of the peak (in time) */
551 rho = ir->userreal4; /* Diameter of the focal spot (nm) */
552 ionize_seed= ir->userint1; /* Random seed for stochastic ionization */
553 ephot = ir->userint2; /* Energy of the photons */
554 mode = ir->userint3; /* Mode of ionizing */
555 interval = 0.001*ir->userint4; /* Interval between pulses (ps) */
556 gaussrand=gmx_rng_init(ionize_seed);
557
558 snew(leg,asize(const_leg));
559 for(i=0;i<asize(const_leg);i++)
560 leg[i]=strdup(const_leg[i]);
561
562 if ((width <= 0) || (nphot <= 0))
563 gmx_fatal(FARGS,"Your parameters for ionization are not set properly\n"
564 "width (userreal3) = %f, nphot (userreal2) = %f",
565 width,nphot);
566
567 if ((mode < 0) || (mode >= eionNR))
568 gmx_fatal(FARGS,"Ionization mode (userint3)"
569 " should be in the range 0 .. %d",eionNR-1);
570
571 switch (mode) {
572 case eionCYL:
573 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
574 break;
575 case eionSURF:
576 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
577 break;
578 }
579 if (ionize_seed == 0)
580 ionize_seed = make_seed();
581 if (PAR(cr)) {
582 for(i=0; (i<cr->nodeid); i++)
583 ionize_seed = INT_MAX*rando(&ionize_seed);
584 fprintf(fp,PREFIX"Modifying seed on parallel processor to %d\n",
585 ionize_seed);
586 }
587
588 for(Eindex=0; (Eindex < NENER) && (Energies[Eindex] != ephot); Eindex++)
589 ;
590 if (Eindex == NENER)
591 gmx_fatal(FARGS,PREFIX"Energy level of %d keV not supported",ephot);
592
593 /* Initiate cross section data etc. */
594 ca = mk_cross_atom(fp,md,mtop,Eindex);
595
596 dq_tot = 0;
597 nkd_tot = 0;
598
599 xvg = xvgropen("ionize.xvg","Ionization Events","Time (ps)","()",oenv);
600 xvgr_legend(xvg,asize(leg),(const char**)leg,oenv);
601 ion = gmx_fio_fopen("ionize.log","w");
602
603 fprintf(fp,PREFIX"Parameters for ionization events:\n");
604 fprintf(fp,PREFIX"Imax = %g, t0 = %g, width = %g, seed = %d\n"
605 PREFIX"# Photons = %g, rho = %g, ephot = %d (keV)\n",
606 imax,t0,width,ionize_seed,nphot,rho,ephot);
607 fprintf(fp,PREFIX"Electron_mass: %10.3e(keV) Atomic_mass: %10.3e(keV)\n"
608 PREFIX"Speed_of_light: %10.3e(nm/ps)\n",
609 ELECTRONMASS_keV,ATOMICMASS_keV,SPEED_OF_LIGHT);
610 fprintf(fp,PREFIX"Interval between shots: %g ps\n",interval);
611 fprintf(fp,PREFIX"Eindex = %d\n",Eindex);
612 fprintf(fp,PREFIX"Doing ionizations for atoms %d - %d\n",start,end);
613
614 fflush(fp);
615
616 bFirst = FALSE;
617 }
618
619 /******************************************************
620 *
621 * H E R E S T A R T S I O N I Z A T I O N
622 *
623 ******************************************************/
624
625 /* Calculate probability */
626 tmax = t0;
627 if (interval > 0)
628 while (t > (tmax+interval*0.5))
629 tmax += interval;
630 /* End when t <= t0 + (N+0.5) interval */
631
632 pt = imax*ir->delta_t*exp(-0.5*sqr((t-tmax)/width));
633 dq = 0;
634 nkdecay = 0;
635
636 hboxx = 0.5*box[XX][XX];
637 hboxy = 0.5*box[YY][YY];
638 rho2 = sqr(rho);
639
640 /* Arrays for ionization statistics */
641 snew(nionize,md->nr);
642 snew(nkhole,md->nr);
643 snew(ndecay,md->nr);
644
645 /* Loop over atoms */
646 for(i=start; (i<end); i++) {
647 /* Loop over collision types */
648 bKHole = FALSE;
649 bIonize = FALSE;
650 for(k=0; (k<ecollNR); k++)
651 /* Determine cross section for this collision type */
652 pcoll[k]= pt*xray_cross_section(k,&(ca[i]));
653
654 /* Total probability of ionisation */
655 ptot = 1 - (1-pcoll[ecollPHOTO])*(1-pcoll[ecollINELASTIC]);
656 if (debug && (i==0))
657 fprintf(debug,PREFIX"Ptot = %g, t = %g\n",ptot,t);
658
659 /* Check whether to ionize this guy */
660 bDOIT = FALSE;
661 switch (mode) {
662 case eionCYL:
663 bDOIT = (((rando(&ionize_seed) < ptot) && (ca[i].n < ca[i].z)) &&
664 ((sqr(x[i][XX] - hboxx) + sqr(x[i][YY] - hboxy)) < rho2));
665 break;
666 case eionSURF:
667 bDOIT = FALSE;
668 break;
669 default:
670 gmx_fatal(FARGS,"Unknown ionization mode %d (%s, line %d)",mode,
671 __FILE__,__LINE__);
672 }
673
674 if (bDOIT) {
675 clear_rvec(dv);
676
677 /* The relative probability for a photoellastic event is given by: */
678 pphot = pcoll[ecollPHOTO]/(pcoll[ecollPHOTO]+pcoll[ecollINELASTIC]);
679
680 if (rando(&ionize_seed) < pphot)
681 k = ecollPHOTO;
682 else
683 k = ecollINELASTIC;
684
685 /* If a random number is smaller than the probability for
686 * an L ionization than do that. Note that the probability
687 * may be zero (H, He), but the < instead of <= covers that.
688 */
689 nK = number_K(&ca[i]);
690 nL = number_L(&ca[i]);
691 bL = (nK == 0) || ( (nL > 0) && (rando(&ionize_seed) > prob_K(k,&(ca[i]))));
692
693 switch (k) {
694 case ecollPHOTO: {
695 /* Select which one to take by yet another random numer */
696 real theta,phi;
697
698 /* Get parameters for photoelestic effect */
699 /* Note that in the article this is called 2 theta */
700 theta = DEG2RAD*gmx_rng_gaussian_table(gaussrand)*26.0+70.0;
701
702 phi = 2*M_PI*rando(&ionize_seed);
703
704 if (bL)
705 E_lost = ephot-recoil[ca[i].z].E_L*(ca[i].n+1);
706 if (ephot == 2) E_lost += 0.5; /* Real energy should be 2.5 KeV*/
707 else {
708 E_lost = ephot-recoil[ca[i].z].E_K;
709 if (ephot == 2) E_lost += 0.5; /* Real energy should be 2.5 KeV*/
710 if ((ca[i].z > 2) && (nL > 0))
711 bKHole = TRUE;
712 }
713 if (debug)
714 fprintf(debug,"i = %d, nK = %d, nL = %d, bL = %s, bKHole = %s\n",
715 i,nK,nL,EBOOL(bL),EBOOL(bKHole));
716 if (E_lost < 0) {
717 E_lost = 0.0;
718 bIonize = FALSE;
719 bKHole = FALSE;
720 }
721 else {
722 /* Compute the components of the velocity vector */
723 factor = ((ELECTRONMASS_keV/(ATOMICMASS_keV*md->massT[i]))*
724 (SPEED_OF_LIGHT*sqrt(2*E_lost/ELECTRONMASS_keV)));
725
726 /* Subtract momentum of recoiling electron */
727 polar2cart(phi,theta,ddv);
728 for(m=0; (m<DIM); m++)
729 dv[m] -= factor*ddv[m];
730
731 if (debug)
732 pr_rvec(debug,0,"ELL",dv,DIM,TRUE);
733
734 bIonize = TRUE;
735 }
736 break;
737 }
738 case ecollINELASTIC: {
739 real theta,Ebind,Eelec;
740
741 if (bL)
742 Ebind = (ca[i].n+1)*recoil[ca[i].z].E_L;
743 else {
744 Ebind = recoil[ca[i].z].E_K;
745 if ((ca[i].z > 2) && (nL > 0))
746 bKHole = TRUE;
747 }
748 theta = DEG2RAD*rand_theta_incoh(Eindex,&ionize_seed);
749 Eelec = (sqr(ephot)/512)*(1-cos(2*theta));
750 if (ephot == 2) Eelec = (sqr(ephot+.5)/512)*(1-cos(2*theta)); /* Real energy should be 2.5 KeV*/
751 bIonize = (Ebind <= Eelec);
752 bKHole = bKHole && bIonize;
753 if (debug)
754 fprintf(debug,PREFIX"Ebind: %g, Eelectron: %g\n",Ebind,Eelec);
755 if (!bIonize) {
756 /* Subtract momentum of recoiling photon */
757 /*phi = 2*M_PI*rando(&ionize_seed);
758 bKHole = FALSE;
759 factor = ephot*438;
760 dv[XX] -= factor*cos(phi)*sin(theta);
761 dv[YY] -= factor*sin(phi)*sin(theta);
762 dv[ZZ] -= factor*cos(theta);
763 */
764 if (debug)
765 pr_rvec(debug,0,"INELL",dv,DIM,TRUE);
766 }
767 break;
768 }
769 default:
770 gmx_fatal(FARGS,"Ga direct naar de gevangenis. Ga niet langs start");
771 }
772 if (bIonize) {
773 /* First increase the charge */
774 if (ca[i].n < ca[i].z) {
775 md->chargeA[i] += 1.0;
776 md->chargeB[i] += 1.0;
777 ca[i].n++;
778 dq ++;
779 }
780 if (debug) {
781 fprintf(debug,"Random-dv[%3d] = %10.3e,%10.3e,%10.3e,"
782 " ephot = %d, Elost=%10.3e\n",
783 i,dv[XX],dv[YY],dv[ZZ],ephot,E_lost);
784 }
785 }
786 /* Now actually add the impulse to the velocities */
787 for(m=0; (m<DIM); m++)
788 v[i][m] += dv[m];
789 if (bKHole) {
790 ca[i].k ++;
791 nkhole[i]++;
792 }
793 else if (bIonize)
794 nionize[i]++;
795 }
796
797 /* Now check old event: Loop over k holes! */
798 nkh = ca[i].k;
799 for (kk = 0; (kk < nkh); kk++)
800 if (khole_decay(fp,&(ca[i]),x,v,i,&ionize_seed,ir->delta_t)) {
801 nkdecay ++;
802 ndecay[i]++;
803 md->chargeA[i] += 1.0;
804 md->chargeB[i] += 1.0;
805 }
806
807 if (debug && (ca[i].n > 0))
808 dump_ca(debug,&(ca[i]),i,__FILE__,__LINE__);
809 }
810
811 /* Sum events for statistics if necessary */
812 if (PAR(cr)) {
813 gmx_sumi(md->nr,nionize,cr);
814 gmx_sumi(md->nr,nkhole,cr);
815 gmx_sumi(md->nr,ndecay,cr);
816 nbuf[0] = dq; nbuf[1] = nkdecay;
817 gmx_sumi(2,nbuf,cr);
818 dq = nbuf[0]; nkdecay = nbuf[1];
819 }
820 /* Now sum global events on this timestep to cumulative numbers */
821 dq_tot += dq;
822 nkd_tot += nkdecay;
823
824 /* Printing time */
825 if (MASTER(cr)) {
826 /* Print data to the file that holds ionization events per atom */
827 fprintf(ion,"%12.8f",t);
828 for(i=0; (i<md->nr); i++) {
829 if (nionize[i])
830 fprintf(ion," I:%d",i+1);
831 if (nkhole[i])
832 fprintf(ion," K:%d",i+1);
833 if (ndecay[i])
834 fprintf(ion," D:%d",i+1);
835 }
836 fprintf(ion,"\n");
837 if (debug)
838 fflush(ion);
839
840 /* Print statictics to file */
841 fprintf(xvg,"%10.5f %10.3e %6d %6d %6d %6d",
842 t,pt,dq,dq_tot,nkdecay,nkd_tot);
843 fprintf(xvg,"\n");
844 if (debug)
845 fflush(xvg);
846 }
847 sfree(nionize);
848 sfree(nkhole);
849 sfree(ndecay);
850 }
851