src/contrib/ehole.c

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