Accelerator Laboratory, University
of Helsinki, P.O. BOX 43, FIN-00014 Helsinki, Finland
(email kai.nordlund@helsinki.fi)
September 1, 1994
Last update: Jarkko Peltola and Kai Nordlund, August 2002.
It describes the programming principles and use of the program, not the physical principles underlying it. The physics is presented in the publication titled Molecular dynamics simulations of ion ranges in the 1-100 keV energy range (Comp. Mat. Sci. 3 (1995) 448-). The reader is encouraged to read the publication first if he wishes to fully understand the forthcoming chapters in this document. A preprint of the paper is available in postscript format (with figures and everything) in here. Related information is available in my PhD thesis and publications.
The program source code files and sample input and output files can be viewed by selecting the links in this document. NOTE: the version here (V1.01A) is not the newest one, and contains a few minor bugs which caused problems on some systems and in some physical cases. Newer versions (where all known bugs have been corrected and a few features added) can be obtained from me.
mdh now has a user interface, mdsetup, which generates the input files needed. mdsetup has a small compound library and has default settings for those parameters for which sensible defaults can be assigned. Bear in mind the perennial caveat about scientific programs in general: garbage in - garbage out. If you attempt to use the program without understanding the physics underlying it, you may obtain more or less insane results, or simply get the program to tilt in some exotic way. However, mdsetup combined with a basic understanding of the physics concerned should be enough to take you through (if you survived TRIM, you are likely to live through mdh ).
The official name of the program is MDRANGE. However, in the actual program files the shorter, more convenient name mdh (abbreviated from Molecular Dynamics High-energy) is used. Both names therefore (at least for now) mean exactly the same program.
The program is a molecular dynamics (MD) simulation program tailored for effective calculation of ion ranges. The word effective used here must be understood in the context of high-energy molecular dynamics calculations; if such a simulation only needs a few hours of CPU time to yield useful results, it can be considered very effective.
Versions prior to 1.8 may have been somewhat inaccurate for collision energies higher than about 100 keV/amu due to potential interpolation problems at very high energies.
Details on the tests can be obtained from the authors.
When compiling in a new environment, the user should set the appropriate compilation options in the file config according to instructions given therein.
If problems are encountered during compilation, the error most likely lies in wrong options in the config file, the script writeheader, or the subroutine eltime() in the file output.c, the only subroutine that contains non-ANSI features.
In longer-term use, it is usually most convenient to start by mdsetup, and after this change the few parameters which need to be changed parameters in the input file param.in by hand.
It works just like mdh, but for the calculation of the range profile it uses a reed-like method. To turn on the reed calculation just use the option
reccalc->reed:= 1and use the executable mdhreed
Note that the other output files (deposited energy, recoil spectra, etc.) are not enhanced by the reed method. So if you are primarily interested in those, using mdhreed will just slow things down, and it is better to use plain mdh.
This package includes mdsetup, a user interface to mdh. You can use mdsetup to create the input files coords.in, elstop.in and param.in for the most common cases of interest. However, mdh does have many options which can not be generated by mdsetup, but need to be handled by editing param.in manually.
mdsetup uses zbl96 to create the elstop - if you have a different version of zbl, just change the content of the string variable (zbl) in mdsetup.c and recompile.
The repulsive potential can be calculated using quantum mechanical software, for example DMol (See notes on using repulsive potential). If you have no suitable programs available, select the default potential (ZBL universal repulsive potential. mdsetup gives you a chance to change the parameters).
All of the input files are of ASCII-format; you can edit them easily.
mdsetup has default values for all parameters for which sensible, "universal" default values can be given. The default value is shown when the value of parameter is requested, e.g.
Please enter the atomic number or chemical symbol of the grid atoms: 22 Ti , Z=22 Is this OK (Y/N) [Y] :has the default value Y. If you are satisfied with the default value, just press enter. mdsetup also has online help where this has been deemed necessary, e.g.
Please enter the average size of the domains (Å, ? for help):
Start mdsetup and select the user level. If you choose to be an expert user, the range of settings you will be asked for is much wider (but you are still offered the defaults). The easiest and fastest way to use mdh is to answer No. First, the coordinates files are created. Next MDSetup creates the electronic stopping file. The electronic stopping is weighted average of the stopping powers of all the elements in the target (averaged over every atom in every layer).
The next step is the creation of a parameter file, which contains most of the settings, e.g. the range of starting angles and positions.To save your time, some of the seldom used parameters have been omitted. You are now asked if you want to create a script. By creating a script you can select the input files used for the run and the output generated. The script is saved; you can either run it immediately (with an adjustable nice level) or later. As in creating the parameter file, some of the seldom used options have been omitted; enter mdh -h for a complete list of options.
If you run the series of simulations using the same substances, e.g. He-Mo with energies 25, 50, 75 and 100 keV, the same reppot.in and coords.in files can be used for all the simulations. As only one line, the energy of the projectile (reccalc->E0) is changed, it's sensible to create the param.in file for the first simulation using mdsetup and edit the file using any ASCII-editor for the subsequent simulations (changing the seed of the random number generator, gen->seed, is also recommended). The elstop file must contain the electronic stopping at sufficiently large velocities - therefore, you might use a 100 keV elstop for an 80 keV simulation but not vice versa.
Apart from outright debugging the debug modes can also be used to obtain a number of physical quantities. For instance, the temperature in the cell as a function of time in the initial displacement calculation can be obtained with the Unix command sequence "mdh -dm |& grep ^ini"
With the options "-r", "-p" and "-pi" output of the recoil/ion or of all atom coordinates during runtime can be selected. By piping the output of mdh to the drawing program dpc this enables making simple on-line animations of atom movement in the program.
There is unfortunately quite some confusion in the use of the terms "recoil" and "primary/secondary recoil spectrum" in the code. Frequently the energetic ion whose motion is followed is called the "recoil" in the code, without making the distinction whether this is an ion starting from outside the sample, or a recoiling atom starting from inside. This is primarily to avoid having to use the term "recoil/ion" all along (which, besides, isn't even possible in variable names). Also, the atoms receiving energy from the energetic ion are sometimes called primary and sometimes secondary recoils. In later versions of the code I try to always use "primary recoil spectrum" but the old term may still hang around in old versions :-(
The basic result output of the program consists of the file range.out, which contains the ion range distribution (in the units "number of atoms") as a function of the depth in Ångströms. The MDRANGE3.0 version produced also a file realrange.out, which contains the range distribution as atom density vs. depth. This is handy, when comparing to the experiments and the experimental fluence is known (ions/cm^2). The fluence is given in the param.in file: physical->dose:= 1e14 (for example).
If the "-FD" command line option is selected the deposited energies are also output to the file depen.out. The first column in this file gives the depth in Ångströms, the second and third the nuclear and electronic deposited energies in units of eV/Å.
To be exact, the angle is chosen using the parameters
reccalc->Fii0:= 0.0 # Azimuth angle of the recoil atom reccalc->Fiimax:= 0.0 # is chosen randomly from [Fii0,Fiimax] reccalc->Theta0:= 8.0 # Polar angle of the recoil atom reccalc->Thetamax:= 8.0 # is chosen randomly from [Theta0,Thetamax]and mdh selects the angle randomly between the min and max. If Theta0 == 0 the theta angle is weighted to give equal weighting of spatial angles.
The exact algorithm used to determine the velocity components (vx,vy,vz) for a given velocity v, theta and phi is:
vx=v*sin(theta)*cos(fii) vy=v*sin(theta)*sin(fii) vz=v*cos(theta)Thus e.g. Fii0=Fiimax=0 and Theta0=Thetamax=0 chooses implantation straight into the 001 channel.
type[0].Z:= 36 type[0].m:= 86.0 type[1].Z:= 29 type[1].m:= 63.546 type[1].prob:= 0.5 type[2].Z:= 47 type[2].m:= 107.87 type[2].prob:= 0.5Of course you have to remember to set appropriate reppots, electronic stopping and Debye temperature.
The attractive potential is only used for getting thermal displacements, and is not used at all in the actual range (recoil event) calculation. This justifies the use of the simple two-body potentials; as long as they yield a stabile lattice and realistic thermal displacements, they should be good enough in this application.
The default setting in mdsetup is to use the Debye method for elements and an attractive potential of the user's choice for compounds. The potential parameters can be modified in the param.in file. For instance, by choosing the potential parameters so that they correspond to the lattice constant, DMol potential and Debye thermal displacement, we have gotten a modified Mazzone potential for InP. It can be used by selecting the Mazzone potential and adding the following parameters in the param.in:
pot->attr.mazD[1][1]:= 1.9838704 pot->attr.mazalpha[1][1]:= 1.0814278 pot->attr.mazr1[1][1]:= 2.8865101 pot->attr.mazK[1][1]:= 1.800 pot->attr.mazr2[1][1]:= 4.150 pot->attr.mazd[1][1]:= 0.330 pot->attr.mazD[1][2]:= 2.5553605 pot->attr.mazalpha[1][2]:= 1.1331184 pot->attr.mazr1[1][2]:= 2.5358716 pot->attr.mazK[1][2]:= 0.0 pot->attr.mazr2[1][2]:= 999.0 pot->attr.mazd[1][2]:= 0.0 pot->attr.mazD[2][2]:= 5.8563391 pot->attr.mazalpha[2][2]:= 1.9528018 pot->attr.mazr1[2][2]:= 1.8773916 pot->attr.mazK[2][2]:= 1.800 pot->attr.mazr2[2][2]:= 4.150 pot->attr.mazd[2][2]:= 0.330
The old, default way is to use a ZBL-type potential [ZBL85],
1 Z_1 Z_2 V(r) = -------------- ------- phi(r), 4 pi epsilon_0 r
-3.2x -0.9423x -0.4029 -0.2016 x phi(x) = 0.1818 e + 0.5099 e + 0.2802 e + 0.02817 ewhere x = r / a_u and
0.8854 a_0 a_u = ------------------- Z_1^0.23 + Z_2^0.23(a_0 is the Bohr atomic radius = 0.529 Å).
The ZBL repulsive potential mode is selected with pot->rep.type[0][1]:= 1 for the atom type pair 0-1. Alternatively, pot->rep.type[0][1]:= 3 does the same thing, but uses an array interpolation in an attempt to speed things up somewhat. Note, however, that type=3 can lead to small interpolation errors.
In mdsetup the exponential terms of phi(x) above can be modified in the param.in file using the pot->rep options. For instance, for He-Ta we have used the parameters
# # DMol parameters for He - Ta # pot->rep.ai[0]:= 0.137986 pot->rep.ai[1]:= 0.575647 pot->rep.ai[2]:= 0.158109 pot->rep.ai[3]:= 0.128257 pot->rep.bi[0]:= -5.13475 pot->rep.bi[1]:= -1.0958 pot->rep.bi[2]:= -0.404614 pot->rep.bi[3]:= -0.404567The new way to treat the repulsive potential is to read it in as a function of r from a file, and by evaluating the potential using spline interpolation. If the target material is a compound, you can use the built-in potential for some atoms and a precalculated one for others.
The file should be an ascii file with data in the format "r V dV", in a file called reppot.type1.type2.in (typically reppot.0.1.in). Note that the derivative should also be given; this is to force the user to do some thinking on how he should derivate the strongly decreasing potential function. Also, from some DFT/LDA codes like DMol the derivative can be directly obtained from the ab initio calculation, which should give a better derivative than any approximative calculation.
The units in the repulsive potential file should be eV, Å (Angstrom) and eV/Å.
The spline repulsive potential mode is selected with pot->rep.type[0][1]:= 2 for the atom type pair 0-1.
I recommend the newer approach, if it is possible to calculate enough points of the repulsive potential to make the spline interpolation reliable. Otherwise, the best approach probably is to fit a ZBL-type potential reliably to the points one has.
Calculating ion ranges using the RIA approximation is a fairly specialized task, which doesn't resemble typical MD simulations very much. Therefore we felt that writing an entirely new program for this task was justified.
The program was written strictly in ANSI C, which allowed for very easy portability between different machine architectures. The sole exception to the use of ANSI C is made in the subroutine eltime(), which uses the Unix system call times() to obtain the CPU time used by the program. The user interface was originally written in C++ and does still contain some of its features - if this is a problem, contact us (it will be rewritten in ANSI C when we have the time).
Version 1.1a of the program has been successfully compiled and run in all Unix architectures (DEC Ultrix, DEC AXP, HP 9000, Sun OS4, Sun OS5, SGI MIPS, IBM RS/6000, Convex and Linux) where compilation has been attempted. Compilation in non-Unix environments has not been tried, but I believe it should be possible to do without major difficulties. Should you want to use the user interface in a non-unix environment, remove the two lines in mdsetup.C that calls chmod (marked /*Remove if using a non.unix environment*/).
Use of the C language had many advantages compared to FORTRAN-77. The program variables were assembled in logical entities by using data structures. This made the program much more readable than the typical FORTRAN program, and eliminated the need to use global variables. Also, the C preprocessor could be used to handle debugging and variable reading in a compact format.
Memory allocation subroutines were used to allocate space for large arrays during runtime in order to yield flexibility in memory usage. The data type real was used instead of float throughout the program to enable easy transformation of single precision to double.
For instance, time variables were collected in the structure struct time. In the beginning of main(), the time variables are defined with the command struct time time. Thus all time variables can be passed on to a subroutine by simply passing time as a subroutine argument.
Examples of other structures are struct unit for program units, struct atom for atom properties, struct reccalc for recoil calculation variables and struct file for input and output files used in the program.
Internal quantity Multiply with to obtain SI ----------------- -------------------------- E 1.6022e-19 J r 1e-10 m F 1.6022e-9 N m 1.6606e-27*m(u) kg t 10.1805e-15*sqrt(m(u)) s v 9822.66/sqrt(m(u)) m/s a 9.6485e17/m(u) m/s^2Here m(u) denotes the atom mass in atomic mass units.
Figure 1. A schematic view of the MDRANGE program. Subroutine calls are indicated with lines connecting the boxes; the starting point of the call is always higher than the subroutine it calls. The calls are shown in the order (left-to-right) in which they are made in the program.
All files, and the most important subroutines and subroutine calls are shown schematically in figure 1.
Each program file contains one or more C subroutines. Them main program main() is in the file mdh.c. It calls all initialization subroutines, and contains the main simulation loops where the recoil calculation is performed.
The initial displacement calculation is performed in another subroutine inistatecalc() in the file inistatecalc.c. Both main() and inistatecalc() contain somewhat differing versions of molecular dynamics loops over time steps (the central part of most MD programs).
Both loops over time steps work essentially by calling just three subroutines. The subroutines solve1() and solve2() are used to solve the equations of motion; two subroutines are needed because of the nature of the Smith-Harrison algorithm used in the solution. The accelerations() subroutine is used to calculate the accelerations the atoms experience. The accelerations are calculated by derivating a two-body potential obtained from the subroutine potential().
In addition to calling these subroutines, the inistatecalc() subroutine also calls borders() to enforce periodic boundary conditions on the simulation cell. The main program also calls moveslice() to perform the cell atom movings, calctsteps() to calculate the size of the time step and subroutines dealing with the recoil atom properties.
During the loop over recoil events, the simulation results are written to files on regular intervals, by default after every 40 recoil events.
All program parameters have default values, and can be omitted from the parameter file if the default value is what is desired. If the parameter file is empty, the program by default simulates 10 keV Si implantation of crystalline silicon. However, although the parameter file can be left empty, it cannot be omitted.
//SPUTTERING
If the user wants to shift the surface of the material by some value,
he can modify the physical->sput_surf variable. The same variable is used
in the surface erosion by sputtering modeling if physical->sputtering:=
1. The model uses Sigmunds linear sputtering theory, which is based on
the amount of nuclear deposited energy in the surface. The program erodes
the surface during the simulation and takes this effect into account in
the range profile calculations. The input term is the surface binding energy
for the atoms in the surface: layers->layer[0].surfbind.
The sputtering model is somewhat simple, so it might give wrong values
for yields. Better way of examing the effect of erosion is to take the
value of the yield from the experiments and feed it to the program as:
physical->Y and using physical->sputtering:= 2.
!Note, the sputtering calculations are under heavy development, so
it might be somewhat buggy.
//PENR MODEL FOR ELECTRONIC STOPPING
The newest version (V3.0) supports Puska, Echenique, Nieminen, Ritchie
(PENR)-local electronic stopping model, which uses the local electronic
density in the target material and phaseshifts for implanted ion for the
calculations. This calculation is very advanced and good at describing
the electronic stopping for ions (V<V(Fermi)) in channeling directions.
The input files are: phaseshift-based electronic stopping vs. one electron
radius stopping.in, and the spherically symmetric electronic density
as a function of radius for the lattice atoms cdens.in. The type
of the electronic stopping model is defined in param.in file: layers->layer[0].elstop:=
1 (PENR) 0 (ZBL=DEFAULT).
! Note, if the PENR model is chosen, the electronic
density at the point of the ion is calculated by summing densities of surrounding
atoms inside a sphere of a radius defined in param.in as potcrit->rfirsov.
It has default value of 1.47 Å for Si. It is the same radius that
is used for Firsov model for inelastic el. stopping, which is set on as:
elstop->Firsov:= 1. If the PENR model is used, the Firsov model is not
necessary. This radius should be set with some physical thinging involved
if changing the target material.
//BK MODEL FOR ELECTRONIC STOPPING
Easier way to calculate local electronic stopping is using Brandt-Kitakawa
(BK)-theory and correction terms. For this user has to give charge density
file for the lattice atom type,
as above, but the program itself precalculates the stopping vs. one
electron radius function by using linear theory and a correction function.
The type of the electronic stopping model is defined in param.in file:
layers->layer[0].elstop:= 2.
reccalc->Ncalc:= 10000 layers->Natomsmax:= 100 reccalc->Theta0:= 6.0The simulation thus simulated implantation of 10 keV Si in c-Si, with the implantation direction deflected 6° from the <100> crystal axis. The ZBL interatomic potential and an experimental electronic stopping measured at our laboratory were used in the simulations.
The resulting range and deposited energy distributions are shown in figure 2. All the program input and output files can be viewed by clicking on the file names in the list above. The standard output of the run is available here.
Figure 2. The range and deposited energy distributions of 10 keV
Si implantation of crystalline Si, beam aligned 6° from the <100>
crystal direction. The surface peak depends on the interpretation of deposited
energy at the surface.
# # Input parameters for mdh # # # # Atom types # type[0].Z:= 7 # Projectile type[0].m:= 15.00 # N-15 type[1].Z:= 14 # Target type[1].m:= 28.08 # Si # # Recoil calculation parameters # reccalc->E0:= 100000 reccalc->Atype:= 0 # Recoil atom type reccalc->Ncalc:= 10000 # Number of histories reccalc->Trange:= 0 # Range type: 0=sz, 1=proj, 2=chord, 3=path reccalc->type:=1 # Recoil calculation type: 1=normal, 2=negative ranges allowed reccalc->Fii0:= 0.0 # Azimuth angle of the recoil atom reccalc->Fiimax:= 0.0 # is chosen randomly from [Fii0,Fiimax] reccalc->Theta0:= 8.0 # Polar angle of the recoil atom reccalc->Thetamax:= 8. # is chosen randomly from [Theta0,Thetamax] # All angles are measured in degrees. reccalc->Startmin.x:= 2.715 # Starting position x coordinate chosen randomly between reccalc->Startmax.x:= 8.145 # [Startmin.x,Startmax.x]. reccalc->Startmin.y:= 2.715 # Likewise y and z. Å. Unless you want to reccalc->Startmax.y:= 8.145 # emphasize any region of the unit cell, you reccalc->Startmin.z:= -2.715 # should set Startmax.x - Startmin.y = l, reccalc->Startmax.z:= -2.715 # Likewise y. Z should be negative to # ensure that the recoil atom does not start # from a position inside the grid. # # Repulsive potential # pot->rep.type[0][1]:= 1 # Exponential screening used for rep.pot. pot->rep.type[0][1]:= 2 # Potential read from file (reppot.0.1.in) pot->rep.type[0][1]:= 3 # Exponential screening with array speedup # # Attractive potential - These settings can be omitted if thermal # displacements are calculated using the Debye temperature # pot->attr.type:= 1 # Mazzone potential for silicon (and # other diamond structures pot->attr.type:= 2 # Morse potential # # Layer parameters # layers->Nlayers:=1 # Number of layers layers->Natomsmax:= 100 # Maximum number of atoms in a layer layers->Ninistate:= 10 # Number of inistate calculations layers->layer[0].minz:=-1.0e9 # Minimum z of layer 0 layers->layer[0].maxz:=1.0e9 # Maximum z of layer 0 layers->layer[0].ltype:=1 # Layer type (i.e. coordinates read in from coords.in.ltype). Must be > 0 layers->layer[0].Timeinior:=8.0 # Time of inistate calculation # # Misc params # box->Max.x:= 10.8 # Size of MD cell, Å box->Max.y:= 10.8 # box->Max.z:= 10.8 # box->movelim.x:= 2.715 # Shift used in generating fresh crystal box->movelim.y:= 2.715 # in front of the recoil. box->movelim.z:= 2.715 # 0.5*lattice const. is a good value, Å. physical->Tini:= 300.0 # Temperature of the inistate calculation physical->debye:= 1 # Debye temp. used in inistate calculation physical->tdebye:= 375 # The Debye temperature of lattice atoms, K time->Timeini:=1000.0 # Time of the inistate calculation gen->seed:= 897234121 # Random number generator seed potcrit->Rmini:= 6.0 # Neighbour list construction cut-off # radius in rec.calc potcrit->R0ini:= 5.0 # Potential cut-off radius in ini.calc. potcrit->Rmrec:= 3.3 # Neighbour list construction cut-off # radius in rec.calc potcrit->R0rec:= 2.7 # Potential cut-off radius in rec.calc. # (Must be smaller than movelimit!) # All cut-offs in Å. elstop->scale:= 1.12 # Electronic stopping power scaled by this pot->rep.scale:= 1.0 # Repulsive potential scaled by this
An example for data in this format:
dpc nticks 19 19 coldcolors ntdown sd 890 890 d 1 contcolor col 0 3000 colcolsoft 3 ticks xy w x 0 90 y 0 90 sort 2 1 3 4 results_channeling_range_interpolated.datAlso plotting of theta, phi in polar coordinates is now possible. An example:
dpc nticks 19 19 coldcolors ntdown sd 890 890 d 1 contcolor col 0 3000 colcolsoft 3 ticks xy w x 0 90 y 0 90 sort 2 1 3 4 results_channeling_range_interpolated.dat
| $ | The start of a new record |
| CaF | The name of the compound, max 70 chars |
| DIA | Grid type. Supported type are DIA, HEX, BCC, FCC and MUU for all others. ("muu" if Finnish for "other" - thanks/blame to Jussi Sillanpää on this feature...) |
| 250 | Debye temperature, zero if N/A |
| 2 | Number of elements in the compound |
| 4.22 | grid length, followed by height (HEX, MUU), width (MUU) and number of atoms in the simulation cell |
| 1
1 etc. |
Atom types |
| 29
etc. |
Atomic number of atom type 1...n |
| 58.61
etc. |
Mass of atom type 1...n, followed by the coordinates of the atoms in the simulation cell (MUU;x1,y1,z1,...) |
| # | End of the record |
Should you add any new compounds, please post me your compound file so that it can be distributed to other users.
which give the corners of a box within which the position of the recoil atom is selected randomly. There is now weighting on the positions, all positions are selected with equal probability. Usually it is advisable to have the inital position somewhat outside the cell, for instance reccalc->Startmin.z:= -3.0, reccalc->Startmax.z:= -3.0 . The x and y box limits are usually well advised to be set to the size of one unit cell in the crystal.
To find out what the default values of these are, look into the source code of readin.c, and at the standard output file startdata.out.
The initial direction selection process is somewhat complex. If the parameters reccalc->Fiimax and reccalc->Thetamax are not given, or are equal to zero, the inital direction fii and theta are set to reccalc->Fii0 and reccalc->Theta0. Theta is the angle from the normal of the simulation cell surface (i.e. from the z axis), and fii is the rotation angle. Both are input in units of degrees, but are in radians internally in the program.
If reccalc->Fiimax and reccalc->Thetamax are given and are nonzero, fii and theta are selected randomly. Fii is selected between Fii0 and Fiimax, Theta between -Thetamax and Theta0 or between Theta0 and Thetamax, using the correct weighting. To find out exactly how this is done, look into the source code of recoil.c
[Maz91] A. M. Mazzone, Interatomic Potentials in Silicon, Phys. Stat. Sol. (b) 165, 395 (1991)
[ZBL] J. F. Ziegler, J. P. Biersack, and U. Littmark, The Stopping and Range of Ions in Matter, 1985,
[Sil00] J. Sillanpää, K. Nordlund, and J. Keinonen, Electronic stopping of Silicon from a 3D Charge Distribution, Phys. Rev. B 62, 3109 (2000).
[Cai96] D. Cai, N. Grønbech-Jensen, C. M. Snell, and K. M. Beardmore, Phenomenological electronic stopping-power model for molecular dynamics and Monte Carlo simulation of ion implantation into silicon, Phys. Rev. B 54, 17147 (1996)
[B98] K. M. Beardmore and N. Gronbech-Jensen, An Efficient Molecular Dynamics Scheme for the Calculation of Dopant Profiles due to Ion Implantation, Phys. Rev. E 57, 7278 (1998).
[Sil99] J. Sillanpää, K. Nordlund, and J. Keinonen, Channeling in manufacturing sharp junctions: a molecular dynamics simulation study, Physica Scripta T79, 272 (1999).
[5] K. Nordlund, Molecular dynamics simulation of ion ranges in the 1 --
100 keV energy range, Comput. Mater. Sci. 3, 448 (1995).
[9] P. Haussalo, K. Nordlund, and J. Keinonen, Stopping of 5 -- 100 keV
helium in tantalum, niobium, tungsten, and AISI 316L steel, Nucl. Instr.
Meth. Phys. Res. B 111, 1 (1996).
[58] J. Sillanpää, J. Peltola, K. Nordlund, J. Keinonen, and M. J. Puska,
Electronic stopping calculated using explicit phase shift factors, Phys.
Rev. B 63, 134113 (2000).
[78] J. Peltola, K. Nordlund, and J. Keinonen, Effects of damage build-up
in range profiles in crystalline Si; molecular dynamics simulations, Nucl.
Instr. Meth. Phys. Res. B (2002), accepted for publication.
[92] J. Peltola, K. Nordlund, and J. Keinonen, Molecular dynamics
simulation method for calculating fluence-dependent range profiles, Nucl.
Instr. Meth. Phys. Res. B (2002), COSIRES 2002 conference paper, submitted
for publication.