?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / local / ionize.c
blob3c8362ccea4902972b98aa8ccaccb200f8b90eb5
1 /*
2 * $Id$
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 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * Great Red Oystrich Makes All Chemists Sane
28 */
29 static char *SRCID_ionize_c = "$Id$";
30
31 #include <string.h>
32 #include "smalloc.h"
33 #include "typedefs.h"
34 #include "macros.h"
35 #include "random.h"
36 #include "physics.h"
37 #include "xvgr.h"
38 #include "vec.h"
39 #include "txtdump.h"
40 #include "ionize.h"
41 #include "names.h"
42 #include "futil.h"
43 #include "ion_data.h"
44
45 #define PREFIX "IONIZE: "
46
47 enum { eionCYL, eionSURF, eionNR };
48
49 static int Energies[] = { 6, 8, 10, 12, 15, 20 };
50 #define NENER asize(Energies)
51
52 enum { ecollPHOTO, ecollINELASTIC, ecollNR };
53
54 typedef struct {
55 int z,n,k;
56 real fj,sigPh,sigIn,vAuger;
57 } t_cross_atom;
58
59 void dump_ca(FILE *fp,t_cross_atom *ca,int i,char *file,int line)
60 {
61 fprintf(fp,PREFIX"atom %d, z = %d, n = %d, k = %d\n",i,ca->z,ca->n,ca->k);
62 }
63
64 t_cross_atom *mk_cross_atom(FILE *log,t_mdatoms *md,
65 char **atomname[],int Eindex)
66 {
67 int elem_index[] = { 0, 0, 0, 0, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 4 };
68 t_cross_atom *ca;
69 int *elemcnt;
70 char *cc;
71 int i,j;
72
73 fprintf(log,PREFIX"Filling data structure for ionization\n");
74 fprintf(log,PREFIX"Warning: all fj values set to 0.95 for now\n");
75 snew(ca,md->nr);
76 snew(elemcnt,NELEM+1);
77 for(i=0; (i<md->nr); i++) {
78 ca[i].n = 0;
79 ca[i].k = 0;
80 cc = *(atomname[i]);
81 for(j=0; (j<NELEM); j++)
82 if (strncmp(cc,element[j].name,strlen(element[j].name)) == 0) {
83 ca[i].z = element[j].nel;
84 break;
85 }
86 if (j == NELEM)
87 fatal_error(0,PREFIX"Don't know number of electrons for %s",
88 *atomname[i]);
89 elemcnt[j]++;
90
91 ca[i].sigPh = element[elem_index[ca[i].z]].cross[Eindex].photo;
92 ca[i].sigIn = element[elem_index[ca[i].z]].cross[Eindex].incoh;
93 ca[i].fj = recoil[ca[i].z].Prob_K;
94 switch (ca[i].z) {
95 case 6:
96 ca[i].vAuger = 0.904;
97 break;
98 case 7:
99 ca[i].vAuger = 0.920;
100 break;
101 case 8:
102 ca[i].vAuger = 0.929;
103 break;
104 case 16:
105 ca[i].vAuger = 1.0;
106 break;
107 default:
108 ca[i].vAuger= -1;
109 }
110 }
111
112 fprintf(log,PREFIX"You have the following elements in your system (%d atoms):\n"PREFIX,md->nr);
113 for(j=0; (j<NELEM); j++)
114 if (elemcnt[j] > 0)
115 fprintf(log," %s: %d",element[j].name,elemcnt[j]);
116 fprintf(log," atoms\n");
117
118 sfree(elemcnt);
119
120 return ca;
121 }
122
123 int number_K(t_cross_atom *ca)
124 {
125 if (ca->z <= 2)
126 return ca->z-ca->n;
127 else
128 return 2-ca->k;
129 }
130
131 int number_L(t_cross_atom *ca)
132 {
133 return ca->k-2+ca->z-ca->n;
134 }
135
136 real cross(int eColl,t_cross_atom *ca)
137 {
138 real c=0;
139 int nK,nL;
140
141 switch (eColl) {
142 case ecollPHOTO:
143 nK = number_K(ca);
144 nL = number_L(ca);
145 if (ca->z == 1)
146 c = ca->sigPh;
147 else if (ca->z == 2)
148 c = ca->sigPh*0.5;
149 else
150 c = (nK*0.5*ca->fj + nL/(ca->z-2)*(1-ca->fj))*ca->sigPh;
151 break;
152 case ecollINELASTIC:
153 c = (ca->z-ca->n)*ca->sigIn/ca->z;
154 break;
155 default:
156 fatal_error(0,"No such collision type %d\n",eColl);
157 }
158 return c;
159 }
160
161 real prob_K(int eColl,t_cross_atom *ca)
162 {
163 real Pl,Pk,P=0;
164
165 if ((ca->z <= 2) || (ca->z == ca->n))
166 return 0;
167
168 switch (eColl) {
169 case ecollPHOTO:
170 Pl = (ca->k-2+ca->z-ca->n)*(1-ca->fj)/(ca->z-2);
171 Pk = (2-ca->k)*ca->fj*0.5;
172 P = Pk/(Pl+Pk);
173 break;
174 case ecollINELASTIC:
175 P = (2-ca->k)/(ca->z-ca->n);
176 break;
177 default:
178 fatal_error(0,"No such collision type %d\n",eColl);
179 }
180 return P;
181 }
182
183 double myexp(double x)
184 {
185 if (x < -70)
186 return 0.0;
187 else
188 return exp(x);
189 }
190
191 real ptheta_incoh(int Eindex,real theta)
192 /* theta should be in degrees */
193 {
194 /* These numbers generated by fitting 5 gaussians to the real function
195 * that describes the probability for theta.
196 * We use symmetry in the gaussian (see 180-angle) therefore there
197 * are fewer parameters (only 8 per energylevel).
198 */
199 static double ppp[NENER][8] = {
200 { -0.00295169, 10.4847, 0.0341099, /*-*/43.1963,
201 -0.0164054, 30.2452, 71.0311, 2.50282 },
202 { -0.00370852, 9.02037, 0.100559, /*-*/42.9962,
203 -0.0537891, 35.5077, 71.4305, 1.05515 },
204 { -0.00427039, 7.86831, 0.118042, /*-*/45.9846,
205 -0.0634505, 38.6134, 70.3857, 0.240082 },
206 { -0.004514, 7.0728, 0.13464, /*-*/48.213,
207 -0.0723, 41.06, 69.38, -0.02 },
208 { -0.00488796, 5.87988, 0.159574, /*-*/51.5556,
209 -0.0855767, 44.7307, 69.0251, -0.414604 },
210 { -0.00504604, 4.56299, 0.201064, /*-*/54.8599,
211 -0.107153, 48.7016, 68.8212, -0.487699 }
212 };
213 double g1,g2,g3,g4,g5,ptheta;
214
215 g1 = myexp(-0.5*sqr((theta-ppp[Eindex][7])/ppp[Eindex][1]));
216 g2 = myexp(-0.5*sqr((theta-180+ppp[Eindex][7])/ppp[Eindex][1]));
217 g3 = myexp(-0.5*sqr((theta-90)/ppp[Eindex][3]));
218 g4 = myexp(-0.5*sqr((theta-ppp[Eindex][6])/ppp[Eindex][5]));
219 g5 = myexp(-0.5*sqr((theta-180+ppp[Eindex][6])/ppp[Eindex][5]));
220
221 ptheta = ppp[Eindex][0]*(g1+g2) + ppp[Eindex][2]*g3 + ppp[Eindex][4]*(g4+g5);
222
223 return ptheta;
224 }
225
226 real rand_theta_incoh(int Eindex,int *seed)
227 {
228 #define NINTP 450
229 #define prev (1-cur)
230 static bool bFirst = TRUE;
231 static real **intp;
232 static int i,j,cur=1;
233 real theta,sum,rrr,dx;
234 real g[NENER],y[2];
235
236 dx = 90.0/(real)NINTP;
237 if (bFirst) {
238 /* Compute cumulative integrals of all probability distributions */
239 snew(intp,NENER);
240 for(i=0; (i<NENER); i++) {
241 snew(intp[i],NINTP+1);
242 y[prev] = ptheta_incoh(i,0.0);
243 /*sum = y[prev];*/
244 for(j=1; (j<=NINTP); j++) {
245 y[cur] = ptheta_incoh(i,j*dx);
246 /*sum += y[cur];*/
247 intp[i][j] = intp[i][j-1] + (y[cur]+y[prev])*dx;
248 cur = prev;
249 }
250 }
251 if (debug) {
252 fprintf(debug,"Integrated probability functions for theta incoherent\n");
253 for(j=0; (j<NINTP); j++) {
254 fprintf(debug,"%10f",dx*j);
255 for(i=0; (i<NENER); i++)
256 fprintf(debug," %10f",intp[i][j]);
257 fprintf(debug,"\n");
258 }
259 }
260 bFirst = FALSE;
261 }
262
263 rrr = rando(seed);
264 for(j=0; (j<NINTP) && (rrr > intp[Eindex][j]); j++)
265 ;
266
267 return (j-1+(rrr-intp[Eindex][j-1])/(intp[Eindex][j]-intp[Eindex][j-1]))*dx;
268 }
269
270 static void polar2cart(real phi,real theta,rvec v)
271 {
272 v[XX] = cos(phi)*sin(theta);
273 v[YY] = sin(phi)*sin(theta);
274 v[ZZ] = cos(theta);
275 }
276
277 void rand_vector(rvec v,int *seed)
278 {
279 real theta,phi;
280
281 theta = 180.0*rando(seed);
282 phi = 360.0*rando(seed);
283 polar2cart(phi,theta,v);
284 }
285
286 bool khole_decay(FILE *log,t_cross_atom *ca,rvec v,int *seed,real dt,int atom)
287 {
288 rvec dv;
289 real ndv,factor;
290 int m;
291
292 if ((ca->vAuger < 0) || (recoil[ca->z].tau == 0)) {
293 dump_ca(stderr,ca,atom,__FILE__,__LINE__);
294 exit(1);
295 }
296 if (rando(seed) < dt/recoil[ca->z].tau) {
297 fprintf(log,"DECAY: Going to decay a k hole\n");
298 ca->n++;
299 ca->k--;
300 /* Generate random vector */
301 rand_vector(dv,seed);
302
303 factor = ca->vAuger;
304 fprintf(log,"DECAY: factor=%10g, dv = (%8.3f, %8.3f, %8.3f)\n",
305 factor,dv[XX],dv[YY],dv[ZZ]);
306 for(m=0; (m<DIM); m++)
307 v[m] += dv[m]*factor;
308 return TRUE;
309 }
310 else
311 return FALSE;
312 }
313
314 void ionize(FILE *log,t_mdatoms *md,char **atomname[],real t,t_inputrec *ir,
315 rvec x[],rvec v[],matrix box)
316 {
317 static FILE *xvg,*ion;
318 static char *leg[] = { "Probability", "Primary Ionization", "Integral over PI", "KHole-Decay", "Integral over KD" };
319 static bool bFirst = TRUE,bImpulse = TRUE,bExtraKinetic=FALSE;
320 static real t0,imax,width,inv_nratoms,rho,nphot,nkdecay,nkd_tot,
321 ztot,protein_radius;
322 static int seed,dq_tot,ephot,mode;
323 static t_cross_atom *ca;
324 static int Eindex=-1;
325 real r,factor,ndv,E_lost=0,cross_atom,dvz,rrc;
326 real pt,ptot,pphot,pcoll[ecollNR],mtot,delta_ekin;
327 real incoh,incoh_abs,sigmaPincoh,hboxx,hboxy,rho2;
328 rvec dv,ddv;
329 bool bIonize=FALSE,bKHole,bL,bDOIT;
330 int nK,nL;
331 char *cc;
332 int i,j,k,kk,m,dq,nkhole;
333
334 if (bFirst) {
335 /* Get parameters for gaussian photon pulse from inputrec */
336 t0 = ir->userreal1; /* Peak of the gaussian pulse */
337 nphot = ir->userreal2; /* Intensity */
338 width = ir->userreal3; /* Width of the peak (in time) */
339 rho = ir->userreal4; /* Diameter of the focal spot (nm) */
340 seed = ir->userint1; /* Random seed for stochastic ionization */
341 ephot = ir->userint2; /* Energy of the photons */
342 mode = ir->userint3; /* Mode of ionizing */
343
344 if ((width <= 0) || (nphot <= 0))
345 fatal_error(0,"Your parameters for ionization are not set properly\n"
346 "width (userreal3) = %f, nphot (userreal2) = %f",
347 width,nphot);
348
349 if ((mode < 0) || (mode >= eionNR))
350 fatal_error(0,"Ionization mode (userint3)"
351 " should be in the range 0 .. %d",eionNR-1);
352
353 switch (mode) {
354 case eionCYL:
355 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
356 break;
357 case eionSURF:
358 imax = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
359 break;
360 }
361 if (seed == 0)
362 seed = 1993;
363
364 for(Eindex=0; (Eindex < NENER) && (Energies[Eindex] != ephot); Eindex++)
365 ;
366 if (Eindex == NENER)
367 fatal_error(0,PREFIX"Energy level of %d keV not supported",ephot);
368
369 /* Initiate cross section data etc. */
370 ca = mk_cross_atom(log,md,atomname,Eindex);
371
372 dq_tot = 0;
373 nkd_tot = 0;
374 inv_nratoms = 1.0/md->nr;
375
376 bExtraKinetic = (getenv("EXTRAEKIN") != NULL);
377 bImpulse = (getenv("NOIMPULSE") == NULL);
378
379 /* compute total charge of the system */
380 ztot = 0;
381 mtot = 0;
382 for(i=0; (i<md->nr); i++) {
383 mtot += md->massA[i];
384 ztot += md->chargeA[i];
385 }
386 protein_radius = pow(mtot*AVOGADRO*1e-27*0.75/M_PI,1.0/3.0);
387
388 xvg = xvgropen("ionize.xvg","Ionization Events","Time (ps)","()");
389 xvgr_legend(xvg,asize(leg),leg);
390 ion = ffopen("ionize.log","w");
391
392 fprintf(log,PREFIX"Parameters for ionization events:\n");
393 fprintf(log,PREFIX"Imax = %g, t0 = %g, width = %g, seed = %d\n"
394 PREFIX"# Photons = %g, rho = %g, ephot = %d (keV), Impulse = %s\n",
395 imax,t0,width,seed,nphot,rho,ephot,yesno_names[bImpulse]);
396 fprintf(log,PREFIX"Electron_mass: %10.3e(keV) Atomic_mass: %10.3e(keV)\n"
397 "Speed_of_light: %10.3e(nm/ps)\n",
398 ELECTRONMASS_keV,ATOMICMASS_keV,SPEEDOFLIGHT);
399 fprintf(log,PREFIX"Eindex = %d\n",Eindex);
400 fprintf(log,PREFIX"Total charge on system: %g e. Total mass: %g u\n",
401 ztot,mtot);
402 fprintf(log,PREFIX"Estimated system radius to be %g nm\n",protein_radius);
403 fprintf(log,PREFIX
404 "Adding extra kinetic energy because of leaving electrons: %s\n",
405 bool_names[bExtraKinetic]);
406 fflush(log);
407 bFirst = FALSE;
408 }
409
410 /******************************************************
411 *
412 * H E R E S T A R T S I O N I Z A T I O N
413 *
414 ******************************************************/
415
416 /* Calculate probability */
417 pt = imax*ir->delta_t*exp(-0.5*sqr((t-t0)/width));
418 dq = 0;
419 nkdecay = 0;
420
421 hboxx = 0.5*box[XX][XX];
422 hboxy = 0.5*box[YY][YY];
423 rho2 = sqr(rho);
424
425 /* Width of gaussian for probability of incoherent scattering */
426 sigmaPincoh = 1/sqrt(44.0);
427
428 /* Print time to the file that holds ionization events per atom */
429 fprintf(ion,"%12.8f",t);
430 where();
431 /* Loop over atoms */
432 for(i=0; (i<md->nr); i++) {
433 /* Loop over collision types */
434 bKHole = FALSE;
435 for(k=0; (k<ecollNR); k++)
436 /* Determine cross section for this collision type */
437 pcoll[k]= pt*cross(k,&(ca[i]));
438
439 /* Total probability of ionisation */
440 ptot = 1 - (1-pcoll[ecollPHOTO])*(1-pcoll[ecollINELASTIC]);
441 if (debug && (i==0))
442 fprintf(debug,PREFIX"Ptot = %g, t = %g\n",ptot,t);
443
444 /* Check whether to ionize this guy */
445 bDOIT = FALSE;
446 switch (mode) {
447 case eionCYL:
448 bDOIT = (((rando(&seed) < ptot) && (ca[i].n < ca[i].z)) &&
449 ((sqr(x[i][XX] - hboxx) + sqr(x[i][YY] - hboxy)) < rho2));
450 break;
451 case eionSURF:
452 bDOIT = FALSE;
453 break;
454 default:
455 fatal_error(0,"Unknown ionization mode %d (%s, line %d)",mode,
456 __FILE__,__LINE__);
457 }
458
459 if (bDOIT) {
460 clear_rvec(dv);
461
462 /* The relative probability for a photoellastic event is given by: */
463 pphot = pcoll[ecollPHOTO]/(pcoll[ecollPHOTO]+pcoll[ecollINELASTIC]);
464
465 if (rando(&seed) < pphot)
466 k = ecollPHOTO;
467 else
468 k = ecollINELASTIC;
469
470 /* If a random number is smaller than the probability for
471 * an L ionization than do that. Note that the probability
472 * may be zero (H, He), but the < instead of <= covers that.
473 */
474 nK = number_K(&ca[i]);
475 nL = number_L(&ca[i]);
476 bL = (nK == 0) || ( (nL > 0) && (rando(&seed) > prob_K(k,&(ca[i]))));
477
478 switch (k) {
479 case ecollPHOTO: {
480 /* Select which one to take by yet another random numer */
481 real theta,phi;
482
483 /* Get parameters for photoelestic effect */
484 /* Note that in the article this is called 2 theta */
485 theta = DEG2RAD*gauss(70.0,26.0,&seed);
486 phi = 2*M_PI*rando(&seed);
487
488 if (bL)
489 E_lost = ephot-recoil[ca[i].z].E_L*(ca[i].n+1);
490 else {
491 E_lost = ephot-recoil[ca[i].z].E_K;
492 if ((ca[i].z > 2) && (nL > 0))
493 bKHole = TRUE;
494 }
495 if (debug)
496 fprintf(debug,"i = %d, nK = %d, nL = %d, bL = %s, bKHole = %s\n",
497 i,nK,nL,BOOL(bL),BOOL(bKHole));
498 if (E_lost < 0) {
499 E_lost = 0.0;
500 bIonize = FALSE;
501 bKHole = FALSE;
502 }
503 else {
504 /* Compute the components of the velocity vector */
505 factor = ((ELECTRONMASS_keV/(ATOMICMASS_keV*md->massT[i]))*
506 (SPEEDOFLIGHT*sqrt(2*E_lost/ELECTRONMASS_keV)));
507
508 /* Subtract momentum of recoiling electron */
509 polar2cart(phi,theta,ddv);
510 for(m=0; (m<DIM); m++)
511 dv[m] -= factor*ddv[m];
512
513 if (debug)
514 pr_rvec(debug,0,"ELL",dv,DIM);
515
516 bIonize = TRUE;
517 }
518 break;
519 }
520 case ecollINELASTIC: {
521 real theta,phi,Ebind,Eelec;
522
523 if (bL)
524 Ebind = (ca[i].n+1)*recoil[ca[i].z].E_L;
525 else {
526 Ebind = recoil[ca[i].z].E_K;
527 if ((ca[i].z > 2) && (nL > 0))
528 bKHole = TRUE;
529 }
530 theta = DEG2RAD*rand_theta_incoh(Eindex,&seed);
531 Eelec = (sqr(ephot)/512)*(1-cos(2*theta));
532 bIonize = (Ebind <= Eelec);
533 bKHole = bKHole && bIonize;
534 if (debug)
535 fprintf(debug,PREFIX"Ebind: %g, Eelectron: %g\n",Ebind,Eelec);
536 if (!bIonize) {
537 /* Subtract momentum of recoiling photon */
538 /*phi = 2*M_PI*rando(&seed);
539 bKHole = FALSE;
540 factor = ephot*438;
541 dv[XX] -= factor*cos(phi)*sin(theta);
542 dv[YY] -= factor*sin(phi)*sin(theta);
543 dv[ZZ] -= factor*cos(theta);
544 */
545 if (debug)
546 pr_rvec(debug,0,"INELL",dv,DIM);
547 }
548 break;
549 }
550 default:
551 fatal_error(0,"Ga direct naar de gevangenis. Ga niet langs start");
552 }
553 if (bIonize) {
554 /* First increase the charge */
555 if (ca[i].n < ca[i].z) {
556 md->chargeA[i] += 1.0;
557 md->chargeB[i] += 1.0;
558 ca[i].n++;
559 dq ++;
560 dq_tot ++;
561 }
562 if (debug) {
563 fprintf(debug,"Random-dv[%3d] = %10.3e,%10.3e,%10.3e,"
564 " ephot = %d, Elost=%10.3e\n",
565 i,dv[XX],dv[YY],dv[ZZ],ephot,E_lost);
566 }
567 }
568 /* Now actually add the impulse to the velocities */
569 if (bImpulse)
570 for(m=0; (m<DIM); m++)
571 v[i][m] += dv[m];
572 if (bKHole) {
573 ca[i].k ++;
574 fprintf(ion," K:%d",i+1);
575 }
576 else if (bIonize)
577 fprintf(ion," I:%d",i+1);
578 }
579
580 /* Now check old event: Loop over k holes! */
581 nkhole = ca[i].k;
582 for (kk = 0; (kk < nkhole); kk++)
583 if (khole_decay(log,&(ca[i]),v[i],&seed,ir->delta_t,i)) {
584 nkdecay ++;
585 nkd_tot ++;
586 fprintf(ion," D:%d",i+1);
587 }
588
589 if (debug && (ca[i].n > 0))
590 dump_ca(debug,&(ca[i]),i,__FILE__,__LINE__);
591 }
592 /* Now add random velocities corresponding to the energy depositied by
593 * the leaving electrons.
594 */
595 delta_ekin = 0;
596 if (bExtraKinetic && (dq > 0)) {
597 for(i=0; (i<dq); i++) {
598 delta_ekin += ONE_4PI_EPS0*ztot/protein_radius;
599 ztot+=1;
600 }
601 delta_ekin=fabs(delta_ekin);
602 if (debug)
603 fprintf(debug,
604 "%d leaving electrons deposited %g kJ/mol in the protein. Z=%.0f\n",
605 dq,delta_ekin,ztot);
606 delta_ekin /= md->nr;
607 for(i=0; (i<md->nr); i++) {
608 rand_vector(dv,&seed);
609 svmul(sqrt(2*delta_ekin/md->massA[i]),dv,dv);
610 rvec_inc(v[i],dv);
611 }
612 }
613
614 fprintf(ion,"\n");
615 fprintf(xvg,"%10.5f %10.3e %6d %6d %6d %6d",
616 t,pt,(int)dq,(int)dq_tot,(int)nkdecay,(int)nkd_tot);
617 if (bExtraKinetic)
618 fprintf(xvg,"%10g\n",delta_ekin);
619 else
620 fprintf(xvg,"\n");
621 fflush(ion);
622 fflush(xvg);
623 }
624
got merge a long time ago