The MDRANGE (mdh) program

Authors of the program V1.0 - V1.83b

Kai Nordlund

Accelerator Laboratory, University of Helsinki, P.O. BOX 43, FIN-00014 Helsinki, Finland

Author of the user interface & some minor modifications

Jussi Sillanpää,

September 1, 1994

Author of the program development from V1.X -> V3.0

Jarkko Peltola,

Last update: Jarkko Peltola and Kai Nordlund, August 2002.



This document contains a brief description of the MDRANGE program and its setup program, mdsetup. This document evolves slowly as I write in more stuff every now and then. It was originally written to describe MDRANGE V1.01a, but now also covers some newer stuff.

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.

Summary of physical principles

This is a brief summary of the physical principles used and features incorporated in version (1.52b) of the program. For more info, see Molecular dynamics simulations of ion ranges in the 1-100 keV energy range (Comp. Mat. Sci. 3 (1995) 448-), as well as the other references cited at the end of this document. .

What it does


Supported sample structures

Treatment of interactions

Tests of simulation accuracy

The accuracy of the integration accuracy of mdh V1.8 has been tested by comparing a single binary collision with a BCA integration. To very briefly summarize the results, using the default parameters the code treats the scattering of the test cases (1 GeV He -> Si, 5 MeV Cu -> Cu and 5 MeV Au->Au and lower energies) with an accuracy of 0.1 degree or better for all impact parameters tested.

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.

Using mdh / mdhreed / mdsetup

Program compilation

mdh, mdhreed and mdsetup are compiled by uncompressing mdrange.tar.gz and issuing the command make. Compilation requires an ANSI C compiler (by default gcc is used) and of course the make (1) system command.

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.

Running mdh

To run mdh you need a least 3 input files, and often 4:
  • the coordinates of the atoms in the simulation cell
  • the electronic stopping (m/s, eV/Å)
  • the simulation parameters
  • the repulsive potential (Å, eV, eV/Å) [Optional].. If not provided, mdh can generate the ZBL repulsive potential automatically. All input (and output) are just plain ascii files. To get these you have basically 4 ways to go:

  • Generate all of them yourself from scratch
  • Start with some of the examples provided in the program package and modify them to suit your needs
  • Use the mdsetup program, see below

    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 by hand.

    Running mdhreed

    mdh now also has a special version mdhreed, a version of mdh which uses a rare-event algorithm similar to the one introduced by Beardmore and Gronbech-Jensen [B98,Sil99]. However, because the reed method holds a U.S. patent, No. 6,366,873, we do not distribute the reed method with the main package, but only to users outside the U.S. on special request.

    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:= 1 
    and 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.

    Running mdsetup

    This package includes mdsetup, a user interface to mdh. You can use mdsetup to create the input files, and 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 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 and 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 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.

    Debug output

    The easiest way to manipulate the output is to make a script using mdsetup. To enable easy debugging two debug options are incorporated in the program. Ordinary debug mode is implemented with the command line option "-d". It prints out to standard output and standard error the values of a number of important variables during run-time. Variables from within the loop over time steps are not printed out for each time step to prevent the number of output data from becoming too large. In maximum debug mode (selected with "-dm") variables from within the loop over time steps along with a large number of other information is printed out.

    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 :-(

    Regular output

    In the beginning of the simulation data about program parameters is written to the file startdata.out. The data includes values for all parameters the program attempts to read in from, including a note on whether the parameter was found in or given its default value.

    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 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/Å.

    Incident angle selection

    The incident angle/crystal orientation is selected using the parameters theta and phi.The crystal is always oriented the same way, but the 'beam' direction is changed by changing the initial direction of each incident ion.

    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.
  • theta determines the polar (tilt) angle, i.e. how much the incident direction is deflected from the z axis.
  • phi determines the azimuthal (twist) angle, i.e. how much the incident direction is deflected from the x axis
  • The exact algorithm used to determine the velocity components (vx,vy,vz) for a given velocity v, theta and phi is:

    Thus e.g. Fii0=Fiimax=0 and Theta0=Thetamax=0 chooses implantation straight into the 001 channel.
    For <110>-implantations in FCC-materials, user can use normal 100-oriented coordinate file (which mdsetup creates) and set Fii0=Fiimax=0 and Theta0=Thetamax=45.
    For <111>-implantations in FCC-materials, user can use normal 100-oriented coordinate file (which mdsetup creates) and set Fii0=Fiimax=45 and Theta0=Thetamax=54.
    ! Note that if the <110> or <111> profiles calculated as above, one has to use  projection obtion (Trange=1) to have the profiles in those directions.

    Creating random alloys

    To create a random alloy, you can use the type[i].prob parameter. For example, to create a random CuAg alloy, first create a Cu FCC cell as usual, but with a lattice constant corresponding to CuAg. Then set the atom types as:
    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.5
    Of course you have to remember to set appropriate reppots, electronic stopping and Debye temperature.

    Notes on using attractive potential

    During the initial state calculation, the program uses either a 2-body potential to model equilibrium interactions to give atoms realistic initial thermal displacements or calculates the displacements using the Debye temperature of the substance. For Si (and other systems with the tetrahedral diamond structure) the program uses a modified Morse potential potential by Mazzone [Maz91]. For metals, the ordinary Morse potential is used.

    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 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

    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

    Notes on using repulsive potential

    From MDRANGE V1.6 forward, the repulsive potential can be handled in two entirely different manners.

    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 e
    where 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 file using the pot->rep options. For instance, for He-Ta we have used the parameters

    #  DMol parameters for He - Ta
    pot->[0]:= 0.137986
    pot->[1]:= 0.575647
    pot->[2]:= 0.158109
    pot->[3]:= 0.128257
    pot->[0]:= -5.13475
    pot->[1]:= -1.0958
    pot->[2]:= -0.404614
    pot->[3]:= -0.404567
    The 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 (typically 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.

    Programming principles

    We have previously calculated ion ranges using an old molecular dynamics code ("MOLDY") originally written to simulate properties of vacancies in semiconductors. Because this program was originally written in the 1970's in FORTRAN-66, and had since been used for many different purposes, it carried with it lots of awkward and redundant features. In particular, it used an integer time step, which made implementation of a variable time step difficult.

    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.

    Use of data structures

    As mentioned above, data structures were used throughout the program. They structures are all defined in the header file structs.h. The actual variables using the structures are defined in the beginning of the main program main().

    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 unit system

    To reduce the amount of arithmetic operations needed during the simulation, the program uses an internal unit system where all masses equal to 1. The internal energy unit is eV and the length unit Ångström. In this system all internal units involving time will depend on the atom mass. These units include the time, velocity and acceleration units. Some important conversion factors to obtain quantities in SI units from the internal units are:
    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^2
    Here m(u) denotes the atom mass in atomic mass units.

    Program organization

    The program consists currently (V1.01a) of 17 ".c" program files, five ".h" header files and a makefile for the unix command make (1). The main header file common.h loads in all the other MDRANGE header files and a number of necessary ANSI C system library headers. common.h is then included in the beginning of all MDRANGE program files, except readvar.c. This subroutine (which is used to read in program parameters from a file) is also used in the FORTRAN program MOLDY. To retain maximum portability of the subroutine, it does not contain any properties that would be specific to either program. Hence it does not use the MDRANGE include files.

    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.


    MDRANGE V1.01a program files

    Note: many of these have bugs which have since been corrected

    More about the input files 

    The parameters of the simulation run are contained in a file with the default name In the parameters file the parameters are marked with the variable name in the program, ":=" after the name and the value after the equal sign. Between ":" and "=" comments can be placed. The parameters can be in arbitrary order in the parameters file. A list of the parameters is in the appendices.

    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.

    coordinates file

    The atom coordinates are read in from the file The three first numbers on each line give the initial atom coordinates and the fourth number gives the atom type. The atom types are defined in the parameters file. If multiple layers are used in the simulations, coordinates for each layer are read in from the files,, etc.
        MDRANGE3.0: In the newest version (V3.0) the coordinate files can have some variation in size. For example, the simulation box does not have to be the same size as the coordinate files for each layer. The only restriction is that, the simulation box (defined in the file) has to be smaller or equal size as the layer coordinate boxes. The size of layer structure in the coordinate file is defined in the first line as: V a b c na nb nc, where 'V' marks the line, 'a,b,c' are the unitcellsizes for the layer and 'na,nb,nc' are the number of unitcells in each direction.  The newest version of MDSetup does this automatically.
    During the reading of the coordinates, the program moves them to a box of (0,0,0,a*na,b*nb,c*nc) so user don't have to think what the offsets for the coordinates in each layer are, if only the offsets are not multiple unitcells away from the box that they are moved to!

    electronic stopping

    The electronic stopping is read in from the file The first column in this file is interpreted as the ion velocity in m/s, and the second as the electronic stopping in units of eV/Å. By using the option elstop->mode:= 1 in the parameter file, the straggling of electronic stopping can be included in the calculations (in V1.5a or bigger).
        MDRANGE3.0 supports multiple layer structures, so that the elstop files are named, etc. if there is more layers than one.

    Advanced options for MDRANGE 3.0

        There is also a possibility to describe the amorphization of the target material for high fluence implantations. This is done by making multiple coordinate files for each layer where the degree of damage changes as a function of dose (eV/atom) in that layer. The coordinate files are named as, etc. where the first number is the number of the layer and the second is the number of amorphization level. The maximum number of the amorphization levels for each layer is defined in file as: layers->layer[1].Namorph:= 2 (in this case) and setting physical->Damage:= 1. Because the volume of the layer changes during the amorphization, there is also an option in file: physical->scaling:= 1. This takes care of the volume increases during the simulation. The method is still under development.

    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.

    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, and the spherically symmetric electronic density as a function of radius for the lattice atoms The type of the electronic stopping model is defined 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 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.

    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 file: layers->layer[0].elstop:= 2.

    List of header and other files

    Default input and output file names

    Clicking on the file names show you sample parameter files obtained from the test simulation described below

    Sample run

    A sample run of MDRANGE was performed using the following parameters file:
    reccalc->Ncalc:= 10000
    layers->Natomsmax:= 100
    reccalc->Theta0:= 6.0
    The 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.

    List of the parameters

    (for a complete list of parameters, run mdh and read the startdata.out file)
    # 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 (
    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 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

    Some examples

  • mdh can also be used to make color scale 2D data plots, simply by setting the dot size to 1 (or 2) and using a regular grid of points at the same time (in the simple x y h 0 format).

    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.dat
    Also 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

    More about the compound library of mdsetup

    The compound data is stored in compounds.txt . The structure of the file is following:

    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 
    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 
    Atom types 
    Atomic number of atom type 1...n 
    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.

    Notes on incidence of the atom

    For experienced mdh users only: There are a number of ways in which the atom can enter a crystal. The initial position of the atom is determined by six parameters:
    reccalc->Startmin.x, reccalc->Startmax.x,
    reccalc->Startmin.y, reccalc->Startmax.y,
    reccalc->Startmin.z, reccalc->Startmax.z,

    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

    Back to beginning of document

    References cited

    [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).

  • The original reed paper

    [Sil99] J. Sillanpää, K. Nordlund, and J. Keinonen, Channeling in manufacturing sharp junctions: a molecular dynamics simulation study, Physica Scripta T79, 272 (1999).

  • Our first implementation of a rare-event algorithm

    Other relevant references

    [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.