src/contrib/ionize.c

   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 3.0
  11  *
  12  * Copyright (c) 1991-2001
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  32  *
  33  * And Hey:
  34  * Gromacs Runs On Most of All Computer Systems
  35  */
  36 static char *SRCID_ionize_c = "$Id$";
  37 #include <string.h>
  38 #include "smalloc.h"
  39 #include "typedefs.h"
  40 #include "macros.h"
  41 #include "random.h"
  42 #include "physics.h"
  43 #include "xvgr.h"
  44 #include "vec.h"
  45 #include "pbc.h"
  46 #include "txtdump.h"
  47 #include "ionize.h"
  48 #include "names.h"
  49 #include "futil.h"
  50 #include "network.h"
  51 #include "ion_data.h"
  52
  53 #define PREFIX "IONIZE: "
  54
  55 enum { eionCYL, eionSURF, eionGAUSS, eionNR };
  56
  57 enum { ecollPHOTO, ecollINELASTIC, ecollNR };
  58
  59 typedef struct {
  60   int  z,n,k;
  61   real fj,sigPh,sigIn,vAuger;
  62 } t_cross_atom;
  63
  64 typedef struct {
  65   int nelec,maxelec,elec0,elmin_type;
  66 } t_electron_db;
  67
  68 /* BEGIN GLOBAL VARIABLES */
  69 static int   Energies[] = { 6, 8, 10, 12, 15, 20 };
  70 static int   ionize_seed = 1993;
  71 #define NENER asize(Energies)
  72
  73 static t_electron_db edb;
  74
  75 /* END GLOBAL VARIABLES */
  76
  77 void dump_ca(FILE *fp,t_cross_atom *ca,int i,char *file,int line)
  78 {
  79   fprintf(fp,PREFIX"(line %d) atom %d, z = %d, n = %d, k = %d\n",
  80           line,i,ca->z,ca->n,ca->k);
  81 }
  82
  83 t_cross_atom *mk_cross_atom(FILE *log,t_mdatoms *md,
  84                             char **atomname[],int Eindex)
  85 {
  86   int elem_index[] = { 0, 0, 0, 0, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 5, 6 };
  87   t_cross_atom *ca;
  88   int  *elemcnt;
  89   char *cc;
  90   int  i,j;
  91
  92   fprintf(log,PREFIX"Filling data structure for ionization\n");
  93   fprintf(log,PREFIX"Warning: all fj values set to 0.95 for now\n");
  94   snew(ca,md->nr);
  95   snew(elemcnt,NELEM+1);
  96   for(i=0; (i<md->nr); i++) {
  97     ca[i].n = 0;
  98     ca[i].k = 0;
  99     cc = *(atomname[i]);
 100     for(j=0; (j<NELEM); j++)
 101       if (strncmp(cc,element[j].name,strlen(element[j].name)) == 0) {
 102         ca[i].z = element[j].nel;
 103         break;
 104       }
 105     if (j == NELEM)
 106       fatal_error(0,PREFIX"Don't know number of electrons for %s",
 107                   *atomname[i]);
 108     elemcnt[j]++;
 109
 110     ca[i].sigPh = element[elem_index[ca[i].z]].cross[Eindex].photo;
 111     ca[i].sigIn = element[elem_index[ca[i].z]].cross[Eindex].incoh;
 112     ca[i].fj    = recoil[ca[i].z].Prob_K;
 113     switch (ca[i].z) {
 114     case 6:
 115       ca[i].vAuger  = 0.904;
 116       break;
 117     case 7:
 118       ca[i].vAuger  = 0.920;
 119       break;
 120     case 8:
 121       ca[i].vAuger  = 0.929;
 122       break;
 123     case 16:
 124     case 20:
 125       ca[i].vAuger = 1.0;
 126       break;
 127     default:
 128       ca[i].vAuger= -1;
 129     }
 130   }
 131
 132   fprintf(log,PREFIX"You have the following elements in your system (%d atoms):\n"PREFIX,md->nr);
 133   for(j=0; (j<NELEM); j++)
 134     if (elemcnt[j] > 0)
 135       fprintf(log,"  %s: %d",element[j].name,elemcnt[j]);
 136   fprintf(log," atoms\n");
 137
 138   sfree(elemcnt);
 139
 140   return ca;
 141 }
 142
 143 int number_K(t_cross_atom *ca)
 144 {
 145   if (ca->z <= 2)
 146     return ca->z-ca->n;
 147   else
 148     return 2-ca->k;
 149 }
 150
 151 int number_L(t_cross_atom *ca)
 152 {
 153   return ca->k-2+ca->z-ca->n;
 154 }
 155
 156 real xray_cross_section(int eColl,t_cross_atom *ca)
 157 {
 158   real c=0;
 159   int  nK,nL;
 160
 161   switch (eColl) {
 162   case ecollPHOTO:
 163     nK = number_K(ca);
 164     nL = number_L(ca);
 165     if (ca->z == 1)
 166       c = ca->sigPh;
 167     else if (ca->z == 2)
 168       c = ca->sigPh*0.5;
 169     else
 170       c = (nK*0.5*ca->fj + nL/(ca->z-2)*(1-ca->fj))*ca->sigPh;
 171     break;
 172   case ecollINELASTIC:
 173     c = (ca->z-ca->n)*ca->sigIn/ca->z;
 174     break;
 175   default:
 176     fatal_error(0,"No such collision type %d\n",eColl);
 177   }
 178   return c;
 179 }
 180
 181 real prob_K(int eColl,t_cross_atom *ca)
 182 {
 183   real Pl,Pk,P=0;
 184
 185   if ((ca->z <= 2) || (ca->z == ca->n))
 186     return 0;
 187
 188   switch (eColl) {
 189   case ecollPHOTO:
 190     Pl = (ca->k-2+ca->z-ca->n)*(1-ca->fj)/(ca->z-2);
 191     Pk = (2-ca->k)*ca->fj*0.5;
 192     P  = Pk/(Pl+Pk);
 193     break;
 194   case ecollINELASTIC:
 195     P = (2-ca->k)/(ca->z-ca->n);
 196     break;
 197   default:
 198     fatal_error(0,"No such collision type %d\n",eColl);
 199   }
 200   return P;
 201 }
 202
 203 double myexp(double x)
 204 {
 205   if (x < -70)
 206     return 0.0;
 207   else
 208     return exp(x);
 209 }
 210
 211 real ptheta_incoh(int Eindex,real theta)
 212      /* theta should be in degrees */
 213 {
 214   /* These numbers generated by fitting 5 gaussians to the real function
 215    * that describes the probability for theta.
 216    * We use symmetry in the gaussian (see 180-angle) therefore there
 217    * are fewer parameters (only 8 per energylevel).
 218    */
 219   static double ppp[NENER][8] = {
 220     { -0.00295169, 10.4847, 0.0341099, /*-*/43.1963,
 221       -0.0164054,  30.2452, 71.0311,    2.50282 },
 222     { -0.00370852, 9.02037, 0.100559,  /*-*/42.9962,
 223       -0.0537891,  35.5077, 71.4305,    1.05515 },
 224     { -0.00427039, 7.86831, 0.118042,  /*-*/45.9846,
 225       -0.0634505,  38.6134, 70.3857,    0.240082 },
 226     { -0.004514,   7.0728,  0.13464,  /*-*/48.213,
 227       -0.0723,     41.06,   69.38,     -0.02 },
 228     { -0.00488796, 5.87988, 0.159574,  /*-*/51.5556,
 229       -0.0855767,  44.7307, 69.0251,   -0.414604 },
 230     { -0.00504604, 4.56299, 0.201064,  /*-*/54.8599,
 231       -0.107153,   48.7016, 68.8212,   -0.487699 }
 232   };
 233   double g1,g2,g3,g4,g5,ptheta;
 234
 235   g1 = myexp(-0.5*sqr((theta-ppp[Eindex][7])/ppp[Eindex][1]));
 236   g2 = myexp(-0.5*sqr((theta-180+ppp[Eindex][7])/ppp[Eindex][1]));
 237   g3 = myexp(-0.5*sqr((theta-90)/ppp[Eindex][3]));
 238   g4 = myexp(-0.5*sqr((theta-ppp[Eindex][6])/ppp[Eindex][5]));
 239   g5 = myexp(-0.5*sqr((theta-180+ppp[Eindex][6])/ppp[Eindex][5]));
 240
 241   ptheta = ppp[Eindex][0]*(g1+g2) + ppp[Eindex][2]*g3 + ppp[Eindex][4]*(g4+g5);
 242
 243   return ptheta;
 244 }
 245
 246 real rand_theta_incoh(int Eindex,int *seed)
 247 {
 248 #define NINTP 450
 249 #define prev (1-cur)
 250   static bool bFirst = TRUE;
 251   static real **intp;
 252   static int  i,j,cur=1;
 253   real theta,sum,rrr,dx;
 254   real g[NENER],y[2];
 255
 256   dx = 90.0/(real)NINTP;
 257   if (bFirst) {
 258     /* Compute cumulative integrals of all probability distributions */
 259     snew(intp,NENER);
 260     for(i=0; (i<NENER); i++) {
 261       snew(intp[i],NINTP+1);
 262       y[prev]    = ptheta_incoh(i,0.0);
 263       /*sum        = y[prev];*/
 264       for(j=1; (j<=NINTP); j++) {
 265         y[cur]     = ptheta_incoh(i,j*dx);
 266         /*sum       += y[cur];*/
 267         intp[i][j] = intp[i][j-1] + (y[cur]+y[prev])*dx;
 268         cur        = prev;
 269       }
 270     }
 271     if (debug) {
 272       fprintf(debug,"Integrated probability functions for theta incoherent\n");
 273       for(j=0; (j<NINTP); j++) {
 274         fprintf(debug,"%10f",dx*j);
 275         for(i=0; (i<NENER); i++)
 276           fprintf(debug,"  %10f",intp[i][j]);
 277         fprintf(debug,"\n");
 278       }
 279     }
 280     bFirst = FALSE;
 281   }
 282
 283   rrr = rando(seed);
 284   for(j=0; (j<NINTP) && (rrr > intp[Eindex][j]); j++)
 285     ;
 286
 287   return (j-1+(rrr-intp[Eindex][j-1])/(intp[Eindex][j]-intp[Eindex][j-1]))*dx;
 288 }
 289
 290 static void polar2cart(real phi,real theta,rvec v)
 291 {
 292   v[XX] = cos(phi)*sin(theta);
 293   v[YY] = sin(phi)*sin(theta);
 294   v[ZZ] = cos(theta);
 295 }
 296
 297 void rand_vector(rvec v,int *seed)
 298 {
 299   real theta,phi;
 300
 301   theta = 180.0*rando(seed);
 302   phi   = 360.0*rando(seed);
 303   polar2cart(phi,theta,v);
 304 }
 305
 306 real electron_cross_section(FILE *fp,rvec v,real mass,int nelec)
 307 {
 308   /* Compute cross section for electrons */
 309   real T,B,U,S,Q,R,N,t,u,lnt,sigma;
 310   real a0 = 0.05292; /* nm */
 311
 312   /* Have to determine T (kinetic energy of electron) */
 313   T = 0.5*mass*iprod(v,v);
 314
 315   /* R is the binding energy of the electron in hydrogen */
 316   R = 13.61*ELECTRONVOLT;
 317
 318   /* Have to determine the binding energy B, differs per orbital of course */
 319   B = R;
 320
 321   /* Have to determine the orbital kinetic energy U */
 322   U = R;
 323
 324   /* Have to know number of electrons */
 325   N = nelec;
 326
 327   /* Magic constant Q */
 328   Q = 1;
 329
 330   /* Some help variables */
 331   t     = T/B;
 332   u     = U/B;
 333   S     = 4*M_PI*sqr(a0)*N*sqr(R/B);
 334   lnt   = log(t);
 335
 336   /* Resulting variable */
 337   sigma = (S/(t+u+1))*( 0.5*Q*lnt*(1-1/sqr(t)) + (2-Q)*(1-1/t-lnt/(t+1)) );
 338
 339   return sigma;
 340 }
 341
 342 static bool analyze_electrons(FILE *fp,t_electron_db *edb,
 343                               int natom,char **atomname[])
 344 {
 345   int  i,etp;
 346   char *cc;
 347
 348   cc = getenv("GENERATE_ELECTRONS");
 349   if (sscanf(cc,"%d",&etp) == 1) {
 350     for(i=0; (i<natom); i++) {
 351       if (strcmp(*atomname[i],"EL") == 0)
 352         break;
 353     }
 354     edb->elec0      = i;
 355     edb->nelec      = 0;
 356     edb->maxelec    = natom-i;
 357     edb->elmin_type = etp;
 358     fprintf(fp,PREFIX"There are %d possible electrons\n",edb->maxelec);
 359
 360     return TRUE;
 361   }
 362   else {
 363     fprintf(fp,PREFIX"No electron features today.\n");
 364     return FALSE;
 365   }
 366 }
 367
 368 void add_electron(FILE *fp,t_mdatoms *md,t_electron_db *edb,int ion,
 369                   rvec x[],rvec v[],rvec dv,real dt)
 370 {
 371   int  m,ee;
 372   real nv;
 373
 374   if (edb->nelec < edb->maxelec) {
 375     ee = edb->elec0+edb->nelec++;
 376     md->chargeA[ee] = md->chargeB[ee] = md->chargeT[ee] = -1;
 377     md->typeA[ee]   = md->typeB[ee]   = edb->elmin_type;
 378
 379     /* Velocity! */
 380     svmul(-md->massA[ion]*md->invmass[ee],dv,v[ee]);
 381     /* Do a first step to prevent the electron from being on top of the
 382      * nucleus, move it 0.05 nm from the nucleus
 383      */
 384     nv = 1.0/norm(v[ee]);
 385     for(m=0; (m<DIM); m++)
 386       x[ee][m] = x[ion][m] + v[ee][m]*nv*0.05;
 387   }
 388   else
 389     fatal_error(0,PREFIX"No more particles to turn into electrons\n");
 390 }
 391
 392 bool khole_decay(FILE *fp,t_cross_atom *ca,rvec x[],rvec v[],int ion,
 393                  int *seed,real dt,bool bElectron,
 394                  t_mdatoms *md,t_electron_db *edb)
 395 {
 396   rvec dv;
 397   real ndv,factor;
 398   int  m;
 399
 400   if ((ca->vAuger < 0) || (recoil[ca->z].tau == 0)) {
 401     dump_ca(stderr,ca,ion,__FILE__,__LINE__);
 402     exit(1);
 403   }
 404   if (rando(seed) < dt/recoil[ca->z].tau) {
 405     if (debug)
 406       fprintf(debug,"DECAY: Going to decay a k hole\n");
 407     ca->n++;
 408     ca->k--;
 409     /* Generate random vector */
 410     rand_vector(dv,seed);
 411
 412     factor = ca->vAuger;
 413     if (debug)
 414       fprintf(debug,"DECAY: factor=%10g, dv = (%8.3f, %8.3f, %8.3f)\n",
 415               factor,dv[XX],dv[YY],dv[ZZ]);
 416     svmul(factor,dv,dv);
 417     rvec_inc(v[ion],dv);
 418
 419     /* Now put the electron in place */
 420     if (bElectron)
 421       add_electron(fp,md,edb,ion,x,v,dv,dt);
 422
 423     return TRUE;
 424   }
 425   else
 426     return FALSE;
 427 }
 428
 429 real electron_atom_interactions(FILE *fp,t_mdatoms *md,t_inputrec *ir,
 430                                 int start,int end,
 431                                 rvec x[],rvec v[],rvec f[],matrix box)
 432 {
 433   /* Compute what the name says... */
 434   int  i,j,m,elec1,e1;
 435   rvec dx;
 436   real mindist2,vc,vtot,fscal,fc,dx2,dx_1,qi,*q;
 437
 438   mindist2 = sqr(0.05);
 439   vtot     = 0;
 440
 441   if (edb.nelec > 0) {
 442     /* Do a search... */
 443     q = md->chargeT;
 444     if (ir->ePBC != epbcNONE)
 445       init_pbc(box,FALSE);
 446     /* the end variable usually includes electrons */
 447     e1 = min(end,edb.elec0);
 448     for(i=start; (i<e1); i++) {
 449       elec1 = edb.elec0 + edb.nelec;
 450       qi = q[i]*ONE_4PI_EPS0;
 451       for(j=edb.elec0; (j<elec1); j++) {
 452         if (ir->ePBC == epbcNONE)
 453           rvec_sub(x[i],x[j],dx);
 454         else
 455           pbc_dx(x[i],x[j],dx);
 456         dx2 = iprod(dx,dx);
 457         if (dx2 < mindist2) {
 458           /* Model collision */
 459         }
 460         else {
 461           /* Do normal coulomb */
 462           dx_1  = invsqrt(dx2);
 463           vc    = qi*q[j]*dx_1;
 464           vtot += vc;
 465           fscal = vc*dx_1*dx_1;
 466           for(m=0; (m<DIM); m++) {
 467             fc       = fscal*dx[m];
 468             f[i][m] += fc;
 469             f[j][m] -= fc;
 470           }
 471         }
 472       }
 473     }
 474   }
 475   return vtot;
 476 }
 477
 478 void ionize(FILE *fp,t_mdatoms *md,char **atomname[],real t,t_inputrec *ir,
 479             rvec x[],rvec v[],int start,int end,matrix box,t_commrec *cr)
 480 {
 481   static FILE  *xvg,*ion;
 482   static char  *leg[] = { "Probability", "Primary Ionization", "Integral over PI", "KHole-Decay", "Integral over KD" };
 483   static bool  bFirst = TRUE,bElectron = FALSE;
 484   static real  t0,imax,width,inv_nratoms,rho,nphot;
 485   static real  interval;
 486   static int   dq_tot,nkd_tot,ephot,mode;
 487   static t_cross_atom *ca;
 488   static int   Eindex=-1;
 489
 490   real r,factor,ndv,E_lost=0,cross_atom,dvz,rrc;
 491   real pt,ptot,pphot,pcoll[ecollNR],tmax;
 492   real incoh,incoh_abs,sigmaPincoh,hboxx,hboxy,rho2;
 493   rvec dv,ddv;
 494   bool bIonize=FALSE,bKHole,bL,bDOIT;
 495   char *cc;
 496   int  i,j,k,kk,m,nK,nL,dq,nkh,nkdecay,elmin_type;
 497   int  *nionize,*nkhole,*ndecay,nbuf[2];
 498
 499   if (bFirst) {
 500     /* Get parameters for gaussian photon pulse from inputrec */
 501     t0       = ir->userreal1;  /* Peak of the gaussian pulse            */
 502     nphot    = ir->userreal2;  /* Intensity                             */
 503     width    = ir->userreal3;  /* Width of the peak (in time)           */
 504     rho      = ir->userreal4;  /* Diameter of the focal spot (nm)       */
 505     ionize_seed = ir->userint1;   /* Random seed for stochastic ionization */
 506     ephot    = ir->userint2;   /* Energy of the photons                 */
 507     mode     = ir->userint3;   /* Mode of ionizing                      */
 508     interval = 0.001*ir->userint4;   /* Interval between pulses (ps)    */
 509
 510     if ((width <= 0) || (nphot <= 0))
 511       fatal_error(0,"Your parameters for ionization are not set properly\n"
 512                   "width (userreal3) = %f,  nphot (userreal2) = %f",
 513                   width,nphot);
 514
 515     if ((mode < 0) || (mode >= eionNR))
 516       fatal_error(0,"Ionization mode (userint3)"
 517                   " should be in the range 0 .. %d",eionNR-1);
 518
 519     switch (mode) {
 520     case eionCYL:
 521       imax  = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
 522       break;
 523     case eionSURF:
 524       imax  = (nphot/(M_PI*sqr(rho/2)))*1e-10*1.0/(width*sqrt(2.0*M_PI));
 525       break;
 526     }
 527     if (ionize_seed == 0)
 528       ionize_seed = make_seed();
 529     if (PAR(cr)) {
 530       for(i=0; (i<cr->nodeid); i++)
 531         ionize_seed = INT_MAX*rando(&ionize_seed);
 532       fprintf(fp,PREFIX"Modifying seed on parallel processor to %d\n",
 533               ionize_seed);
 534     }
 535
 536     for(Eindex=0; (Eindex < NENER) && (Energies[Eindex] != ephot); Eindex++)
 537       ;
 538     if (Eindex == NENER)
 539       fatal_error(0,PREFIX"Energy level of %d keV not supported",ephot);
 540
 541     /* Initiate cross section data etc. */
 542     ca      = mk_cross_atom(fp,md,atomname,Eindex);
 543
 544     dq_tot  = 0;
 545     nkd_tot = 0;
 546     inv_nratoms = 1.0/md->nr;
 547
 548     xvg   = xvgropen("ionize.xvg","Ionization Events","Time (ps)","()");
 549     xvgr_legend(xvg,asize(leg),leg);
 550     ion   = ffopen("ionize.log","w");
 551
 552     fprintf(fp,PREFIX"Parameters for ionization events:\n");
 553     fprintf(fp,PREFIX"Imax = %g, t0 = %g, width = %g, seed = %d\n"
 554             PREFIX"# Photons = %g, rho = %g, ephot = %d (keV), Electrons = %s\n",
 555             imax,t0,width,ionize_seed,nphot,rho,ephot,yesno_names[bElectron]);
 556     fprintf(fp,PREFIX"Electron_mass: %10.3e(keV) Atomic_mass: %10.3e(keV)\n"
 557             PREFIX"Speed_of_light: %10.3e(nm/ps)\n",
 558             ELECTRONMASS_keV,ATOMICMASS_keV,SPEED_OF_LIGHT);
 559     fprintf(fp,PREFIX"Interval between shots: %g ps\n",interval);
 560     fprintf(fp,PREFIX"Eindex = %d\n",Eindex);
 561     fprintf(fp,PREFIX"Doing ionizations for atoms %d - %d\n",start,end);
 562
 563     bElectron = analyze_electrons(fp,&edb,md->nr,atomname);
 564
 565     fflush(fp);
 566
 567     bFirst = FALSE;
 568   }
 569
 570   /******************************************************
 571    *
 572    *    H E R E    S T A R T S   I O N I Z A T I O N
 573    *
 574    ******************************************************/
 575
 576   /* Calculate probability */
 577   tmax        = t0;
 578   if (interval > 0)
 579     while (t > (tmax+interval*0.5))
 580       tmax += interval;
 581   /*  End when t <= t0 + (N+0.5) interval */
 582
 583   pt          = imax*ir->delta_t*exp(-0.5*sqr((t-tmax)/width));
 584   dq          = 0;
 585   nkdecay     = 0;
 586
 587   hboxx       = 0.5*box[XX][XX];
 588   hboxy       = 0.5*box[YY][YY];
 589   rho2        = sqr(rho);
 590
 591   /* Width of gaussian for probability of incoherent scattering */
 592   sigmaPincoh = 1/sqrt(44.0);
 593
 594   /* Arrays for ionization statistics */
 595   snew(nionize,md->nr);
 596   snew(nkhole,md->nr);
 597   snew(ndecay,md->nr);
 598
 599   /* Loop over atoms */
 600   for(i=start; (i<end); i++) {
 601     /* Loop over collision types */
 602     bKHole = FALSE;
 603     for(k=0; (k<ecollNR); k++)
 604       /* Determine cross section for this collision type */
 605       pcoll[k]= pt*xray_cross_section(k,&(ca[i]));
 606
 607     /* Total probability of ionisation */
 608     ptot = 1 - (1-pcoll[ecollPHOTO])*(1-pcoll[ecollINELASTIC]);
 609     if (debug && (i==0))
 610       fprintf(debug,PREFIX"Ptot = %g, t = %g\n",ptot,t);
 611
 612     /* Check whether to ionize this guy */
 613     bDOIT = FALSE;
 614     switch (mode) {
 615     case eionCYL:
 616       bDOIT = (((rando(&ionize_seed) < ptot) && (ca[i].n < ca[i].z)) &&
 617                ((sqr(x[i][XX] - hboxx) + sqr(x[i][YY] - hboxy)) < rho2));
 618       break;
 619     case eionSURF:
 620       bDOIT = FALSE;
 621       break;
 622     default:
 623       fatal_error(0,"Unknown ionization mode %d (%s, line %d)",mode,
 624                   __FILE__,__LINE__);
 625     }
 626
 627     if (bDOIT) {
 628       clear_rvec(dv);
 629
 630       /* The relative probability for a photoellastic event is given by: */
 631       pphot = pcoll[ecollPHOTO]/(pcoll[ecollPHOTO]+pcoll[ecollINELASTIC]);
 632
 633       if (rando(&ionize_seed) < pphot)
 634         k = ecollPHOTO;
 635       else
 636         k = ecollINELASTIC;
 637
 638       /* If a random number is smaller than the probability for
 639        * an L ionization than do that. Note that the probability
 640        * may be zero (H, He), but the < instead of <= covers that.
 641        */
 642       nK = number_K(&ca[i]);
 643       nL = number_L(&ca[i]);
 644       bL = (nK == 0) || ( (nL > 0) && (rando(&ionize_seed) > prob_K(k,&(ca[i]))));
 645
 646       switch (k) {
 647       case ecollPHOTO: {
 648         /* Select which one to take by yet another random numer */
 649         real theta,phi;
 650
 651         /* Get parameters for photoelestic effect */
 652         /* Note that in the article this is called 2 theta */
 653         theta = DEG2RAD*gauss(70.0,26.0,&ionize_seed);
 654         phi   = 2*M_PI*rando(&ionize_seed);
 655
 656         if (bL)
 657           E_lost = ephot-recoil[ca[i].z].E_L*(ca[i].n+1);
 658         else {
 659           E_lost = ephot-recoil[ca[i].z].E_K;
 660           if ((ca[i].z > 2) && (nL > 0))
 661             bKHole = TRUE;
 662         }
 663         if (debug)
 664           fprintf(debug,"i = %d, nK = %d, nL = %d, bL = %s, bKHole = %s\n",
 665                   i,nK,nL,BOOL(bL),BOOL(bKHole));
 666         if (E_lost < 0) {
 667           E_lost  = 0.0;
 668           bIonize = FALSE;
 669           bKHole  = FALSE;
 670         }
 671         else {
 672           /* Compute the components of the velocity vector */
 673           factor = ((ELECTRONMASS_keV/(ATOMICMASS_keV*md->massT[i]))*
 674                     (SPEED_OF_LIGHT*sqrt(2*E_lost/ELECTRONMASS_keV)));
 675
 676           /* Subtract momentum of recoiling electron */
 677           polar2cart(phi,theta,ddv);
 678           for(m=0; (m<DIM); m++)
 679             dv[m] -= factor*ddv[m];
 680
 681           add_electron(fp,md,&edb,i,x,v,dv,ir->delta_t);
 682
 683           if (debug)
 684             pr_rvec(debug,0,"ELL",dv,DIM);
 685
 686           bIonize = TRUE;
 687         }
 688         break;
 689       }
 690       case ecollINELASTIC: {
 691         real theta,phi,Ebind,Eelec;
 692
 693         if (bL)
 694           Ebind = (ca[i].n+1)*recoil[ca[i].z].E_L;
 695         else {
 696           Ebind  = recoil[ca[i].z].E_K;
 697           if ((ca[i].z > 2) && (nL > 0))
 698             bKHole = TRUE;
 699         }
 700         theta      = DEG2RAD*rand_theta_incoh(Eindex,&ionize_seed);
 701         Eelec      = (sqr(ephot)/512)*(1-cos(2*theta));
 702         bIonize    = (Ebind <= Eelec);
 703         bKHole     = bKHole && bIonize;
 704         if (debug)
 705           fprintf(debug,PREFIX"Ebind: %g, Eelectron: %g\n",Ebind,Eelec);
 706         if (!bIonize) {
 707           /* Subtract momentum of recoiling photon */
 708           /*phi     = 2*M_PI*rando(&ionize_seed);
 709             bKHole  = FALSE;
 710             factor  = ephot*438;
 711             dv[XX] -= factor*cos(phi)*sin(theta);
 712             dv[YY] -= factor*sin(phi)*sin(theta);
 713             dv[ZZ] -= factor*cos(theta);
 714           */
 715           if (debug)
 716             pr_rvec(debug,0,"INELL",dv,DIM);
 717         }
 718         break;
 719       }
 720       default:
 721         fatal_error(0,"Ga direct naar de gevangenis. Ga niet langs start");
 722       }
 723       if (bIonize) {
 724         /* First increase the charge */
 725         if (ca[i].n < ca[i].z) {
 726           md->chargeA[i] += 1.0;
 727           md->chargeB[i] += 1.0;
 728           ca[i].n++;
 729           dq ++;
 730         }
 731         if (debug) {
 732           fprintf(debug,"Random-dv[%3d] = %10.3e,%10.3e,%10.3e,"
 733                   " ephot = %d, Elost=%10.3e\n",
 734                   i,dv[XX],dv[YY],dv[ZZ],ephot,E_lost);
 735         }
 736       }
 737       /* Now actually add the impulse to the velocities */
 738       for(m=0; (m<DIM); m++)
 739         v[i][m] += dv[m];
 740       if (bKHole) {
 741         ca[i].k ++;
 742         nkhole[i]++;
 743       }
 744       else if (bIonize)
 745         nionize[i]++;
 746     }
 747
 748     /* Now check old event: Loop over k holes! */
 749     nkh = ca[i].k;
 750     for (kk = 0; (kk < nkh); kk++)
 751       if (khole_decay(fp,&(ca[i]),x,v,i,&ionize_seed,ir->delta_t,
 752                       bElectron,md,&edb)) {
 753         nkdecay ++;
 754         ndecay[i]++;
 755       }
 756
 757     if (debug && (ca[i].n > 0))
 758       dump_ca(debug,&(ca[i]),i,__FILE__,__LINE__);
 759   }
 760
 761   /* Sum events for statistics if necessary */
 762   if (PAR(cr)) {
 763     gmx_sumi(md->nr,nionize,cr);
 764     gmx_sumi(md->nr,nkhole,cr);
 765     gmx_sumi(md->nr,ndecay,cr);
 766     nbuf[0] = dq; nbuf[1] = nkdecay;
 767     gmx_sumi(2,nbuf,cr);
 768     dq = nbuf[0]; nkdecay = nbuf[1];
 769   }
 770   /* Now sum global events on this timestep to cumulative numbers */
 771   dq_tot  += dq;
 772   nkd_tot += nkdecay;
 773
 774   /* Printing time */
 775   if (MASTER(cr)) {
 776     /* Print data to the file that holds ionization events per atom */
 777     fprintf(ion,"%12.8f",t);
 778     for(i=0; (i<md->nr); i++) {
 779       if (nionize[i])
 780         fprintf(ion,"  I:%d",i+1);
 781       if (nkhole[i])
 782         fprintf(ion,"  K:%d",i+1);
 783       if (ndecay[i])
 784         fprintf(ion,"  D:%d",i+1);
 785     }
 786     fprintf(ion,"\n");
 787     if (debug)
 788       fflush(ion);
 789
 790     /* Print statictics to file */
 791     fprintf(xvg,"%10.5f  %10.3e  %6d  %6d  %6d  %6d",
 792             t,pt,dq,dq_tot,nkdecay,nkd_tot);
 793     fprintf(xvg,"\n");
 794     if (debug)
 795       fflush(xvg);
 796   }
 797   sfree(nionize);
 798   sfree(nkhole);
 799   sfree(ndecay);
 800 }
 801