src/kernel/ionize.c

   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