src/contrib/ehole.c

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