?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Remove some boost::scoped_ptr uses
[gromacs/AngularHB.git] / src / contrib / ehole.c
blobd5510ecb5d3a3fc88d1b06d2ee6b3aac84ae1431
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 "gromacs/math/units.h"
52 #include "gromacs/math/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,0,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 }