?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Rename ffopen and ffclose to gmx_ff*.
[gromacs/AngularHB.git] / src / contrib / ehole.c
blob740546cbcef63e1e23388b7d5268327f0ef5bd71
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 }