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