src/contrib/ehole.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.3.3
  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-2008, 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  * Groningen Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdlib.h>
  40 #include <math.h>
  41 #include <string.h>
  42 #include "typedefs.h"
  43 #include "gromacs/utility/smalloc.h"
  44 #include "macros.h"
  45 #include "copyrite.h"
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/utility/fatalerror.h"
  48 #include "random.h"
  49 #include "gromacs/fileio/pdbio.h"
  50 #include "gromacs/utility/futil.h"
  51 #include "physics.h"
  52 #include "vec.h"
  53 #include "names.h"
  54 #include "ehdata.h"
  55
  56 typedef struct {
  57   int  maxparticle;
  58   int  maxstep;
  59   int  nsim;
  60   int  nsave;
  61   int  nana;
  62   int  seed;
  63   int  nevent;
  64   gmx_bool bForce;
  65   gmx_bool bScatter;
  66   gmx_bool bHole;
  67   real dt;
  68   real deltax;
  69   real epsr;
  70   real Alj;
  71   real Eauger;
  72   real Efermi;
  73   real Eband;
  74   real rho;
  75   real matom;
  76   real evdist;
  77   real size;
  78 } t_eh_params;
  79
  80 #define ELECTRONMASS 5.447e-4
  81 /* Resting mass of electron in a.m.u. */
  82 #define HOLEMASS (0.8*ELECTRONMASS)
  83 /* Effective mass of a hole! */
  84 #define HBAR (PLANCK/2*M_PI)
  85
  86 static void calc_forces(int n,rvec x[],rvec f[],real q[],real ener[],real Alj)
  87 {
  88   const real facel = FACEL;
  89   int  i,j,m;
  90   rvec dx;
  91   real qi,r2,r_1,r_2,fscal,df,vc,vctot,vlj,vljtot;
  92
  93   for(i=0; (i<n); i++)
  94     clear_rvec(f[i]);
  95
  96   vctot = vljtot = 0;
  97   for(i=0; (i<n-1); i++) {
  98     qi = q[i]*facel;
  99     for(j=i+1; (j<n); j++) {
 100       rvec_sub(x[i],x[j],dx);
 101       r2      = iprod(dx,dx);
 102       r_1     = 1.0/sqrt(r2);
 103       r_2     = r_1*r_1;
 104       vc      = qi*q[j]*r_1;
 105       vctot  += vc;
 106       vlj     = Alj*(r_2*r_2*r_2);
 107       vljtot += vlj;
 108       fscal   = (6*vlj+vc)*r_2;
 109       for(m=0; (m<DIM); m++) {
 110         df = fscal*dx[m];
 111         f[i][m] += df;
 112         f[j][m] -= df;
 113       }
 114     }
 115   }
 116   ener[eCOUL]   = vctot;
 117   ener[eREPULS] = vljtot;
 118   ener[ePOT]    = vctot+vljtot;
 119 }
 120
 121 static void calc_ekin(int nparticle,rvec v[],rvec vold[],
 122                       real q[],real m[],real ener[],real eparticle[])
 123 {
 124   rvec vt;
 125   real ekh=0,eke=0,ee;
 126   int  i;
 127
 128   for(i=0; (i<nparticle); i++) {
 129     rvec_add(v[i],vold[i],vt);
 130     ee = 0.125*m[i]*iprod(vt,vt);
 131     eparticle[i] = ee/ELECTRONVOLT;
 132     if (q[i] > 0)
 133       ekh += ee;
 134     else
 135       eke += ee;
 136   }
 137   ener[eHOLE]     = ekh;
 138   ener[eELECTRON] = eke;
 139   ener[eKIN]      = ekh+eke+ener[eLATTICE];
 140 }
 141
 142 static void polar2cart(real amp,real phi,real theta,rvec v)
 143 {
 144   real ss = sin(theta);
 145
 146   v[XX] = amp*cos(phi)*ss;
 147   v[YY] = amp*sin(phi)*ss;
 148   v[ZZ] = amp*cos(theta);
 149 }
 150
 151 static void rand_vector(real amp,rvec v,int *seed)
 152 {
 153   real theta,phi;
 154
 155   theta = M_PI*rando(seed);
 156   phi   = 2*M_PI*rando(seed);
 157   polar2cart(amp,phi,theta,v);
 158 }
 159
 160 static void rotate_theta(rvec v,real nv,real dth,int *seed,FILE *fp)
 161 {
 162   real   dphi,theta0,phi0,cc,ss;
 163   matrix mphi,mtheta,mphi_1,mtheta_1;
 164   rvec   vp,vq,vold;
 165
 166   copy_rvec(v,vold);
 167   theta0 = acos(v[ZZ]/nv);
 168   phi0   = atan2(v[YY],v[XX]);
 169   if (fp)
 170     fprintf(fp,"Theta = %g  Phi = %g\n",theta0,phi0);
 171
 172   clear_mat(mphi);
 173   cc = cos(-phi0);
 174   ss = sin(-phi0);
 175   mphi[XX][XX] = mphi[YY][YY] = cc;
 176   mphi[XX][YY] = -ss;
 177   mphi[YY][XX] = ss;
 178   mphi[ZZ][ZZ] = 1;
 179   m_inv(mphi,mphi_1);
 180
 181   clear_mat(mtheta);
 182   cc = cos(-theta0);
 183   ss = sin(-theta0);
 184   mtheta[XX][XX] = mtheta[ZZ][ZZ] = cc;
 185   mtheta[XX][ZZ] = ss;
 186   mtheta[ZZ][XX] = -ss;
 187   mtheta[YY][YY] = 1;
 188   m_inv(mtheta,mtheta_1);
 189
 190   dphi   = 2*M_PI*rando(seed);
 191
 192   /* Random rotation */
 193   polar2cart(nv,dphi,dth,vp);
 194
 195   mvmul(mtheta_1,vp,vq);
 196   mvmul(mphi_1,vq,v);
 197
 198   if (fp) {
 199     real cold = cos_angle(vold,v);
 200     real cnew = cos(dth);
 201     if (fabs(cold-cnew) > 1e-4)
 202       fprintf(fp,"cos(theta) = %8.4f  should be %8.4f  dth = %8.4f  dphi = %8.4f\n",
 203               cold,cnew,dth,dphi);
 204   }
 205 }
 206
 207 static int create_electron(int index,rvec x[],rvec v[],rvec vold[],rvec vv,
 208                            real m[],real q[],
 209                            rvec center,real e0,int *seed)
 210 {
 211   m[index] = ELECTRONMASS;
 212   q[index] = -1;
 213
 214   clear_rvec(v[index]);
 215   svmul(sqrt(2*e0/m[index]),vv,v[index]);
 216   copy_rvec(v[index],vold[index]);
 217   copy_rvec(center,x[index]);
 218
 219   return index+1;
 220 }
 221
 222 static int create_pair(int index,rvec x[],rvec v[],rvec vold[],
 223                        real m[],real q[],
 224                        rvec center,real e0,t_eh_params *ehp,rvec dq)
 225 {
 226   static real massfactor = 2*HOLEMASS/(ELECTRONMASS*(ELECTRONMASS+HOLEMASS));
 227   rvec x0;
 228   real ve,e1;
 229
 230   m[index]        = ELECTRONMASS;
 231   m[index+1]      = HOLEMASS;
 232   q[index]        = -1;
 233   q[index+1]      = 1;
 234
 235   rand_vector(0.5*ehp->deltax,x0,&ehp->seed);
 236   rvec_sub(center,x0,x[index]);
 237   rvec_add(center,x0,x[index+1]);
 238
 239   ve = sqrt(massfactor*e0)/(0.5*ehp->deltax);
 240   svmul(-ve,x0,v[index]);
 241   svmul(ELECTRONMASS*ve/HOLEMASS,x0,v[index+1]);
 242   copy_rvec(v[index],vold[index]);
 243   copy_rvec(v[index+1],vold[index+1]);
 244   e1 = 0.5*(m[index]*iprod(v[index],v[index])+
 245             m[index+1]*iprod(v[index+1],v[index+1]));
 246   if (fabs(e0-e1)/e0 > 1e-6)
 247     gmx_fatal(FARGS,"Error in create pair: e0 = %f, e1 = %f\n",e0,e1);
 248
 249   return index+2;
 250 }
 251
 252 static int scatter_all(FILE *fp,int nparticle,int nstep,
 253                        rvec x[],rvec v[],rvec vold[],
 254                        real mass[],real charge[],real ener[],real eparticle[],
 255                        t_eh_params *ehp,int *nelec,int *nhole,t_ana_scat s[])
 256 {
 257   int  i,m,np;
 258   real p_el,p_inel,ptot,nv,ekin,omega,theta,costheta,Q,e0,ekprime,size2,fac;
 259   rvec dq,center,vv;
 260
 261   size2 = sqr(ehp->size);
 262   np    = nparticle;
 263   for(i=0; (i<nparticle); i++) {
 264     /* Check cross sections, assume same cross sections for holes
 265      * as for electrons, for elastic scattering
 266      */
 267     if ((size2 == 0) || (iprod(x[i],x[i]) < size2)) {
 268       nv   = norm(v[i]);
 269       ekin = eparticle[i];
 270       p_el = cross_el(ekin,ehp->rho,NULL)*nv*ehp->dt;
 271
 272       /* Only electrons can scatter inelasticlly */
 273       if (charge[i] < 0)
 274         p_inel = cross_inel(ekin,ehp->rho,NULL)*nv*ehp->dt;
 275       else
 276         p_inel = 0;
 277
 278       /* Test whether we have to scatter at all */
 279       ptot = (1 - (1-p_el)*(1-p_inel));
 280       if (debug && 0)
 281         fprintf(debug,"p_el = %10.3e  p_inel = %10.3e ptot = %10.3e\n",
 282                 p_el,p_inel,ptot);
 283       if (rando(&ehp->seed) < ptot) {
 284         /* Test whether we have to scatter inelastic */
 285         ptot = p_inel/(p_el+p_inel);
 286         if (rando(&ehp->seed) < ptot) {
 287           add_scatter_event(&(s[i]),x[i],TRUE,ehp->dt*nstep,ekin);
 288           /* Energy loss in inelastic collision is omega */
 289           if ((omega = get_omega(ekin,&ehp->seed,debug,NULL)) >= ekin)
 290             gmx_fatal(FARGS,"Energy transfer error: omega = %f, ekin = %f",
 291                         omega,ekin);
 292           else {
 293             /* Scattering angle depends on energy and energy loss */
 294             Q = get_q_inel(ekin,omega,&ehp->seed,debug,NULL);
 295             costheta = -0.5*(Q+omega-2*ekin)/sqrt(ekin*(ekin-omega));
 296
 297             /* See whether we have gained enough energy to liberate another
 298              * hole-electron pair
 299              */
 300             e0      = band_ener(&ehp->seed,debug,NULL);
 301             ekprime = e0 + omega - (ehp->Efermi+0.5*ehp->Eband);
 302             /* Ouput */
 303             fprintf(fp,"Inelastic %d: Ekin=%.2f Omega=%.2f Q=%.2f Eband=%.2f costheta=%.3f\n",
 304                     i+1,ekin,omega,Q,e0,costheta);
 305             if ((costheta < -1) || (costheta > 1)) {
 306               fprintf(fp,"Electron/hole creation not possible due to momentum constraints\n");
 307               /* Scale the velocity according to the energy loss */
 308               svmul(sqrt(1-omega/ekin),v[i],v[i]);
 309               ener[eLATTICE] += omega*ELECTRONVOLT;
 310             }
 311             else {
 312               theta = acos(costheta);
 313
 314               copy_rvec(v[i],dq);
 315               /* Rotate around theta with random delta phi */
 316               rotate_theta(v[i],nv,theta,&ehp->seed,debug);
 317               /* Scale the velocity according to the energy loss */
 318               svmul(sqrt(1-omega/ekin),v[i],v[i]);
 319               rvec_dec(dq,v[i]);
 320
 321               if (ekprime > 0) {
 322                 if (np >= ehp->maxparticle-2)
 323                   gmx_fatal(FARGS,"Increase -maxparticle flag to more than %d",
 324                               ehp->maxparticle);
 325                 if (ehp->bHole) {
 326                   np = create_pair(np,x,v,vold,mass,charge,x[i],
 327                                    ekprime*ELECTRONVOLT,ehp,dq);
 328                   (*nhole)++;
 329                 }
 330                 else {
 331                   copy_rvec(x[i],center);
 332                   center[ZZ] += ehp->deltax;
 333                   rand_vector(1,vv,&ehp->seed);
 334                   np = create_electron(np,x,v,vold,vv,mass,charge,
 335                                        x[i],ekprime*ELECTRONVOLT,&ehp->seed);
 336                 }
 337                 ener[eLATTICE] += (omega-ekprime)*ELECTRONVOLT;
 338                 (*nelec)++;
 339               }
 340               else
 341                 ener[eLATTICE] += omega*ELECTRONVOLT;
 342             }
 343           }
 344         }
 345         else {
 346           add_scatter_event(&(s[i]),x[i],FALSE,ehp->dt*nstep,ekin);
 347           if (debug)
 348             fprintf(debug,"Elastic scattering event\n");
 349
 350           /* Scattering angle depends on energy only */
 351           theta = get_theta_el(ekin,&ehp->seed,debug,NULL);
 352           /* Rotate around theta with random delta phi */
 353           rotate_theta(v[i],nv,theta,&ehp->seed,debug);
 354         }
 355       }
 356     }
 357   }
 358   return np;
 359 }
 360
 361 static void integrate_velocities(int nparticle,rvec vcur[],rvec vnext[],
 362                                  rvec f[],real m[],real dt)
 363 {
 364   int i,k;
 365
 366   for(i=0; (i<nparticle); i++)
 367     for(k=0; (k<DIM); k++)
 368       vnext[i][k] = vcur[i][k] + f[i][k]*dt/m[i];
 369 }
 370
 371 static void integrate_positions(int nparticle,rvec x[],rvec v[],real dt)
 372 {
 373   int i,k;
 374
 375   for(i=0; (i<nparticle); i++)
 376     for(k=0; (k<DIM); k++)
 377       x[i][k] += v[i][k]*dt;
 378 }
 379
 380 static void print_header(FILE *fp,t_eh_params *ehp)
 381 {
 382   fprintf(fp,"Welcome to the electron-hole simulation!\n");
 383   fprintf(fp,"The energies printed in this file are in eV\n");
 384   fprintf(fp,"Coordinates are in nm because of fixed width format\n");
 385   fprintf(fp,"Atomtypes are used for coloring in rasmol\n");
 386   fprintf(fp,"O: electrons (red), N: holes (blue)\n");
 387   fprintf(fp,"Parametes for this simulation\n");
 388   fprintf(fp,"seed = %d maxstep = %d dt = %g\n",
 389           ehp->seed,ehp->maxstep,ehp->dt);
 390   fprintf(fp,"nsave = %d nana = %d Force = %s Scatter = %s Hole = %s\n",
 391           ehp->nsave,ehp->nana,gmx_bool_names[ehp->bForce],
 392           gmx_bool_names[ehp->bScatter],gmx_bool_names[ehp->bHole]);
 393   if (ehp->bForce)
 394     fprintf(fp,"Force constant for repulsion Alj = %g\n",ehp->Alj);
 395 }
 396
 397 static void do_sim(FILE *fp,char *pdbfn,t_eh_params *ehp,
 398                    int *nelec,int *nhole,t_ana_struct *total,
 399                    t_histo *hmfp,t_ana_ener *ae,int serial)
 400 {
 401   FILE         *efp;
 402   int          nparticle[2];
 403   rvec         *x,*v[2],*f,center,vv;
 404   real         *charge,*mass,*ener,*eparticle;
 405   t_ana_struct *ana_struct;
 406   t_ana_scat   *ana_scat;
 407   int          step,i,cur = 0;
 408 #define next (1-cur)
 409
 410   /* Open output file */
 411   fprintf(fp,"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
 412   fprintf(fp,"Simulation %d/%d\n",serial+1,ehp->nsim);
 413
 414   ana_struct = init_ana_struct(ehp->maxstep,ehp->nana,ehp->dt,
 415                                ehp->maxparticle);
 416   /* Initiate arrays. The charge array determines whether a particle is
 417    * a hole (+1) or an electron (-1)
 418    */
 419   snew(x,ehp->maxparticle);          /* Position  */
 420   snew(v[0],ehp->maxparticle);       /* Velocity  */
 421   snew(v[1],ehp->maxparticle);       /* Velocity  */
 422   snew(f,ehp->maxparticle);          /* Force     */
 423   snew(charge,ehp->maxparticle);     /* Charge    */
 424   snew(mass,ehp->maxparticle);       /* Mass      */
 425   snew(eparticle,ehp->maxparticle);  /* Energy per particle */
 426   snew(ana_scat,ehp->maxparticle);   /* Scattering event statistics */
 427   snew(ener,eNR);                    /* Eenergies */
 428
 429   clear_rvec(center);
 430   /* Use first atom as center, it has coordinate 0,0,0 */
 431   if (ehp->bScatter) {
 432     /* Start with an Auger electron */
 433     nparticle[cur]=0;
 434     for(i=0; (i<ehp->nevent); i++) {
 435       if (ehp->nevent == 1) {
 436         clear_rvec(vv);
 437         vv[ZZ] = 1;
 438       }
 439       else
 440         rand_vector(1,vv,&ehp->seed);
 441       nparticle[cur]  = create_electron(nparticle[cur],x,v[cur],v[next],
 442                                         vv,mass,charge,center,
 443                                         ehp->Eauger*ELECTRONVOLT,&ehp->seed);
 444       rand_vector(ehp->evdist*0.1,vv,&ehp->seed);
 445       rvec_inc(center,vv);
 446     }
 447   }
 448   else if (ehp->bForce) {
 449     /* Start with two electron and hole pairs */
 450     nparticle[cur]  = create_pair(0,x,v[cur],v[next],mass,charge,center,
 451                                   0.2*ehp->Eauger*ELECTRONVOLT,ehp,center);
 452     center[ZZ] = 0.5; /* nm */
 453     (*nelec)++;
 454     (*nhole)++;
 455   }
 456   else {
 457     fprintf(fp,"Nothing to do. Doei.\n");
 458     return;
 459   }
 460   nparticle[next] = nparticle[cur];
 461   for(step=0; (step<=ehp->maxstep); step++) {
 462     if (ehp->bScatter)
 463       nparticle[next] = scatter_all(fp,nparticle[cur],step,x,v[cur],v[next],
 464                                     mass,charge,ener,eparticle,ehp,
 465                                     nelec,nhole,ana_scat);
 466
 467     if (ehp->bForce)
 468       calc_forces(nparticle[cur],x,f,charge,ener,ehp->Alj);
 469
 470     integrate_velocities(nparticle[next],v[cur],v[next],f,mass,ehp->dt);
 471
 472     calc_ekin(nparticle[next],v[cur],v[next],charge,mass,ener,eparticle);
 473     ener[eTOT] = ener[eKIN] + ener[ePOT];
 474
 475     /* Produce output whenever the user says so, or when new
 476      * particles have been created.
 477      */
 478     if ((step == ehp->maxstep) ||
 479         ((ehp->nana != 0) && ((step % ehp->nana) == 0))) {
 480       analyse_structure(ana_struct,(step*ehp->dt),center,x,
 481                         nparticle[next],charge);
 482       add_ana_ener(ae,(step/ehp->nana),ener);
 483     }
 484     cur = next;
 485
 486     integrate_positions(nparticle[cur],x,v[cur],ehp->dt);
 487   }
 488   for(i=0; (i<nparticle[cur]); i++) {
 489     analyse_scatter(&(ana_scat[i]),hmfp);
 490     done_scatter(&(ana_scat[i]));
 491   }
 492   sfree(ener);
 493   sfree(ana_scat);
 494   sfree(eparticle);
 495   sfree(mass);
 496   sfree(charge);
 497   sfree(f);
 498   sfree(v[1]);
 499   sfree(v[0]);
 500   sfree(x);
 501   dump_as_pdb(pdbfn,ana_struct);
 502   add_ana_struct(total,ana_struct);
 503   done_ana_struct(ana_struct);
 504   sfree(ana_struct);
 505 }
 506
 507 void do_sims(int NFILE,t_filenm fnm[],t_eh_params *ehp)
 508 {
 509   t_ana_struct *total;
 510   t_ana_ener   *ae;
 511   t_histo      *helec,*hmfp;
 512   int          *nelectron;
 513   int          i,imax,ne,nh;
 514   real         aver;
 515   FILE         *fp,*logfp;
 516   char         *pdbbuf,*ptr,*rptr;
 517
 518   ptr  = ftp2fn(efPDB,NFILE,fnm);
 519   rptr = strdup(ptr);
 520   if ((ptr = strstr(rptr,".pdb")) != NULL)
 521     *ptr = '\0';
 522   snew(pdbbuf,strlen(rptr)+10);
 523
 524   total = init_ana_struct(ehp->maxstep,ehp->nana,ehp->dt,1);
 525   hmfp  = init_histo((int)ehp->Eauger,0,(int)ehp->Eauger);
 526   helec = init_histo(500,0,500);
 527   snew(ae,1);
 528
 529   logfp = gmx_ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
 530   print_header(logfp,ehp);
 531
 532   for(i=0; (i<ehp->nsim); i++) {
 533     nh = ne = 0;
 534     sprintf(pdbbuf,"%s-%d.pdb",rptr,i+1);
 535     do_sim(logfp,pdbbuf,ehp,&ne,&nh,total,hmfp,ae,i);
 536     add_histo(helec,ne,1);
 537     fprintf(stderr,"\rSim: %d/%d",i+1,ehp->nsim);
 538   }
 539   fprintf(stderr,"\n");
 540   gmx_ffclose(logfp);
 541
 542   sfree(rptr);
 543   sfree(pdbbuf);
 544   dump_ana_struct(opt2fn("-maxdist",NFILE,fnm),opt2fn("-nion",NFILE,fnm),
 545                   opt2fn("-gyr_com",NFILE,fnm),opt2fn("-gyr_origin",NFILE,fnm),
 546                   total,ehp->nsim);
 547   dump_ana_ener(ae,ehp->nsim,ehp->dt*ehp->nana,
 548                 opt2fn("-ener",NFILE,fnm),total);
 549   done_ana_struct(total);
 550
 551   dump_histo(helec,opt2fn("-histo",NFILE,fnm),
 552              "Number of cascade electrons","N","",enormFAC,1.0/ehp->nsim);
 553   dump_histo(hmfp,opt2fn("-mfp",NFILE,fnm),
 554              "Mean Free Path","Ekin (eV)","MFP (nm)",enormNP,1.0);
 555 }
 556
 557 int main(int argc,char *argv[])
 558 {
 559   const char *desc[] = {
 560     "[TT]ehole[tt] performs a molecular dynamics simulation of electrons and holes",
 561     "in an implicit lattice. The lattice is modeled through scattering cross",
 562     "sections, for elastic and inelastic scattering.",
 563     "A detailed description of the scatterning processes simulated in ehole",
 564     "can be found in Timneanu et al. Chemical Physics 299 (2004) 277-283",
 565     "The paper also includes a description how to calculate the input files.[PAR]",
 566     "Description of the input files for [TT]ehole[tt]:[BR]",
 567     "[TT]-sigel.dat[tt]: elastic cross section (per atom). Two columns: Impact electron energy (eV) vs Elastic cross section (A2).[BR]",
 568     "[TT]-siginel.dat[tt]: inelastic cross section (per atom). Two columns: Impact electron energy (eV) vs Inelastic cross section (A2).[BR]",
 569     "[TT]-band-ener.dat[tt]: Probability of finding an electron in the valence band.",
 570     "Two columns: Impact electron energy (eV) vs Probability[BR]",
 571     "[TT]-eloss.dat[tt]: Probability of energy loss due to inelastic scattering. Three columns: Impact electron energy (eV) vs  Integrated probability vs Energy loss in inelastic scattering (eV).[BR]",
 572     "[TT]-theta-el.dat[tt]: Probability of elastic scattering angle. Three columns: Impact electron energy (eV) vs Integrated probability vs Scattering angle (rad).[BR]",
 573     "[TT]-qtrans.dat[tt]: Four columns: Impact electron energy (eV) vs Inelastic energy loss (eV) vs Integrated probability vs Scattering angle (rad).[PAR]",
 574     "The program produces a number of output files. It is important that",
 575     "the actual content is well-defined, sucht that no misunderstanding can",
 576     "occur (famous last words...). Anyway, the program does a number of",
 577     "simulations, and averages results over these. Here is a list of each of",
 578     "the results and how they are computed:[BR]",
 579     "[TT]-histo[tt] Distribution of nuber of liberated secondary electrons per simulation.[BR]",
 580     "[TT]-maxdist[tt] The maximum distance from the origin that any electron in any simulation reaches.[BR]",
 581     "[TT]-gyr_com[tt] The radius of gyration of the electron cloud with respect to its center of mass (contains 4 columns).[BR]",
 582     "[TT]-gyr_origin[tt] The radius of gyration of the electron cloud with respect to the origin (contains 4 columns).[BR]",
 583     "[TT]-mfp[tt] The mean free path of the electrons as a function of energy. If this is not a smooth curve you need to increase the number of simulations.[BR]",
 584     "[TT]-nion[tt] The number of ions as a function of time, averaged over simulations.[BR]",
 585     "[TT]-ener[tt] The energy terms in the simulation (note that there are multiple columns, so use [TT]xmgrace -nxy[tt]). This shows important information about the stability of the simulation, that is the total energy should be conserved. In this figure you can also inspect the kinetic energy per electron in order to check whether the electrons have thermalized.[BR]"
 586   };
 587   static t_eh_params ehp = {
 588     100,    /* Max number of particles. Is a parameter but should be dynamic */
 589     100000, /* Number of integration steps */
 590     1,      /* nsave */
 591     1,      /* nana */
 592     1,      /* nsim */
 593     1993,   /* Random seed */
 594     1,      /* Number of events */
 595     FALSE,  /* Use forces */
 596     TRUE,   /* Use scattering */
 597     FALSE,  /* Creat holes */
 598     1e-5,   /* Time step */
 599     0.05,   /* Distance (nm) between electron and hole when creating them */
 600     1.0,    /* Dielectric constant */
 601     0.1,    /* Force constant for repulsion function */
 602     250,    /* Starting energy for the first Auger electron */
 603     28.7,   /* Fermi level (eV) of diamond. */
 604     5.46,   /* Band gap energy (eV) of diamond */
 605     3.51,   /* Density of the solid */
 606     12.011, /* (Average) mass of the atom */
 607     10000.0,/* Distance between events */
 608     0.0     /* Size of the system */
 609   };
 610   static gmx_bool bTest    = FALSE;
 611   t_pargs pa[] = {
 612     { "-maxparticle", FALSE, etINT,  {&ehp.maxparticle},
 613       "Maximum number of particles" },
 614     { "-maxstep",     FALSE, etINT,  {&ehp.maxstep},
 615       "Number of integration steps" },
 616     { "-nsim",        FALSE, etINT,  {&ehp.nsim},
 617       "Number of independent simulations writing to different output files" },
 618     { "-nsave",       FALSE, etINT,  {&ehp.nsave},
 619       "Number of steps after which to save output. 0 means only when particles created. Final step is always written." },
 620     { "-nana",        FALSE, etINT,  {&ehp.nana},
 621       "Number of steps after which to do analysis." },
 622     { "-seed",        FALSE, etINT,  {&ehp.seed},
 623       "Random seed" },
 624     { "-dt",          FALSE, etREAL, {&ehp.dt},
 625       "Integration time step (ps)" },
 626     { "-rho",         FALSE, etREAL, {&ehp.rho},
 627       "Density of the sample (kg/l). Default is for Diamond" },
 628     { "-matom",       FALSE, etREAL, {&ehp.matom},
 629       "Mass (a.m.u.) of the atom in the solid. Default is C" },
 630     { "-fermi",       FALSE, etREAL, {&ehp.Efermi},
 631       "Fermi energy (eV) of the sample. Default is for Diamond" },
 632     { "-gap",         FALSE, etREAL, {&ehp.Eband},
 633       "Band gap energy (eV) of the sample. Default is for Diamond" },
 634     { "-auger",       FALSE, etREAL, {&ehp.Eauger},
 635       "Impact energy (eV) of first electron" },
 636     { "-dx",          FALSE, etREAL, {&ehp.deltax},
 637       "Distance between electron and hole when creating a pair" },
 638     { "-test",        FALSE, etBOOL, {&bTest},
 639       "Test table aspects of the program rather than running it for real" },
 640     { "-force",       FALSE, etBOOL, {&ehp.bForce},
 641       "Apply Coulomb/Repulsion forces" },
 642     { "-hole",        FALSE, etBOOL, {&ehp.bHole},
 643       "Create electron-hole pairs rather than electrons only" },
 644     { "-scatter",     FALSE, etBOOL, {&ehp.bScatter},
 645       "Do the scattering events" },
 646     { "-nevent",      FALSE, etINT,  {&ehp.nevent},
 647       "Number of initial Auger electrons" },
 648     { "-evdist",      FALSE, etREAL, {&ehp.evdist},
 649       "Average distance (A) between initial electronss" },
 650     { "-size",        FALSE, etREAL, {&ehp.size},
 651       "Size of the spherical system. If 0, then it is infinite" }
 652   };
 653 #define NPA asize(pa)
 654   t_filenm fnm[] = {
 655     { efLOG, "-g",          "ehole",      ffWRITE },
 656     { efDAT, "-sigel",      "sigel",      ffREAD },
 657     { efDAT, "-sigin",      "siginel",    ffREAD },
 658     { efDAT, "-eloss",      "eloss",      ffREAD },
 659     { efDAT, "-qtrans",     "qtrans",     ffREAD },
 660     { efDAT, "-band",       "band-ener",  ffREAD },
 661     { efDAT, "-thetael",    "theta-el",   ffREAD },
 662     { efPDB, "-o",          "ehole",      ffWRITE },
 663     { efXVG, "-histo",      "histo",      ffWRITE },
 664     { efXVG, "-maxdist",    "maxdist",    ffWRITE },
 665     { efXVG, "-gyr_com",    "gyr_com",    ffWRITE },
 666     { efXVG, "-gyr_origin", "gyr_origin", ffWRITE },
 667     { efXVG, "-mfp",        "mfp",        ffWRITE },
 668     { efXVG, "-nion",       "nion",       ffWRITE },
 669     { efXVG, "-ener",       "ener",       ffWRITE }
 670   };
 671 #define NFILE asize(fnm)
 672   int seed;
 673
 674   CopyRight(stdout,argv[0]);
 675   parse_common_args(&argc,argv,PCA_BE_NICE,NFILE,fnm,
 676                     NPA,pa,asize(desc),desc,0,NULL);
 677   please_cite(stdout,"Timneanu2004a");
 678
 679   if (ehp.deltax <= 0)
 680     gmx_fatal(FARGS,"Delta X should be > 0");
 681   ehp.Alj = FACEL*pow(ehp.deltax,5);
 682   ehp.rho = (ehp.rho/ehp.matom)*AVOGADRO*1e-21;
 683
 684   init_tables(NFILE,fnm,ehp.rho);
 685
 686   if (bTest)
 687     test_tables(&ehp.seed,ftp2fn(efPDB,NFILE,fnm),ehp.rho);
 688   else
 689     do_sims(NFILE,fnm,&ehp);
 690
 691   gmx_thanx(stdout);
 692
 693   return 0;
 694 }