Home

FHI98MD Computer code for density-functional theory calculations

1. e To avoid folding of the k point set all components of the i_facs parameter in start inp should be set to 1 e Remember to take enough of empty states parameter nempty in start inp if you want to get the lowest unoccupied bands as well Parameter Value i_facs 1 3 i i il k points folding factors no folding t_kpoint_rel false k points are given in Cartesian coordi nates in units 27 Qjat true frame of reference for k points is spanned by the reciprocal lattice vectors i e k ky by kobe kg b3 e Projected bulk band structure Calculate the bulk band structure for the same set of kj points as for the surface for which you want to project the bulk bands but with k going from 2r a_ to 2r a1 where a is the unit cell size in the perpendicular to the surface direction If the symmetry of your system includes a mirror plane parallel to the surface then it is sufficient to take only half of the interval e g kq 0 27 a_ e Other parameters epsekinc inp mod use this parameter to control the eigenvalues convergence d Analysis e The k point coordinates together with the corresponding sorted Kohn Sham eigenvalues in eV and occupation numbers are listed at the end of the output file fort 6 see Example section For example the following command line on UNIX systems would extract the eigenvalues to the file eigenvalues dat grep gt gt eig fort 6 sed s gt eig g gt eigenvalues dat
2. 2 2 The input files inp mod and start inp 13 Parameter ekt Parameter ibrav Parameter ion_damp Parameter ion_fac Type Value Type real integer oOo WwW N Type real temperature of the artificial Fermi smear ing of the electrons in eV if tmetal true this parameter has to be opti mized for each system in order to ob tain optimal convergence For instance in case of semiconductor systems smallness of this parameter would ensure that the semiconducting character is unaffected Commonly used values lie in the range 0 01 0 10 eV connected to tmetal cell type the cell types are specified in the routine latgen and can be specified individually the most common cell types are available structure is kept as supplied in the file start inp simple cubic lattice fec lattice bec lattice orthorhombic connected to celldm see also How to set up atomic geometries damping parameter 0 1 only active if tfor is set to true and tsdp is set to false connected to tfor and tsdp see also How to set up a structural relaxation run Type real mass parameter if tfor is set to true if tdyn is set to true ion_fac specifies the ionic mass in a u connected to tdyn and tfor 2 2 The input files inp mod and start inp 14 Parameter Type i_facs 1 3 integer k point folding factors 1 1 1 corre sponds to no folding connected to tfor and tsdp see also Choice of the k point mesh
3. Parameter Type lmaz integer highest angular momentum of the pseu dopotential 1 s 2 gt p 3 gt d connected to lloc Parameter Type Lloc integer angular momentum of the local pseudo potential lloc lt lmaz connected to Imax Parameter Type na integer number of atoms of speciesis has to be specified for each element Parameter Type nel_exc real number of excess electrons required for a calculation with a charged supercell for bulk defect calculations Parameter Type nempty integer number of empty states Parameter Type only for MD simulations if tdyn is true nfi_rescale integer number of ionic moves before velocities are rescaled only active if nthm is 1 connected to nthm Parameter Type nkpt integer number of k points 2 2 The input files inp mod and start inp 15 Parameter Type only relevant for parallel code npes integer number of processor elements PE s minpes integer minimal number of PE s currently inactive ngrpx integer number of PE s in a group performing the FFT currently inactive Parameter Type Value only for MD simulations if tdyn is true npos integer set up of initial ionic positions and veloc ities only active if tdyn is set to true 1 coordinates taul and velocities vau0 are read in from file start inp 2 coordinates are read in from file fort 70 velocities from file fort 20 3 coordinates are read in from file fort 70 the velocities are set according to the ini ti
4. Under certain circumstances it might be desirable to continue a previous calculation because e the run time in the queuing system of your hardware was too short to fully converge the calculation e the calculation terminated prematurely due to some computer problem e the plane wave cut off or the number of k points needs to be increased to ascertain convergence with respect to these parameters Restart files in binary format are written every iprint self consistency cycles c f inp mod The file fort 71 contains the wave functions and the coordinates of the ions as well as all other information required for a restart The electron density is stored in the file fort 72 11 General continuation run or run with increased cut off The following parameters need to be reset for performing continuation runs Notice that some parameters may be set in start inp or likewise in inp ini In the latter case running the fhi98start program can be omitted e Rename the file fort 71 into fort 70 e In the file inp mod set the parameter nbeg to 2 e In the file inp ini or start inp set nrho 2 which causes the program to constructs the initial electron density from fort 70 When nrho 3 it reads the initial electron density from fort 72 instead e Usually one would also set coordwave true which causes the pro gram to read the last available ionic positions from fort 70 Otherwise the program reads the ionic positions from inp ini Th
5. 000048 a u gaussian energy 5 801080 a u amp amp s atomic positions and local nl forces on ions Galium Arsenic time elapsed for vofrho t 0600 time elapsed for n x nkpt x graham ortho 2200 nel dfmax fscale efermi ekt seq sneq 8 00000 030 40000 2 60750 005 00000 00000 gt gt gt n_it nfi Ekinc Etot Eharr Ezero mForce mChange Seq Sneq Efermi Dvolt Wa gt gt gt 1 0 2 70633 8 30572 8 45169 8 30572 00000 000 0000 0000 2 6075 0000 0000 I finished storing wavefunctons and data on fort 71 gt gt gt OK I stop after timestep nr 10000 time elapsed per electronic time step t 1 9200 time in queue 14400 max number of steps 7111 gt gt gt 2 0 2 83522 8 30155 8 44343 8 30155 00000 000 0000 0000 2 6057 0000 0000 gt gt gt 22 0 09358 8 29527 8 44641 8 29527 00000 000 0000 0000 2 6057 0000 0000 LOOP n_it 23 phfac is n_ideal 1 0 phfac is n_ideal 2 0 RHOE read from unit7 in RHOOFR Number of electrons in real space 7 99999630944237516 internal energy at zero temperature 8 295270 a u non equillibrium entropy 000000 kB equillibrium entropy 000000 kB kT energy 005 eV non eq free energy 8 295270 a u non eq total energy 8 295270 a u Harris energy 8 446409 a u kinetic energy 3 225369 a u electrostatic energy 10 000821 a u real hartree energy 763345 a u pseudopotential energy 1 472645 a u n l pseudopotential energy 62243
6. 4019882 2 2 8405414 3 2 8405414 4 2 8405414 oniskb is 2 wnl 1 3798314 2 1 6212706 3 1 6212706 4 1 6212706 2 6705 2 6705 2 6705 formfa rho of atom in 3 000000 formfa rho of atom in 5 000000 ELECTRON DENSITY WILL NOT BE UPDATED SYMMETRY OPERATIONS gt s isym in latt coord 1 0 0 0 1 0 gt sym in cart coord 1 00000 00000 00000 00000 1 00000 00000 gt sym in cart coord 1 00000 00000 00000 00000 1 00000 00000 gt Center of symmetry sym0 000000 000000 000000 Table of symmetry relations of atoms is ia iasym isym taud xneu 1 1 1 1 00000 00000 00000 00000 00000 00000 2 67050 2 67050 2 67050 2 67050 2 67050 2 67050 14555 initialization loop ik 1 Number of plane waves 40 Number of localized orbitals 8 Total number of states for initialization 48 time consumption for a single k point time used for mapping of G G 00 time used for creation of BO 35 time used for lt Gt k H_n1 G k gt 03 time used for lt Gtk H_loc G k gt 00 time used for lt BG H BG gt 04 time used for diagonalization 02 ik 40 Number of plane waves 51 Number of localized orbitals 8 Total number of states for initialization 59 nel dfmax fscale efermi ekt seq sneq 8 00000 030 1 00000 2 60750 005 00000 00000 k point 500 500 500 eigenvalues and occupation numbers eig 8 584 4 211 991 991 3 318 7 045 7 045 10 305 12
7. Alternatively the k points coordinates kj weights wz and the cor responding eigenvalues together with some additional information are saved in the file report txt in the format report txt number of k points number of eigenvalues per k point N ki kiz kis Wk Eil EiN VIZ COD COCORRFRRFRRFRRFRFROODOOOOCOOCCOCCOCOCOCOCOCOCCSO 682 0 0 75 4422650 3845300 3267950 2690600 2113250 1535900 0958548 0381198 oono BP Une 9292890 8585790 7878680 75 6792890 6085790 5378680 3 6 How to set up a band structure calculation 51 The parallel version of fhi98md writes also an output file fort 3 containing the eigenvalues e The value of the Fermi energy in eV respectively the Fermi level location can be traced during the calculation and picked from the Efermi field of the information record written to fort 6 at each electronic iteration gt gt gt n_it nfi Etot Eharr Ezero mForce mChange Seq Sneq Efermi gt gt gt 23 0 8 29527 8 44641 8 29527 00000 000 0000 0000 2 6057 It is important to remember that the Fermi energy should be taken from a real self consistent calculation e g the one used to generate the fort 72 file and not from the band structure run see also Fig 3 8 1 Example e The sample input file start inp below sets up a bulk band structure calcu lation for GaAs along the standard symmetry lines of the Brillouin zone inset in Fig 3 8 of the face c
8. Wiliams Soler algorithm first order damped Joannopoulos algorithm second order exchange correlation functional LDA Ceperley Alder 18 Perdew Zunger 19 Becke Perdew BP 20 21 Perdew Wang PW91 22 Becke Lee Yang Par BLYP 23 Perdew Burke Ernzerhof PBE 24 only relevant for MD simulations scheme for solving the ionic equation of motion only active if tdyn is true predictor corrector predictor Verlet algorithm connected to nOrder and tdyn type of basis set used in the initialization if nbeg is set to 1 il plane wave basis set 2 LCAO basis set 3 mixed basis set LCAO and plane waves connected to ecuti and tinit_basis in file start inp Type integer number of electronic iterations between a detailed output of i energies and eigenvalues in file fort 6 and i restart files fort 70 wave functions and iii fort 72 electron density 2 2 The input files inp mod and start inp Parameter max_no_force Parameter Type integer maximum number of electronic iterations for which no local forces shall be calcu lated per ionic step Type mesh_accuracy real Parameter nbeg Parameter nomore nomore_imt Value degree to which the sampling theorem shall be satisfied The sampling theorem implies that the size of the Fourier mesh nrx 1 xnrx 2 xnrx 3 which determines the accuracy of the charge den sity should obey nra 1 gt ai WEcut likewise f
9. Choice of the k point mesh 34 Parameter Value nkpt 1 number of k points supplied wk 1 8 wkpt 0 5 0 5 0 5 1 0 k points and weights i_facs 1 3 444 k point folding factors frame of reference for k t_kpoint_rel true points zk 1 3 e k point set for a slab calculation For a surface calculation with the z axis as the surface normal you want the k point mesh to lie in the xy plane There is no dispersion of the electronic band structure of the slab in z direction to sample If there would be it just means that the repeated slabs are not decoupled as they should be i e the vacuum region was chosen too thin Therefore the z coordinate of all k points should be zero The input typically looks like Parameter Value nkpt 1 number of k points supplied wk 1 8 0 50 5 0 0 1 0 k points and weights ifacs 1 8 881 k point folding factors t_kpoint_rel true frame of reference for k points xk 1 3 Figure 3 1 2D Brillouin zone of a surface with cubic symmetry with a 8 x 8 Monkhorst Pack grid The thin square indicates the conventional first Brillouin zone the thick square marks the Brillouin zone as realized in the fhi98md code The location of one special k point out of 64 within its tile is marked by the cross VIZ 3 2 Choice of the k point mesh 35 Note We recommend to use even numbers for the folding parameters As a general rule one should avoid using high symmetry points in
10. achieved or when the preset number of iterations or the allowed CPU time is exceeded Output is generated at the last block of the self consistency loop and by the routines fiondyn fionsc fermi and vofrho The routines fiondyn and o_wave generate restart files for MD simulations and total energy calculations In the mixed basis set initialization the self consistency loop closely follows the organization of that discussed above First of all the initial electron density is obtained either from a superposition of contracted atomic pseudo densities or from an electron density of a previous calculation file fort 72 The local contributions to the potential and the energy are calculated by the routine vofrho The localized orbitals to construct the mixed basis set are set up by the routine project The non local contributions to the potential and the energy in the localized basis set are calculated by the routine nlrhkb_b In the second step the Hamiltonian is constructed with the help of routine dforce_b The Hamiltonian is diagonalized by standard diagonalization routines The new electron density is calculated routine rho_psi Finally the new electron density is mixed with the old density by a Broyden mixing routine broyd 2 2 The input files inp mod and start inp 6 2 2 The input files inp mod and start inp Two input files are required as input for the start utility fhi98start The file inp mod contains the control parameters for the run
11. instead in order to ensure the stability of convergence Note that the values for delta2 and gamma2 should be smaller than that ones for delta and gamma 3 3 Total Energy Minimization Schemes 40 Parameter Type Value delt real gt 1 Step length of electronic iteration gamma 0 lt real lt 1 Ift_edyn 2 see below damping param eter delt2 see delt c f eps chg dlt gammaz2 see gamma c f eps chg dlt eps_chg_dlt real gt 0 If the total energy varies less than eps chg dlt delt2 and gammaz2 replace delt and gamma e Choose the desired scheme to iterate the wave functions in the input file inp mod setting up the parameter i edyn according to the following numbers shown below Parameter Value i_edyn Schemes to iterate the wave functions 0 steepest descent 1 Williams Soler 2 damped Joannopoulos 1 Example e As an example we consider a total energy calculation for GaAs bulk in the zinc blende structure The Williams Soler minimization scheme is used to iterate the wave functions and different values for the electronic time step delt were used The results are shown below Note how the convergence can be accelerated by changing the electronic time step delt The input files start inp and inp mod used for these calculations are shown below start inp 1 number of processors 1 number of minimal processors 1 number of processors per group 2 number of species 0 excess electrons 4 number of empty states 2 0 ibrav p
12. space agree in both calculations to check whether this is actually the case inspect the list of k points appearing in the inp ini file This goal can be achieved in a simple way by choosing appropriate i_facs As an example imagine you want to compare two slab calculations one with a 2 x 1 the other with a 4 x 2 unit cell In this case use Parameter Value i_facs 1 3 481 k point folding factors in the first case and Parameter Value i_facs 1 3 241 k point folding factors in the second case leaving the other parameters unchanged Note b is orthogonal to the real lattice vectors ay and ag If a is the long edge of your real space unit cell b spans the short edge of your Brillouin zone Therefore the k point sampling mesh has fewer points in the b direction and more points in the b direction in the above example Chadi Cohen mesh Another convention for choosing a k point mesh has been proposed by Chadi and Cohen 31 and has been applied to slab calculations by Cunningham 32 In contrast to Monkhorst and Pack the refinement of the k point mesh to ob tain higher sampling density is based on a recursive scheme However for cubic symmetry the outcome of this algorithm can also be interpreted as a special Monkhorst Pack grid To discuss differences between the schemes we resort to the simple case of a two dimensional mesh for a slab calcula tion An example where Cunningham s scheme leads to results differen
13. the Brillouin zone as sampling points because this would result in an inferior sampling quality at comparable numerical effort compared to a similar number of off axis k points Conventionally in contrast to our above definition the Brillouin zone is chosen to have the origin in its center For odd numbers of the folding parameters and the setting 0 5 0 5 some of the unfolded k points will fall on the zone boundary of the conventional Brillouin zone which is often a symmetry plane Likewise the k point set may contain a periodic image of the T point This is normally undesirable The concept of equivalent k points Usually one is not interested in the total energies themselves but in comparing different structures i e accurate energy differences are required If the two structures have the same unit cell the comparison should always be done using the same k point set so that possible errors from a non converged k point sampling tend to cancel out A similar strategy can also be applied when comparing structures with different unit cells We allude to this concept here as equivalent k point sampling The structure with a large unit cell has a smaller Brillouin zone associated with it The k points sampling along this smaller Brillouin zone should be chosen as a subset of the k point mesh in the larger Brillouin zone such that the position of the k points in this subset expressed in Cartesian coordinates in reciprocal
14. units depend on ibrav L Setting up parameter values e First the user specifies whether he she wants to set up the crystallo graphic vectors a a2 a3 spanning the unit cell directly ibrav 0 in start inp or prefers to select the type of the super cell to be used in the calculation from a pre defined list ibrav gt 0 In the first case the crystallographic vectors a a2 a3 should be spec ified in three input lines at the end of start inp In the latter case for ibrav gt 0 the fhi98start utility calls the latgen routine latgen f that sets up the crystallographic vectors a a2 a3 In both cases their reciprocal counterparts b1 b2 b3 are calculated The following lattice symmetries are implemented Parameter Value ibrav 0 user supplied cell 1 simple cubic sc 2 face centered cubic fcc 3 body centered cubic bcc 4 hexagonal Zn bulk hcp A3 structures 8 orthorhombic 10 rhombohedral A7 structures 12 base centered orthorhombic 411 A20 structures e Specify how fhi98start should determine point group symmetries parameter pgind in start inp With pgind 0 an automatic search for the point group symmetries and the symmetry center is performed This is the preferred setting No tice however that all input coordinates may get shifted if this enables to enhance the number of symmetries In successive restart runs if there is the risk that the number of point group symmetries variable nrot may change
15. unwantedly between the runs set pgind 1 in order to avoid conflicts pgind values in the range 2 32 could be used for in stance for some specific purposes in structural relaxation runs see also How to set up structural relaxation run NZ 3 1 How to set up atomic geometries 29 Parameter Value pgind 0 automatic search for symmetries 1 no symmetries assumed 1 C1 2 32 point group index pgind gt 1 defines the crystallographic point group as follows for more details consult pgsym routine Group pgind Group pgind Group pgind Group 1 C1 9 3m Cu 17 A mmm Dap 25 222 D2 T G 10 3m Dza 18 6 Ce 26 2mm Cov Ga 11 4 Cs 19 6 Cn 27 mmm D n m Cs 12 4 S4 20 6 m Cen 28 23 T 2 m Czn 13 4 m Can 21 622 De 29 m3 Th 3 C3 14 422 D4 22 6mm Cov 30 432 O 3 Se 15 4mm Ca 23 62m D n 31 43m Ta 32 Ds 16 42m Dza 24 6 mmm Den 32 m3m On Define the lattice parameters of the super cell parameter celldm in start inp The meaning of each celldm component depends on ibrav Usu ally celldm 1 contains the lattice constant a in bohr For hexago nal and rhombohedral super cells ibrav 4 and 10 respectively the c a as ai ratio is specified in celldm 2 When defining systems having A7 structure the common crystal phase of As Sb and Bi ibrav 10 the positions of the two atoms in the basis are given by u 0 0 c where the dimensionless parameter u should be pro
16. vanish Some comments for interested users e For a slab k point set those shells that contain contributions from lat tice vectors with a finite z component cannot vanish thus they must be disregarded when judging the quality of the basis set 3 3 Total Energy Minimization Schemes 38 e The quality assessment only makes sense for systems with a band gap The effect of having a sharp integration boundary the Fermi surface for a metal is not accounted for by Chadi and Cohen s shell analysis Note Even if there are no point group symmetries the vectors k and k are symmetry equivalent by virtue of time inversion symmetry For this reason the number of k points is reduced by at least a factor of two for any reasonable choice of a k point mesh 3 3 Total Energy Minimization Schemes In an electronic structure calculation using a plane wave basis the Hilbert space is typically spanned by a huge number of basis functions up to 10 plane waves Therefore it would be unwise to attempt to diagonalize the Hamiltonian operator in this high dimensional space directly Instead one uses algorithms which only imply vector operations on the wave function vector in Hilbert space rather than matrix operations The wave functions are gradually improved in an iterative process until they eventually converge towards the eigenvectors d Electronic minimization The goal is to minimize the total energy with respect to the wave function Wik sta
17. 0 0 025 pane Gamma 11414 fold parameter false k point coordinates relative or absolute 8 0 4 0 Ecut Ry Ecuti Ry 0 005 true false ekt tmetal tdegen true true 1 tmold tband nrho 5 2 1234 npos nthm nseed 873 0 1400 0 1e8 1 T_ion T_init Q nfi_rescale true true tpsmesh coordwave 1 3 Gallium number of atoms zv name 1 0 3 0 0 73 3 gauss radius mass damping 1_max 1_loc SER SCS ks t_init_basis 0 0 0 0 0 0 f taud tford 1 5 Arsenic number of atoms zv name 1 0 3 0 0 7 3 3 gauss radius mass damping 1_max 1_loc Ge ete ks t_init_basis 0 25 0 25 0 25 f taud tford e A sample output file fort 6 generated by fhi98md using input from the start inp file above fort 6 xkxkkk This is the complex fhi98md program Fo kk kk February 1999 k kk kk k gt gt gt electronic time step 20 0000 gamma 2000 gt gt gt Using delt 8 0000 gamma 2000 when energy changes less than 1000E 03 gt accuracy for convergency epsel 00001 epsfor 00050 epsekinc 10000 gt schultze algorithm for electron dynamics gt ions are not allowed to move normally no mixing of old charge is done N B The same atom should not be involved in more than one center of mass type constraint gt ibrav 2 pgind 0 nrot 24 alat 10 682 omega 304 7177 mesh 14 14 14 gt ecut 8 0 ryd ecuti 4 0 ryd nkpt 40 nel tmetal ekt tdegen 8 00000000000000000 T 0 50000000000
18. 00 1997 11 P Kratzer E Pehlke M Scheffler M B Raschke and U Hofer Phys Rev Lett 81 5596 1998 12 T K Zywietz J Neugebauer and M Scheffler Appl Phys Lett 74 1695 1999 13 G B Bachelet D R Hamann and M Schliiter Phys Rev B 26 4199 1982 14 D R Hamann Phys Rev B 40 2980 1989 15 N Troullier and J L Martins Phys Rev B 43 1993 1991 16 X Gonze R Stumpf and M Scheffler Phys Rev B 44 8503 1991 17 L Kleinman and D M Bylander Phys Rev Lett 48 1425 1982 18 D M Ceperley and B J Alder Phys Rev Lett 45 567 1980 19 J P Perdew and A Zunger Phys Rev B 23 5048 1981 20 A D Becke Phys Rev A 38 3098 1988 21 J P Perdew Phys Rev B 33 8822 1986 22 J P Perdew J A Chevary S H Vosko K A Jackson M R Pederson D J Singh and C Fiolhais 46 6671 1992 23 C Lee W Yang and R G Parr Phys Rev B 37 785 1988 24 J P Perdew K Burke and M Ernzerhof Phys Rev Lett 77 3865 1996 57 BIBLIOGRAPHY 58 25 26 27 28 29 30 31 32 M Bockstedte A Kley J Neugebauer and M Scheffler Comput Phys Com mun 107 187 1997 M Fuchs and M Scheffler Comput Phys Commun 119 67 1999 M C Payne et al Rev Mod Phys 64 1045 1992 M C Payne et al Phys Rev Lett 56 2656 1986 A Williams and J Soler Bull Am Phys Soc 32 5
19. 0000000000000000000000000000000n0000 FHI98md start utility ended normally different systems e Bulk GaAs fcc lattice basis GaAs can be regarded as a fcc lattice with the two point basis 0 and ai a2 az A straightforward way to set up a bulk calculation in this case is to use ibrav 2 and to allow for an automatic search for sym metries by setting pgind 0 The lattice constant a 10 682 a u is given as celldm 1 and the coordinates of the diatomic basis are there fore specified in units of a Rea 0 0 0 and Ras 0 25 0 25 0 25 Here is a sample input file start inp start inp 1 number of processors 1 number of minimal processors 1 number of processors per group 2 number of species 0 excess electrons 5 number of empty states 2 0 ibrav pgind 10 6 0 682 0 0 0 0000 celldm 1 number of k points 0 5 0 5 0 5 1 0 k point coordinates weight 444 fold parameter true k point coordinates relative or absolute 8 0 2 0 Ecut Ry Ecuti Ry 0 005 true false true false 1 5 2 1234 873 0 1400 0 1e8 1 ekt tmetal tdegen tmold tband nrho npos nthm nseed T_ion T_init Q nfi_rescale 3 2 Choice of the k point mesh 33 true true tpsmesh coordwave 1 3 Gallium number of atoms zv name 1 0 3 0 0 733 gauss radius mass damping 1_max 1_loc SG CRA SEs t_init_basis 0 0 0 0 0 0 op taud tford 1 5 Arsenic number of atoms zv name
20. 0000010E 02 F gt alat 10 682000 alat 10 682000 omega 304 717734 Lattice vectors al 5 341000 000000 5 341000 a2 000000 5 341000 5 341000 a3 5 341000 5 341000 000000 Reciprocal lattice vectors bi 1 000000 1 000000 1 000000 b2 1 000000 1 000000 1 000000 b3 1 000000 1 000000 1 000000 Atomic positions tau0 Species Nr x y z Galium 1 0000 0000 0000 Arsenic 1 2 6705 2 6705 2 6705 Number of kpoints nkpt 40 Ratios of FFT mesh dimensions to sampling theorem 1 029 1 029 1 029 gt ps pots as given on radial mesh are used gt gvk ngwx and max nr of plane waves 134 125 k point weight of g vectors 1 50 50 50 0010 120 2 44 44 44 0010 114 40 00 00 00 0010 113 Weigthed number of plane waves npw 116 221 Ratio of actual nr of PWs to ideal nr 99816 of electrons 8 0000 of valencestates 4 of conduction states 5 3 6 How to set up a band structure calculation atomic data for 2 atomic species pseudopotentialparameters for Galium gt nr of atoms 1 valence charge 3 000 force fac 3 0000 speed damp 7000 l_max 3 1_loc 3 rad of gaussian charge 1 000 pseudopotentialparameters for Arsenic gt nr of atoms 1 valence charge 5 000 force fac 3 0000 speed damp 7000 l_max 3 1_loc 3 rad of gaussian charge 1 000 Final starting positions 0000 0000 0000 phfac is n_ideal 1 0 phfac is n_ideal 2 0 phfac is i_kgv n_class 100 phfac is i_kgv n_class 200 oniskb is 1 wnl 1
21. 1 0 3 0 0 733 gauss radius mass damping 1_max 1_loc i ae ae a t_init_basis 0 25 0 25 0 25 f taud tford 3 2 Choice of the k point mesh For a periodic system integrals in real space over the infinitely extended system are replaced by integrals over the finite first Brillouin zone in reciprocal space by virtue of Bloch s theorem In fhi98md such integrals are performed by sum ming the function values of the integrand for instance the charge density at a finite number of points in the Brillouin zone called the k point mesh Choosing a sufficiently dense mesh of integration points is crucial for the convergence of the results and is therefore one of the major objectives when performing convergence tests Here it should be noted that there is no variational principle governing the convergence with respect to the k point mesh This means that the total energy does not necessarily show a monotonous behavior when the density of the k point mesh is increased d Monkhorst Pack mesh In order to facilitate the choice of k points the fhi98md package offers the possibility to choose k points according to the scheme proposed by Monkhorst and Pack 30 This essentially means that the sampling k points are distributed homogeneously in the Brillouin zone with rows or columns of k points running parallel to the reciprocal lattice vectors that span the Brillouin zone This option is enabled by setting t_kpoint_rel to true which should b
22. 5 WS J 412 0 18 0 3 GaAs 96 35 J 6 0 3 GaN 70 In As 40 WS 4 Ag 66 Ag 64 J 4 1 Ag 7 64 J 2 0 3 Ru 21 CO O 54 J 2 0 2 3 4 How to set up a structural relaxation run 43 3 4 How to set up a structural relaxation run 4 Structural optimization schemes The present program provides two methods to find the equilibrium geometry namely the modified steepest descent scheme and the damped Newton scheme In the damped Newton scheme the ions whose coordinates are given in the input file start inp are relaxed according to the iterative scheme geal Nee Sa ae Pa where Az and uzr are the damping and mass parameters respectively Note that these parameters determine whether the ions lose their initial potential energy slowly in an oscillatory like motion or whether they move straight into the closest local minimum as in the modified steepest descent scheme L Pre settings e Structural optimization run The program has switches tfor and tdyn in inp mod for setting up a structural relaxation a molecular dynamics calculation or an electronic structure calculation The following table shows how to set up a struc tural relaxation Parameter Value tfor tdyn true false relaxation of ionic positions c f epsfor force_eps and tford e Choosing the structural optimization scheme Set the parameter tsdp in inp mod to choose which structural optimiza tion scheme will be used Parameter Value tsdp true modified steepest desce
23. 535 occ 2 0000 2 0000 2 0000 2 0000 0000 0000 0000 0000 0000 k point 000 000 000 eigenvalues and occupation numbers eig 10 005 2 030 2 030 2 030 2 781 6 232 6 232 6 232 9 962 occ 2 0000 2 0000 2 0000 2 0000 0000 0000 0000 0000 0000 gt gt gt 1 10 8263451 7 7123788 time per tb loop 18 49 53 0 0 00000 00000 1 00000 00000 00000 1 00000 iasym nr of symmetric at xneu tau0 lt sym Op isym gt 3 6 How to set up a band structure calculation 54 total time for 1 loops 18 49 init stores starting density for mixing in c_fft_store time elapsed for init 26 8600 total number of wave function components 4652 LOOP n_it 1 time elapsed for nlrh t 0300 RHOE read from unit7 in RHOOFR Number of electrons in real space 7 99999630944237516 time elapsed for rhoofr t 0000 internal energy at zero temperature 8 305721 a u non equillibrium entropy 000000 kB equillibrium entropy 000000 kB kT energy 005 eV non eq free energy 8 305721 a u non eq total energy 8 305721 a u Harris energy 8 451691 a u kinetic energy 3 225369 a u electrostatic energy 10 000821 a u real hartree energy 763345 a u pseudopotential energy 1 472645 a u n l pseudopotential energy 632887 a u exchange correlation energy 2 370027 a u exchange correlation potential energy 3 090283 a u kohn sham orbital energy 739312 a u self energy 13 564038 a u esr energy
24. 62 1987 H J Monkhorst and J D Pack Phys Rev B 13 5188 1976 D J Chadi and M L Cohen Phys Rev B 8 5747 1973 S L Cunningham Phys Rev B 10 4988 1974 Index al 19 a2 19 a3 19 alat 19 ampre 6 amprp 6 atom 12 atoms 3 movement of 3 number of 20 42 b1 19 b2 19 b3 19 Brillouin zone 33 37 50 56 integrals in 37 symmetry lines in 50 broyd 5 celldm 12 constraints ini 47 constrw 47 coordwave 12 48 delt 6 delt2 6 delt_ion 6 dforce_b 5 ecut 12 ecuti 12 ekt 13 energy 26 eps_chg_dlt 7 epsekinc 7 epsel 7 epsfor 7 44 fermi 5 fhi98md 1 tiling strategy in 36 algorithms used in 3 structure of the program 3 fhi98pp package 23 fhi98start 1 6 49 fiondyn 5 59 fionsc 5 44 force_eps 1 7 44 force_eps 2 7 fort 1 26 fort 20 3 49 fort 21 26 49 fort 3 51 fort 6 24 26 fort 70 12 fort 71 26 48 fort 72 5 26 48 49 fort 73 26 fort 74 26 fort 80 26 gamma 7 gamma2 7 iedyn 8 i_facs 14 50 irc 8 iatm 47 ibrav 13 idyn 8 ineq_pos 19 init_basis 12 initbasis 8 initialization 3 20 block 3 mixed basis set 3 inp ini 3 18 37 48 inp mod 3 6 39 43 ion_damp 13 44 ion_fac 13 44 iprint 8 48 Lloc 14 24 latgen 31 Imag 14 24 maz_basis_n 19 max_no force 9 44 mesh 19 Chadi Cohen 35 Fourier 21 k point 33 38 INDEX Monkhorst Pack 33 41
25. 7 a u exchange correlation energy 2 370027 a u exchange correlation potential energy 3 090283 a u kohn sham orbital energy 640411 a u self energy 13 564038 a u esr energy 000048 a u gaussian energy 5 801080 a u 3 6 How to set up a band structure calculation 55 amp amp s atomic positions and local nl forces on ions Galium gt amp amp s n 000000 000000 000000 000000 000000 000000 Arsenic gt amp amp s n 2 670500 2 670500 2 670500 000000 000000 000000 gt sum of all local nl forces n_atoms 0000000000 0000000000 0000000000 should 0 nel dfmax fscale efermi ekt seq sneq 8 00000 030 40000 2 60571 005 00000 00000 gt 1 k point 500 500 500 ngw 120 gt Eigenvalues and Occupations zeig 8 164 4 000 1 350 1 350 3 484 7 111 7 111 10 332 12 784 gt occ 2 0000 2 0000 2 0000 2 0000 0000 0000 0000 0000 0000 gt 40 k point 000 000 000 ngw 113 gt Eigenvalues and Occupations gt eig 9 667 2 453 2 453 2 453 3 041 6 239 6 239 6 239 9 894 gt occ 2 0000 2 0000 2 0000 2 0000 0000 0000 0000 0000 0000 gt gt gt n_it nfi Ekinc Etot Eharr Ezero mForce mChange Seq Sneq Efermi gt gt gt 23 0 09107 8 29527 8 44641 8 29527 00000 000 0000 0000 2 6057 I finished storing wavefunctons and data on fort 71 END OF THE MAIN LOOP average time elapsed for nlrh rhoofr vofrho n x nkpt x dforce nkpt x graham ortho rest in main per elec time s
26. 8 unit cell whose size should be tested to satisfy the above condition see also Choice of the k point mesh e Slab calculations The ibrav value in this case should be chosen to reflect the symme try of the surface unit cell employed in the calculation Coordinates of the nuclei are specified as described above The size of the su per cell in the direction perpendicular to the surface should ensure a large enough vacuum region between the periodic slab images see also Choice of the k point mesh 1 Examples e The A7 crystal structure of As A sample start inp file to set up bulk calculation for the A7 structure ibrav 10 of As with lattice parameters atat 7 1 bohr c aja 2 67 and u 0 227 as specified in celldm 1 3 No scaling of the super cell will be performed celldm 4 6 1 0 1 0 1 0 In order to switch on the automatic search for symmetries pgind has been set to 0 see for example R J Needs R Martin and O H Nielsen Phys Rev B 33 3778 1986 3 1 How to set up atomic geometries number of processors number of minimal processors number of processors per group number of species excess electrons number of empty states 10 0 ibrav pgind 7 1 2 67 0 227 1 1 1 celldm 1 number of k points 0 25 0 25 0 25 1 0 k point coordinates weight 555 fold parameter oOonorrrrkr true k point coordinates relative or absolute 10 4 0 Ecut Ry Ecuti Ry 0 1 true f ekt tme
27. FHI98MD Computer code for density functional theory calculations for poly atomic systems alana User s Manual authors of this manual P Kratzer C G Morgan E Penev A L Rosa A Schindlmayr L G Wang T Zywietz program version 1 03 August 1999 Contents 1 Introduction 1 2 Description of the program structure and input files 3 2 1 The program structure aaao 3 2 2 The input files inp mod and start inp 0 0 00 0 000 6 2 3 The input fileinpini 2 0 0 000002 00000 18 2 4 Pseudopotential file s 0 0 0 2 000 23 2 5 Input files for advanced users 2 2 00 2 ee 25 2 5 1 The input file constraints ini 004 25 2 5 2 The input fileinp oce o s ocs s s o adt aE 00000202 ae 25 2 6 Runtime control files aoaaa 2 2 2 0 0 0020004 25 2T Fhe output files s u ee ee pe BS Pe A 26 3 Step by step description of calculational aspects 27 3 1 How to set up atomic geometries 2 2 20004 27 3 2 Choice of the k point mesh aoaaa 33 3 3 Total Energy Minimization Schemes 38 3 4 How to set up a structural relaxation run 43 3 5 How to set up acontinuationrun 48 3 6 How to set up a band structure calculation 49 Bibliography 57 Index 59 iii Chapter 1 Introduction Total energy calculations and molecular dynamics simulations employing density functional theory 1 represent a reliable tool in condensed
28. OP FRPN r OO 2 3 The input file inp ini The input file inp ini is usually generated automatically by the start utility fhi98start from the files inp mod start inp and constraints ini However the pro gram fhi98md also runs individually without the help of the start utility This requires the user to provide the file inp ini in addition to inp mod the pseudopo tentials and possible restart file s 2 3 The input file inp ini 19 QI inp ini The file inp ini contains processed data from start inp plus a copy of con straints ini All new quantities are described below in alphabetical order while unchanged parameters taken from either start inp or constraints ini are not listed again We refer instead to the detailed information given earlier Parameter al 1 3 a2 1 3 a3 1 3 Parameter alat Parameter b1 1 3 b2 1 3 b3 1 3 Parameter ineg_pos 1 3 Parameter max_basis_n Parameter MMNALL Parameter n fft store Type 3 x real 3 x real 3 x real Type real Type 3 x real 3 x real 3 x real Type 3 x integer Type integer Type integer Type Value integer 1 lattice vectors in a u lattice constant in a u reciprocal lattice vectors in a u if t coord auto true number of mesh points along the corresponding lattice vec tor of the super cell the mesh must con tain the origin 0 0 0 of the super cell connected to t coord_auto 1 3 max_
29. The file start inp describes the geometry of the super cell and the configuration of the nuclei It also contains information that is relevant for the MD simulation or the structure optimization and the calculation of the electron ground state 1 inp mod The file inp mod contains mainly the control parameters for the runs like e g the time steps for the electronic and atomic minimization schemes the convergence criteria and maximum number of steps etc In the following the parameters are explained in alphabetical order Parameter Type ampre real an amplitude of a random perturbation is added to the wave function but only if the parameter trane is set to true see also parameter trane Parameter Type amprp real an amplitude of a random perturbation is added to the ionic positions but only if the parameter trane is set to true see also parameter trane Parameter Type delt real step length of the electronic iteration this value has to be individually optimized in order to obtain optimal convergence see also parameter delt2 Parameter Type delt2 real second step length of the electronic itera tion connected to parameter eps_chg_dlt Parameter Type only relevant for MD simulations delt_ion real time step for the integration of the ionic equations of motion in a u 2 2 The input files inp mod and start inp Parameter epsekinc epsel epsfor Parameter eps_chg_dlt Parameter force_eps force_eps 1 fo
30. arameter Type only for MD simulations if tdyn is true T_init real temperature of initial velocities in K only active if npos is either 3 or 5 connected to npos Parameter Type only for MD simulations if tdyn is true T_ion real ionic temperature in K only active if nthm is either 1 or 2 connected to nthm 2 2 The input files inp mod and start inp 17 Parameter Value t_kpoint_rel false k points are given in Cartesian coordi nates in units 27 aja true frame of reference for k points is spanned by the reciprocal lattice vectors i e k ky b kobe k3b3 Parameter Type tau0 1 3 real ionic coordinates units depend on ibrav flag whether atoms may move given by tford connected to tford Parameter Type tband logical if set to true the k point set is not re duced by fhi98start utility the electron density is not recalculated after the ini tialization Parameter Type tdegen logical if set to true the initial occupation numbers are read in from file inp occ kept fixed for the run Parameter Type Value tmetal logical occupation of the electronic states true Fermi distribution see also ekt false step like distribution connected to ekt Parameter Type tmold logical if set to false only the initialization is performed Parameter Type Value tpsmesh logical form of the pseudopotential true tabulated on logarithmic mesh false set up from parameterized form 2 3 The input file i
31. as rough guide as everyone has to find her his optimum parameters e For an example of a geometry optimization calculation see Section How to set_up a structural relaxation run 5 T delt 5 see delt 10 delt 20 delt 40 In E E Hartree UR oO 15 1 Au Ni 1 0 20 40 Electronic time steps Figure 3 4 A total energy calculation for GaAs bulk in the zinc blende structure within LDA approximation by using the Williams Soler minimization scheme to iterate the wave functions Different values for the electronic time step delt are used The Monkhorst Pack k points mesh was taken as 4x4x4 with the initial k point 0 5 0 5 0 5 and an energy cutoff of 8 Ry was used Eo is the limiting value of the total energy approached in a very well converged run 3 3 Total Energy Minimization Schemes 42 Table 3 1 Electronic time step delt damping parameter gamma and minimization schemes used for some bulk and slab with and without adsorbates calculations The labels J WS and SD mean that either the damped Joannopoulos or Williams Soler or Steepest Descent minimization scheme was employed The dimension along the z direction is given by dim z Bulk Material Number of atoms Scheme delt gamma Si 2 J 30 0 4 GaAs 2 J 20 0 3 GaAs 64 defect J 6 0 3 GaN 4 WS 4 Al 1 WS 6 Al 1 SD 12 Ag 1 J 5 0 2 Pt 1 J 2 0 3 Slab Material aoa Adsorbates ae Scheme delt gamma Si 70 H 43 J 12 0 4 GaAs 70 In 3
32. ay that a maximum number of symmetries is met The fhi98md code is distributed together with the fhi98start utility which helps to search for relevant symmetries and to set up crystals and slabs from some standard crystallographic symmetry classes Under UNIX environment for example the fhi98start utility is usually invoked by the same shell script that runs fhi98md below a protocol of the fhi98start run is saved in the file start out bin csh xvf cp FHI98MD fhimd fhi98md cp FHI98MD start fhi98start HEHHEHHHHHRHRRARAHRHR AERA HARA REAR AH RH RAH RBH RH REAR H run fhi98start program and build up inp ini in principle one could create inp ini by hand but fhi98start gives a consistent and optimized input for fhi98md HHEHHHHHHHHEHEHEHEHEHHAHEHEHEHAHBHBH RH RH EHR ER HEH R EHH R ES fhi98start tee start out HEBHRHHHHRRHEHHRAHRHREHAHEB HABER HR EERE EAH EB BH BH ER BBE run f hi98md HEBHRHHHHRRHEHHREHRHREHAHEBEABHAHHR EHR EAH ER EA BHR ERE E hi98md 1 Basic input parameters All input information about crystal structure is produced by the fhi98start utility according to the values of the following parameters specified in the file start inp 27 3 1 How to set up atomic geometries 28 Parameter Type ibrav integer Bravais lattice pgind integer point group index celldm 1 6 6 x real lattice parameters of the super cell depends on ibrav tau0 1 3 3 x real coordinates of the nuclei
33. basis_n max nz nz_basis connected to nz and na_basis maximum size of pseudo potential grid option currently not implemented 2 3 The input file inp ini Parameter nar Parameter nel Parameter ngwix Parameter ngwr Parameter nlmax Parameter nlmaxinit Parameter NNT Parameter nrot Type integer Type real Type integer Type integer Type integer Type integer Type integer Type integer 20 maximum number of atoms per species number of electrons maximum number of plane waves during the initial diagonalization should obey ngwix gt in x omega X ecuti connected to omega and ecuti maximum number of plane waves used to represent the wave functions should obey ngwr gt an x omega x ecut connected to omega and ecut nlmax max l maz 2 x lloc 1 connected to maz and lloc maximum number of atomic orbitals per atom in the initialization nnrz nrz 1 1 x nrz 2 x nrz 3 connected to nra 1 3 number of point symmetries connected to s 3 3 2 3 The input file inp ini Parameter nrz 1 3 Parameter nschitz Parameter NST Parameter nT Parameter nx_basis Parameter ng init Parameter omega Parameter s 8 8 Type integer Value Type integer Type integer Type integer Type integer Type real Type real 21 dimensions of the Fourier mes
34. be smaller than this criterion force_eps 1 0 001 0 1 maximum allowed relative variation in lo cal forces before if tfor true and tdyn false executing a geometry op timization step or if tfor false and tdyn true calculating total forces e Other parameters timequeue is the maximum cpu time in seconds nomore is the maxi mum number of atomic moves nstepe is maximum number of electronic iterations allowed to converge forces If the forces on the ions do not converge to the criterion force_eps 1 after nstepe iterations then the program terminates The maximum number maz_no _force of electronic iterations for which no local forces shall be calculated after each ionic step move should also be set 4 Example e The example presented here is the structure optimization of the noornerre 8 0 GaAs 110 cleavage surface The start atomic coordinates are at the ideal bulk positions The files start inp and inp mod are given below start inp number of processors number of minimal processors number of processors per group number of species excess electrons number of empty states ibrav pgind 7 403408 10 47 51 823856 0 0 O celldm BOR number of k points 0 5 0 0 1 0 k point coordinates weight 5 41 fold parameter 3 4 How to set up a structural relaxation run 45 true 10 4 0 0 004 true f true false 1 k point coordinates relative or absolute Ecut Ry Ecut
35. cation Furthermore the choice of available gradient corrected functionals has been increased The package consists of the program fhi98md and a start util ity fhi98start The program fhi98md can be used to perform static total energy calculation or ab initio molecular dynamics simulations The utility fhi98start assists in generating the input file required to run fhi98md thereby ensuring the lowest possible memory demand for each individual run Thus no recompilations are required a full calculation can be performed by calling the two binary executables fhi98start and fhi98md in sequence This manual consists of two parts The first part is a reference list to look up the function of input parameters which are described in alphabetical order The second part describes typical procedures of how one might use the code Here we describe the parameters in their functional context Also the way in which the input parameters control the algorithms used in the code is described in more detail Finally we provide examples of how to use the code to solve a particular physical problem Chapter 2 Description of the program structure and input files 2 1 The program structure For detailed information about the algorithms used in fhi98md we refer the inter ested reader to the article 25 The structure of the program fhi98md is sketched in Fig 2 1 The self consistent calculation of the electron ground state forms the main body of the program wh
36. ce may be allowed to relax in the ry plane but not along the z direction An adsorbing molecule may be allowed to rotate and to move in the zy plane but its center of mass may not be allowed to move in the z direction When such constraints on the atomic motion are being set up the location of the center of mass of a group of atoms is determined by atomic weights chosen by the user These atomic weights do not have to be proportional to the actual atomic masses of the atoms The information concerning the desired constraints on the atomic motion is given in the file constraints ini which is used by the start program The first line of this file contains the number of single atom constraints nr_constraints One line for each of these constraints must follow containing the vector direc tion along which the atom is not allowed to move vec j j 1 3 followed by the species number is and the atom number ia In the sample constraints ini file below we see that three single atom constraints are being set up Atom number 31 of species number 1 is not allowed to move in the 110 direction and atom number 32 of species number 1 is not allowed to move in the 001 direction or the 010 direction The second part of the file constraints ini must contain a line with the number of center of mass type constraints nconstr followed by a short section for each constraint The section for constraint i begins with a line giving the number of ato
37. d k points and weights k points and weights folding factors frame of reference for k points Here we have made use of the tiling strategy employed in fhi98md An even denser k point set consisting of 54 points in the full Brillouin zone may be obtained by using the 6 k points of the first example but as a pattern 3 2 Choice of the k point mesh 37 Figure 3 3 2D Brillouin zone of a fec 111 surface with hexagonal symmetry with set of 18 special k points following Cunningham The thin polygon indicates the conventional first Brillouin zone the thick polygon marks the Brillouin zone as realized in the fhi98md code repeated in each of the nine tiles i e by setting the folding parameters in the first example to 3 3 1 L User supplied k point sets In some cases like a band structure calculation the user might find it more convenient to specify the k point mesh directly with respect to the coordinate axes in reciprocal space rather then with respect to the reciprocal lattice vectors This can be achieved by setting t_kpoint_rel to false The unit of length on the coordinate axes is 27 djaz in this case The folding parameters can be used as well to enhance the number of k points if desired However one should keep in mind that the k point sets specified in that way might have little symmetry i e their number is not significantly reduced by the built in symmetry reduction algorithm of fhi98start a Reduced k p
38. e the default for total energy calculations The density of k points can be chosen by the folding parameters i_facs 1 3 With these parameters you specify to cover the entire Brillouin zone by a mesh of i_facs 1 x i_facs 2 x i_facs 3 x nkpt points The details of this procedure are as follows In fhi98md the Brillouin zone is spanned by the reciprocal lattice vectors b b and bg attached to the origin of the coordinate system According to this definition one corner of the Brillouin zone rests in the ori gin The entire Brillouin zone is tiled by small polyhedra of the same shape as the Brillouin zone itself The parameters specify how many tiles you have along the b b and bg direction In each tile you specify k points supplied in form of a list The coordinates of these k points are given relative to the spanning vectors of a small polyhedron or tile i e k zk 1 b xk 2 b2 zk 3 bs The supplied k point pattern is then spread out over the whole Brillouin zone by translations of the tile In other words the k point pattern of a smaller Brillouin zone which would correspond to a larger unit cell in real space is unfolded in the Brillouin zone of your system under study Normally the pattern consists only of a single point in the center of the tile leading to the conventional Monkhorst Pack k point sets e k point set for a bulk calculation A k point set typically used in a bulk calculation could look like 3 2
39. enic 1 0000 0000 4 3032 Arsenic 2 0000 0000 4 3032 gt alat 7 100000 alat 7 100000 omega 275 864416 lattice vectors al 4 099187 000000 6 319000 a2 2 049593 3 550000 6 319000 a3 2 049593 3 550000 6 319000 number of symmetries of bravais lattice 12 number of symmetries of bravais latt at basis 12 symmetry matrixes in lattice coordinates gt start out 3 1 How to set up atomic geometries 32 1 0 0 0 1 0 0 0 1 12 1 0 0 0 0 1 0 1 0 centered atomic positions gt is ia positions 1 1 00000 00000 4 30324 1 2 00000 00000 4 30324 reciprocal lattice vectors b1 1 154701 000000 374532 b2 577350 1 000000 374532 b3 577350 1 000000 374532 mesh parameter 1 ideal 15 1639 used 2 ideal 15 1639 used 16 ratio 1 0551 3 ideal 15 1639 used 16 ratio 1 0551 The k point coordinates are assumed to be relative absolute k point coordinates in 2Pi alat 16 ratio 1 0551 k points weight 1 00 00 06 0080 125 00 00 96 0080 Using the existing symmetries the set of k points can be reduced to k points weight 1 0000000 0000000 0561798 0080000 35 0000000 0000000 9550562 0080000 gt Shell analysis of quality of k points after Chadi Cohen gt Number of A_m 0 shells N 63 gt Weighted sum of A_ms 0002112 gt should be small and gives measure to compare List of a_m s O zero x changing n nonzero 00000000000000000000000000000000
40. entered cubic lattice L T X W K T 40 equidistant k points nkpt 40 Note that with the number of empty states nempty 5 the five lowest unoccupied bands will be calculated as well By setting t kpoint rel false one guarantees that the k point coordinates will be used exactly as specified start inp number of processors number of minimal processors number of processors per group number of species excess electrons number of empty states ibrav pgind 0 000 celldm number of k points 0 5 0 5 0 025 t ieaiai 0 4422650 0 4422650 0 025 4 0 3845300 0 3845300 0 025 0 3267950 0 3267950 0 025 P 0 2690600 0 2690600 0 025 a 0 2113250 0 2113250 0 025 t 0 1535900 0 1535900 0 025 h 0 0958548 0 0958548 0 025 0 0381198 0 0381198 0 025 y i 0 0 0 025 Gamma n 0 0 0 025 0 0 0 025 t 0 0 0 025 h 0 0 0 025 e 0 0 0 025 0 0 0 025 B 0 0 0 025 r 0 0 0 025 i 0 0 0 025 1 0 0 0 025 S 1 0 1 0 0 025 o 0 2 0 0 025 i 0 3 0 0 025 n 0 4 0 0 025 0 5 0 0 025 W Z 0 5707110 0 0 025 o 0 6414210 0 0 025 n 0 7121320 0 0 025 3 e 0 75 0 0 025 ese K 0 6792890 0 0 025 0 6085790 0 0 025 0 5378680 0 0 025 Dvolt 0000 3 6 How to set up a band structure calculation 52 0 4671570 0 4671570 0 0 025 0 3964470 0 3964470 0 0 025 0 3257360 0 3257360 0 0 025 0 2550250 0 2550250 0 0 025 0 1843150 0 1843150 0 0 025 0 1136040 0 1136040 0 0 025 0 0428932 0 0428932 0 0 025 0 0
41. fort 11 2 gt fort 12 nsp gt fort Insp In practice the pseudopotential data are copied to the working space by the same shell script batch file under Windows which runs the fhi98md program bin csh xvf HHHHHHHHHHHHHH set directories HHHRHRHHHRHRHREREE set PSEUDO set WORK pseudopotentials directory working directory HHHHHHHHHHH change to the working space HHH cd WORK move pseudopotentials to the working space cp PSEUDO ga lda ham cpi fort 11 cp PSEUDO as lda ham cpi fort 12 HEBHRHHHHRRHEHHRAHRHREHAHER HABER HR BARE EAH ER BH BH ERE RH run f hi98md HEBHRHHHHRRHEHHRAHRHRAHRRHERERBHRBHREHRHEAHER BHR ER RH hi98md Pseudopotentials are read in during the initialization when control is delivered to the init routine The tpsmesh parameter in start inp allows the user to instruct the fhi98md about the format of the pseudopotential data file s Parameter Type Value tpsmesh logical form of the pseudopotential true tabulated on logarithmic mesh false set up from parametrized form Ionic pseudopotentials in the format accepted by the fhi98md program can be generated and tested by means of the fhi98pp package The latter provides the psgen tool which produces as chief output a pseudopotential data file name cpi formated as shown in Fig 2 2 for more details the user is referred to the article 26 1The current version of init supports up to 6 pseudopotential files Customize
42. gence to energy potential and forces move atoms structure optimization or dynamics have atoms been moved compute recalculate dforce G k Arsh O Tk struct structure and phase factors iterate wave functions U k piifac IP k P k calculate nlrhkb non local contributions repeat for all states to energy potential and forces vofrho ewald sum and local contributions perform to energy potential and forces graham Grahm Schmidt orthonormalisation wet ke WY 9 k repeat for all k points calculate e g for metals fermi occupation according to fermi distribution validate iteration parameters c f files delt stopfile and stopprogram terminate on convergence or if number of iterations or cpu time is exceeded KarI terminate Figure 2 1 Flow chart of the program fhi98md Output is generated at the end of each self consistency cycle and by the routines fiondyn fionsc fermi init and vofrho Restart files are written by the routine fiondyn and by a call of routine o_wave in the main program 2 1 The program structure 5 running These parameters are updated from the files stopfile remaining number of electronic iterations and stopprogram remaining number of structure optimization steps If these files are empty the parameters are not changed Finally the conver gence criteria are checked The program terminates when convergence is
43. gind 10 68 0 00 0 111 celldm 1 number of k points 0 5 0 5 0 5 1 0 k point coordinates weight 444 fold parameter true k point coordinates relative or absolute 8 4 0 Ecut Ry Ecuti Ry 0 1 true fi ekt tmetal tdegen true false 1 tmold tband nrho 5 2 1234 npos nthm nseed 873 0 1400 0 10000000 1 T_ion T_init Q nfi_rescale true false tpsmesh coordwave 1 5 Arsenic number of atoms zv name 1 0 74 92 0 7 3 3 gauss radius mass damping 1_max 1_loc Sts Sts ots t_init_basis 0 0 0 00 0 f tau0 tford 1 3 Gallium number of atoms zv name 3 3 Total Energy Minimization Schemes 41 1 0 69 72 0 7 3 3 gauss radius mass damping 1_max 1_loc PE T OE t_init_basis 0 25 0 25 0 25 f tau0 tford inp mod 1 100 1000000 nbeg iprint timequeue 100 1 nomore nomore_init 10 0 0 2 delt gamma 4 0 0 3 0 0001 delt2 gamma2 eps_chg dlt 400 2 delt_ion nOrder 0 0 1 0 pfft_store mesh_accuracy 22 idyn i_edyn 0 false i_xc t_postc F 0 001 F 0 002 trane ampre tranp amprp false false false 1800 tfor tdyn tsdp nstepe false tdipol 0 0001 0 0005 0 3 epsel epsfor epsekinc 0 001 0 001 3 force_eps max_no_force 2 init_basis e We attach a table containing some values to the electronic time step delt the damping parameter gamma and the minimization method adopted which have been successfully used Of course this table provides only a few numbers and should be used only
44. gy fort 1 fort 6 report txt status txt contents number of iterations total energy in Hartree at finite elec tronic temperature Harris energy in Hartree atomic relaxations and forces general output summary of the run current state of the calculation error messages L Output files for data analysis using the EZWAVE graphical user interface fort 80 rhoz for visualization of the potentials charge density along the the z axis i e at x y 0 1 Binary output files used for restarting the program or for analysis of the run with utility pro grams fort 21 fort 71 fort 72 fort 73 fort 74 when performing molecular dynamics positions and ve locities of the atoms along the trajectory complete restart information including all wavefunctions electronic charge density total effective potential electrostatic potential Chapter 3 Step by step description of calculational aspects 3 1 How to set up atomic geometries Before an electronic structure calculation can be performed it is necessary to specify a starting geometry for the atomic structure of the system we want to study If this geometry is invariant under certain discrete crystallographic symmetry operations the computational load can be reduced considerably by exploiting these symmetries Therefore it is recommended to analyze the symmetry of the atomic geometry before starting the calculation and to choose the unit cell in such a w
45. h should obey nre 1 gt 2m7 al ecut etc connected to a 1 3 and ecut control flag for efficient storage allocation should be set to zero unless i_edyn 2 damped Joannopoulos algorithm disabled all iteration schemes enabled connected to i_edyn in inp mod maximum number of atomic species maximum number of electronic states per k point maximum number of atomic orbitals in the initialization ng_init ngwix nx_basis connected to ngwiz and nz_basis super cell volume in a u connected to al 1 3 a2 1 3 a3 1 3 3x3 point symmetry matrices in units of the lattice vectors connected to nrot 2 3 The input file inp ini Parameter t_coord_auto 1 3 1 Example 22 Type logical if tford false enable more efficient treatment of fixed ions on a mesh com mensurate with the super cell connected to tford and ineq_pos 1 3 The input file inp ini below was generated from the sample files start inp and constraints ini for metallic Ga given earlier The parameters are arranged in the order expected by the main program fhi98md inp ini 1 4 11 739 5912 nsx nax nx ngwxt 1 ngwx 8 8 102 118 24 24 24 1 ngwix nx_init nrx 1 nrx 2 nrx 3 nschltz 16 16 4 nx_basis max_basis_n nlmax_init 14400 39 4 570 nnrx nkpt nlmax mmaxx n_fft_store 1 1 1 minpes ngrpx nrpes 12 0 ibrav pgind 12 0000 T 0 10000 F nel tmetal ekt tdegen 20 00000 5 00000 ecut ecuti T F 1 tmold tband n
46. i Ry ekt tmetal tdegen tmold tband nrho 5 2 1234 npos nthm nseed 873 0 1400 0 1e8 1 T_ion T_init Q nfi_rescale t true tpsmesh coordwave 9 3 Gallium 1 0 1 000 0 733 number of atoms zv name gauss radius mass damping 1_max 1_loc eGirGite ct t_init_basis 0 00000 0 00000 14 806816 sb 3 701704 5 23500 11 105112 Stik 0 00000 0 00000 7 403408 st 3 701704 5 23500 3 701704 f 0 00000 0 00000 0 00000 sE 3 701704 5 23500 3 701704 f 0 00000 0 00000 7 403408 a 3 701704 5 23500 11 105112 ste 0 00000 0 00000 14 806816 SEs 95 Arsenic number of atoms zv name 1 0 1 000 0 7 3 3 gauss radius mass damping 1_max 1_loc SGe Steck t_init_basis 3 701704 2 61750 14 806816 st 0 00000 7 85250 11 105112 Sy 3 701704 2 61750 7 403408 te 0 00000 7 85250 3 701704 f 3 701704 2 61750 0 00000 e 0 00000 7 85250 3 701704 a 3 701704 2 61750 7 403408 SEx 0 00000 7 85250 11 105112 at 3 701704 2 61750 14 806816 t 1 20 1000000 inp mod nbeg iprint timequeue 1000 1 nomore nomore_init 2 0 0 2 delt gamma 1 0 0 2 0 0001 delt2 gamma2 eps_chg_dl1t 400 2 delt_ion nOrder 0 0 1 0 pfft_store mesh_accuracy 22 idyn i_edyn 0 false i_xc t_postc F 0 001 F 0 002 true false false false 0 0001 0 0005 0 1 0 001 0 001 3 1 4 Analysis trane ampre tranp amprp 100 tfor tdyn tsdp nstepe tdipol epsel epsfor epsekinc force_eps max_no_force init_bas
47. ich is displayed on the left hand side of Fig 2 1 The movement of the atoms is accomplished in the block move atoms which is sketched on the right hand side of Fig 2 1 Note that the generation of output is not explicitly accounted for in the flow chart and we refer to it at the end of this section The first block in the flow chart is the initialization block where the program reads the input files inp mod inp ini and the pseudopotential data Then the rou tines calculating form factors structure factors and phase factors are called and the initial wave functions are set up either from a restart file or by a few self consistency cycles using the mixed basis set initialization Having obtained the initial wave func tion pw the program enters the self consistency loop First the electron density and the contributions to energy potential and forces are calculated Note that the forces are only calculated during MD simulations and structure optimization when the electrons are sufficiently close to the Born Oppenheimer surface Within the block move atoms the atomic equations of motion EOMs are integrated for one time step in a MD simulation or a structure optimization is performed provided the electrons are sufficiently close to the Born Oppenheimer surface i e the forces are converged The control over the calculation of the forces is handled by this block If the nuclei have been moved i e either the atomic EOMs have been i
48. is After a job is executed you have to check whether it has been finished properly and the program has arrived at the optimized geometry You should check e The difference between the total energy and the Harris energy The difference should be very small e g lt 1075 Hartree e CPU time limit Structural relaxation usually takes a lot of CPU time and in most cases it stops because the CPU time limit is out e The total energy convergence see Fig 3 5 e Mean force on ions see Fig 3 6 e Mean displacement per coordinate see Fig 3 7 3 4 How to set up a structural relaxation run 257 38 257 39 J o 2 wo gt gt oO io 257 4 4 257 41 l l l 0 500 1000 1500 iterations Figure 3 5 Total energy as a function of electronic self consistency iterations 0 01 0 008 p 4 0 006 3 0 004 mean force Hartree a u o oo fb A 0 50 100 150 200 ionic steps Figure 3 6 Mean force on ions as a function of ionic steps 46 3 4 How to set up a structural relaxation run 47 0 008 4 0 006 J 0 004 mean displacement a u o oo2 A ae 0 50 100 150 200 ionic steps Figure 3 7 Mean displacement per coordinate as a function of ionic steps 1 constraint relaxation The present program also allows the user to constrain the motion of the atoms during a structural optimization For example an atom adsorbing on a surfa
49. is option would allow the user to modify the ionic positions slightly between consecutive runs Note that it is very efficient to use a previous wavefunction file as fort 70 when increasing the cut off In this case only the additionally required plane wave coefficients with larger wave number need to be recalculated When doing so however it is inevitable to rerun fhi98start to generate a modified inp ini file 3 6 How to set up a band structure calculation 49 1 Continuation run for a molecular dynamics simulation e Proceed like described above in particular set coordwave true e Rename the file fort 21 into fort 20 e In the file inp ini or start inp set npos to 6 This causes the program to read all required information to continue the ions trajectories from the binary file fort 20 L Continuation run with increased k point set or band structure cal culation e In the file inp ini set the parameter nrho to 3 e Make sure that the old output charge density in fort 72 is available In this case the wavefunctions need to be recalculated therefore nbeg remains to be equal to 1 Also the program fhi98start needs to be rerun 3 6 How to set up a band structure calculation l Pre settings Before running fhi98md in band structure mode do an ordinary total energy run for the system of interest to get a well converged electron density The latter is stored in fort 72 The unit cell should be chosen as follows e B
50. matter physics material science chemical physics and physical chemistry A large variety of applications in systems as different as molecules 2 3 bulk materials 4 5 6 7 and sur faces 8 9 10 11 12 have proven the power of these methods in analyzing as well as in predicting equilibrium and non equilibrium properties Ab initio molec ular dynamics simulations enable the analysis of the atomic motion and allow the accurate calculation of thermodynamic properties such as the free energy diffusion constants and melting temperatures of materials The package fhi98md described in this paper is especially designed to in vestigate the material properties of large systems It is based on an iterative approach to obtain the electronic ground state Norm conserving pseudopoten tials 13 14 15 16 in the fully separable form of Kleinman and Bylander 17 are used to describe the potential of the nuclei and core electrons acting on the valence electrons Exchange and correlation are described by either the local density approximation 18 19 or various generalized gradient approxima tions 20 21 22 23 24 The equations of motion of the nuclei are integrated using standard schemes in molecular dynamics Optionally an efficient structure optimization can be performed by a damped dynamics scheme The package fhi98md is based on a previous version fhi96md 25 The new version however is based on FORTRAN9O and allows dynamic memory allo
51. mesh_accuracy 9 minpes 15 mmaxx 19 nfft store 19 na 14 natm 47 naz 20 nbeg 9 48 nconstr 47 nel 20 nel_exc 14 nempty 14 50 nempty 25 nfi_rescale 14 ngrpx 15 ngwiz 20 ngwr 20 nkpt 14 nlmaz 20 nlmaz_init 20 nlrhkb_b 5 nnrz 20 nomore 9 25 nomore_init 9 nOrder 10 npes 15 npos 15 nr_constraints 47 nrho 12 15 48 49 nrot 20 nrz 21 nschltz 21 nseed 15 nsp 15 23 nstepe 10 44 nsx 21 nthm 16 nz 21 25 nx_basis 21 nx_init 21 o_wave 5 omega 21 pfft store 10 pgind 16 pgsym 29 project 5 pseudopotentials 1 19 23 format of file s 23 ionic 23 psgen tool 23 60 Q 16 report txt 26 rho_psi 5 rhoz 26 run 1 band structure 49 51 56 continuation 25 48 for a MD simulation 49 with increased k point set 49 structural relaxation 43 8 8 3 21 start inp 6 12 30 43 start out 31 37 status txt 26 stopfile 3 25 stopprogram 5 25 t_coord_auto 22 T init 16 tLinit_basis 16 T_ion 16 t_kpoint_rel 17 50 t_postc 10 tau0 17 31 43 tband 17 tdegen 17 25 tdipol 10 tdyn 10 43 tfor 11 13 43 tford 43 timequeue 11 tmetal 17 tmold 17 tpsmesh 17 23 trane 6 11 tranp 11 tsdp 11 13 43 vau0 18 vec 47 vofrho 5 wkpt 18 tk 18 50 zu 18 24
52. ms which are involved in this constraint natm i followed by the vector direction along which the center of mass of this group of atoms is not allowed to move One line follows for each atom in the group giving the weight of this atom in determining the center of mass position constrw i k and the species and atom number which are contained in the array iatm i k 1 l 1 2 In the example file below there is one center of mass type constraint which involves the motion of two atoms The center of mass of this two atom group is not allowed to move in the 001 direction The two atoms in this group are atoms number 29 and 30 of species number 2 and both are 3 5 How to set up a continuation run 48 weighted equally in determining the center of mass position We note that if a particular atom is included in a center of mass type con straint this atom should not be included in any other constraints involving a different group of atoms This has not been a serious restriction in the past since the most useful applications tend to involve at most one center of mass type constraint constraints ini 3 imr_constraints 1 1 0 ali 31 vec j j 1 3 is ia 0 t bE 32 vec j j 1 3 is ia 0 1 0 1 32 vec j j 1 3 is ia 1 mconstr 2 0 0 t natm i constrv i j j 1 3 0 5 2 29 constrw i k iatm i k 1 1 1 2 0 5 2 30 constrw i k iatm i k 1 1 1 2 3 5 How to set up a continuation run
53. np ini 18 Parameter Type only for MD simulations if tdyn is true vau0 1 3 real ionic velocities in a u only active if tdyn is set to true and npos is either 1 or 4 connected to tdyn and npos Parameter Type wkpt real weights of the k points given by xk 1 3 connected to rk 1 3 nkpt i_facs and t_kpoint_rel Parameter Type wk 1 8 real k point coordinates connected to wkpt nkpt ifacs and t_kpoint_rel Parameter Type zu real valence charge of the pseudopotential to be supplied for each element L Example The sample input file start inp below sets up a typical bulk calculation for metallic Ga The parameters are given in the order re quired by the fhi98md program start inp npes number of processors npesmin number of minimal processors npespg number of processors per group number of species excess electrons number of empty states 0 ibrav pgind 682 0 0 0 0000 celldm 6 number of k points 50 5 0 5 0 025 k point coordinates weight 11 fold parameter false k point coordinates relative or absolute 4 0 Ecut Ry Ecuti Ry 05 true false ekt tmetal tdegen true true 1 tmold tband nrho 5 2 1234 npos nthm nseed 873 0 1400 0 1e8 1 T_ion T_init Q nfi_rescale true true tpsmesh coordwave 3 Gallium number of atoms zv name O 3 0 0 7 3 3 gauss radius mass damping 1_max 1_loc CEs Aste ts t_init_basis 0 0 0 0 0 0 Le tau0 tford NORRPR FE P
54. nt scheme cur rently not supported false damped Newton scheme e Initial ionic coordinates and flag The initial ionic coordinates for all ions in the supercell have to be given and the flag which determines whether the ions may move has to be set Parameter Type Value tau0 1 3 3 x real ion coordinates tford true the ions move false the ions do not move e Mass and damping parameters The mass and damping parameters ion_fac and ion_damp determine how the ions lose their initial potential energy c f structural optimiza tion schemes 3 4 How to set up a structural relaxation run 44 Parameter Value ion_fac 0 3 2 mass parameter tfor true Note it is the ionic mass in amu if tdyn true and tfor false ion damp 0 4 0 7 damping parameter if tfor true and tdyn false Note itis used only when tsdp false e Force convergence parameters The calculation of the electron ground state for a fixed configuration of the ions is terminated if the improvement of the total energy and the wave functions per iteration is smaller than epsel and epsekinc respectively However for a structural relaxation calculation the residual forces acting on the ion should also be sufficiently small smaller than epsfor The error of the forces force_eps 1 is monitored in the subroutine fionsc before performing a structure optimization step Parameter Value epsfor 0 0005 forces on ions with tford true must
55. ntegrated for another time step or one structure optimization step has been executed it also recalculates the structure and phase factors and other quantities that explicitly depend on the positions of the nuclei Upon the first call to the routine fiondyn in this block the restart file fort 20 is read if provided and all necessary steps are taken to restart or initialize the dynamics The following two blocks update the wave functions using the damped Joannopoulos or the Williams Soler algorithm see section 3 3 and ortho normalize the wave function by a Grahm Schmidt ortho normalization In the last block within the self consistency loop the occupation numbers are updated e g for a metallic system according to a Fermi distribution using a mixing scheme between occupation numbers of subsequent electronic iterations This block also enables an interactive control over the remaining numbers of iterations while the program is 3 2 1 The program structure 4 initialization header read inp mod init read inp ini and pseudo potentials calculate structure and form factors formf nlskb phfac struct set up initial electron density and wave functions calculate call rhoofr electron density n r fionsc for structure optimization nlrhkb non local contributions fiondyn molecular dynamics to energy potential and forces both routines monitor the error of the forces and vofrho local contributions move the ions after conver
56. oints and symmetry Apart from the translational symmetry of the Bravais lattice the crystal struc ture under investigation may often have additional point group symmetries These can be used to reduce the number of k points which are needed in the actual calculation and thus the memory demand substantially To perform the integrals in the Brillouin zone it is sufficient to sample the contribution from a subset of non symmetry equivalent k points only Therefore the in tegrand e g the charge density is calculated only at these points The integrand with the full symmetry can be recovered from its representation by non symmetry equivalent k points whenever this is required The fhi98start utility is set up to automatically exploit these point group symmetries First the point group symmetry operations applicable to the unit cell are determined and stored in the form of symmetry matrices Sec ondly fhi98start seeks to reduce the elements of the k point mesh to the q u subset which is irreducible under those symmetry matrices Only this subset is forwarded in the inp ini file for further use in the main computations The performance of the reduction procedure can be monitored by inspecting the output in the file start out An estimate for the sampling quality of the k point set is given on the basis of the analysis of shells see Chadi and Cohen 31 for details For a good k point set the contribution from the leading shells should
57. or nrx 2 and nrz 3 Using smaller values for nre means skipping the highest G vectors in n r Se D Daa wi fix aatake A anand results in a faster performance However the applicability of the grid should then be carefully checked For systems with strongly localized or bitals in particular this may be an unacceptable approximation set up of the initial wavefunction the initial wave function is obtained by diagonal ization of the Hamiltonian matrix in the basis set as specified by tinit_basis the initial wavefunction is read in from file fort 70 connected to tinit_basis Type integer maximum number of electronic steps if tfor and tdyn are false or maximum number of atomic moves if tfor and tdyn are true integer maximum number of steps in the initial ization if nbeg is 1 see also nbeg and init_basis 2 2 The input files inp mod and start inp 10 Parameter Type Value nOrder integer order of the scheme for solving the ionic equation of motion if idyn is set to 0 or 1 predictor corrector 0 1 2 il 2 3 2 3 4 3 4 5 connected to idyn Parameter Type nstepe integer if tfor or tdyn is true maximum num ber of electronic iterations allowed to con verge forces the program terminates after nstepe iterations see also force_eps Parameter Type pfft store real fraction of wavefunctions for which a second transformation to real space is avoided currently not implemented Parameter Type Value t
58. postc logical post LDA functional true post LDA with functional i_xc false start with functional i_xc Parameter Type tdipol logical if set to true the surface dipol correc tion is calculated Parameter Type tdyn logical if set to true a molecular dynamics simulation is performed tfor must be set to false connected to force_eps and tford 2 2 The input files inp mod and start inp 11 Parameter Type tfor logical if set to true ionic positions are re laxed tdyn must be set to false connected to force_eps epsfor and tford Parameter Type timequeue integer maximum CPU time in seconds the pro gram writes output and restart files before the limit is exceeded and the program is terminated automatically Parameter Type trane logical if set to true the initial wave functions are perturbed by the value trane connected to ampre Parameter Value tranp logical if set to true the ions which are allowed to move are perturbed by amprp Bohr connected to amprp Parameter Type Value tsdp logical scheme for structural optimization true modified steepest descent scheme false damped dynamics scheme L Example The sample input file inp mod below sets up a typical bulk cal culation for metallic Ga The parameters are given in the order required by the fhi98md program The corresponding start inp file is given in the section below 1 100 1000000 100 2 6 0 0 2 4 0 0 2 0 0001 400 2 0 0 1 0 22 0 fal
59. rce_eps 2 Parameter gamma gamma Type real real real criteria to end self consistent cycle in a u stop if the average change of wave functions for the last three iterations is less than epsekinc and if the variation of the total energy is less than epsel for the last three iterations and if the forces on ions are smaller than epsfor this parameter is only active if tfor and tfordis true Type real Type real real Value real real if the total energy varies less than eps chg_dit the parameters delt and gamma are replaced by delt2 and gamma2 connected to parameter delt delt2 gamma gamma convergence criteria for local and total forces maximum allowed relative variation in lo cal forces before if tfor is true exe cuting a geometry optimization step or if tdyn is true calculating total forces maximum allowed relative variation in to tal forces before moving ions if tdyn is true connected to parameter tfor and tdyn damping parameter for the second or der electronic minimization scheme only used if i_edyn is 2 second damping parameter refer to eps_chg_delt connected to i_edyn and eps_chg_delt 2 2 The input files inp mod and start inp Parameter i_edyn Parameter 1_X C Parameter idyn Parameter initbasis Parameter iprint Value Value gt e N Value Value scheme to iterate the wave functions steepest descent
60. rho 5 2 1234 npos nthm nseed 873 00 1400 00 0 1000E 09 1 T_ion T_init Q nfi_rescale 1 T T nsp tpsmesh coordwave 39 nkpt 0 1000000 0 0000000 0 1182033 0 0160000 xk 1 3 wkpt 8 25000000 0 00000000 0 00000000 lattice vector al 0 00000000 4 11262491 6 97950013 lattice vector a2 0 00000000 4 11262491 6 97950013 lattice vector a3 1 00000000 0 00000000 0 00000000 rec lattice vector bi 0 00000000 1 00300905 0 59101654 rec lattice vector b2 0 00000000 1 00300905 0 59101654 rec lattice vector b3 8 2500000 473 61709094 alat omega gt Gallium 2 4 3 00000 69 72000 name number valence charge ion_fac 0 70000 1 00000 3 3 lon_damp rgauss 1_max 1_loc T T F t_init_basis s p d 0 651750030 0 000000000 2 135727000 F F F T 0 651750030 0 000000000 2 135727000 F F F T 4 776750030 0 000000000 4 843773130 F F F T 3 473249970 0 000000000 9 115227130 F F F T 0 0 0 ineq_pos 4 nrot number of symmetries 1 1 0 0 0 1 0 o o0 1 2 1 0 0 o o0 1 0 1 0 3 1 O 0 oO i o 1 O 4 1 0 0 0 1 Q 0 0 1 0 2 4 Pseudopotential file s 23 2 4 Pseudopotential file s Besides the inp mod and start inp inp ini files considered so far the fhi98md re quires supply of pseudopotentials for each of the nsp atomic species listed in start inp Pseudopotential file s must be provided in the working directory ac cording to the following naming convention species pseudopotential filename 1 gt
61. rting with a trial wave function yd The energy minimization scheme is formulated in terms of an equation of motion for the wave function wi in a fictitious time variable t e Steepest Descent The simplest scheme to iterate the wave functions is the steepest descent approach 27 It can be derived from a first order equation of motion d Ca Wie Ex fiks Ww imposing the ortho normality constraint pf As 2 i j where fks is the Kohn Sham Hamiltonian and k are the Lagrange multipliers introduced to account for the ortho normality constraint In the simplest possible discretization of this differential equation only information from the last step is used G klyfk G klou 8 G kW n G kl Hes with 8 i k t and 7 t However it turns out that this discretization scheme is not very efficient e Damped Joannopoulos A more efficient scheme based on a second order equation of motion might also be used Z yi 27 S169 a i WS e where y is a damping parameter The equation of motion is integrated for a step length t by the Joannopoulos approach 28 which iteratively improves the initial wave functions In this algorithm the new wave function roe Dy is constructed from the wave functions of the last two iteration steps t and t 1 G kyt G kyt be G kly 7a G k na G k xslyt 3 3 Total Energy Minimization Schemes 39 where the coefficients are n he d
62. se F 0 001 F 0 002 false false false 1800 false 0 0001 0 0005 0 2 0 001 0 001 3 3 nbeg inp mod iprint timequeue nomore nomore_init delt gamma delt2 gamma2 eps_chg_dlt delt_ion nOrder pfft_store mesh_accuracy idyn i_edyn i_xc t_postc trane ampre tranp amprp tfor tdyn tsdp nstepe tdipol epsel epsfor epsekinc force_eps 1 force_eps 2 max_no_force init_basis 2 2 The input files inp mod and start inp 12 CI start inp The file start inp contains all structural informations the geometry of the supercell and e g the positions of the nuclei All parameters in the start inp necessary for the fhi98md program are given in alphabetical order below Parameter Type atom character 10 name of the pseudopotential supplied necessary for each species Parameter Type celldm 1 6 real lattice parameters of the supercell celldm 1 gives typically the lattice con stant in bohr Parameter Type coordwave logical if set to true and nrho is set to 2 the ionic positions are read in from the file fort 70 connected to nrho Parameter Type ecut real plane wave energy cutoff in Rydberg the cutoff depends on the pseudopotentials and has to be individually checked for ev ery system Parameter Value ecutt real plane wave energy cutoff for the initializa tion if plane waves or a mixed basis are chosen by the parameter init_basis in the file inp mod connected to init_basis
63. t 1 G k Axs G k e771 Ba 7 ik G k Hxs G k ya e i he dt e77 1 G ik G k Hxs G k with W Aks v The function h t is defined by 2e cos wa dt if wz gt 0 he st 2e cosh Veal t if w2 lt 0 with w G k Aks G k gx Williams Soler Although the damped Joannopoulos algorithm is more efficient than the first order scheme additional storage for the wave function oe Dy is needed Therefore the Williams Soler algorithm 29 is recommended whenever storage requirements do not permit to employ the damped Joannopoulos algorithm The coefficients of this scheme are oan jeer alee i ik G k Hxs G k o in G k Axs G k dt _ 1 NG See a ao aoa ee ik G k Hxs G k with yg 0 Thus the damped Joannopoulos scheme contains the Williams Soler scheme as a limiting case when y gt oo On the other hand the Williams Soler scheme itself approaches the steepest descent scheme if t is sufficiently small The choice of dt and y depends on the atomic species and the configu ration The corresponding parameters in the inp mod input file are delt and gamma Typically delt lies between 1 and 40 and gamma is within the rangeO lt y7 lt 1 Set up the parameters delt gamma delt2 and gammaz2 in the file inp mod If the improvement in the total energy per iteration is less than eps_chg _delt the parameters delt2 and gamma2 are used
64. t from Monkhorst Pack are systems with hexagonal symmetry e g slabs with 111 surface of fec metals Here Cunningham proposes to use a hexagonal k point 3 2 Choice of the k point mesh 36 mesh To realize such meshes in the fhi98md code one has to provide ex plicitly a list of k points forming the desired pattern Cunningham s 6 point pattern in the full Brillouin zone can be obtained as follows Parameter nkpt xk 1 xk 1 xk 1 uk 1 uk 1 uk 1 3 wkpt 3 wkpt 3 wkpt 8 wkpt 8 wkpt 3 wkpt i_facs 1 3 t_kpoint_rel Value 6 0 33333 0 00000 0 0 0 0 66667 0 00000 0 0 0 0 00000 0 33333 0 0 0 0 00000 0 66667 0 0 0 0 66667 0 33333 0 0 0 0 33333 0 66667 0 0 0 111 true 16667 16667 16667 16667 16667 16667 number of k points k points and weights k point folding factors frame of reference for k points Figure 3 2 2D Brillouin zone of a fec 111 surface with hexagonal symmetry with set of 6 special k points following Cunningham The thin polygon indicates the conventional first Brillouin zone the thick polygon marks the Brillouin zone as realized in the fhi98md code When a denser mesh in the same cell is desired Cunningham s 18 point pat tern is obtained from the input Parameter nkpt xk 1 xk 1 3 wkpt 3 wkpt i_fa 1 3 t_kpoint_rel Value 2 0 3333 0 3333 0 0 0 5 0 6666 0 6666 0 0 0 5 3 3 il true number of k points supplie
65. tal tdegen true false 1 tmold tband nrho 5 2 1234 npos nthm nseed 873 0 1400 0 1e8 1 T_ion T_init Q nfi_rescale t true tpsmesh coordwave 2 5 Arsenic number of atoms zv name 1 0 74 92 0 7 3 3 gauss radius mass damping 1_max 1_loc ts oan Sake t_init_basis 0 0 0 0 ifs tau0 tford 0 0 0 0 ats tau0 tford ood 0 0 31 start inp Note that in this example the coordinates of the two nuclei in the basis parameter taul are fictitious parameters latgen calls a special rou tine that uses only information provided in celldm parameter to generate the atomic positions Here is the protocol from the fhi98start run saved in the file start out FARK RE KK fhi98md start utility FR kk kk k Akk kk kk January 1999 FR kk k k k number of species 1 number of excess electrons 0 number of empty states 5 ibrav pgind 10 0 celldm 7 10000 2 67000 22700 1 00000 1 00000 1 00000 number of k points 1 k_point 25000 25000 25000 1 00000 i_fac 5 5 5 t_kpoint_rel T ecut ecuti 10 0000 4 0000 ekt tmetal tdegen 1000 T F tmold tband nrho TF 1 npos nthm nseed 5 2 1234 T_ion T_init Q nfi_rescale 873 0001400 000100000000 00 ul tpsmesh coordwave TT 25 00 Arsenic 1 00 74 92000 70000 3 3 TTF 00000 00000 00000 F 00000 00000 00000 F gt latgen anx any anz 1 1 1 this is atpos_special positions tau0 from uniti0 and atpos tau0 species Nr x y z Ars
66. tal temperature set by T_init 4 set according to 1 but the total momen tum is set to zero 5 set according to 3 but the total momen tum is set to zero 6 restart from file fort 20 connected to tdyn T_init tau0 and vau0 Parameter Type Value nrho integer set up of the initial electron density 1 superposition of atomic electron densities 2 constructed from file fort 70 3 read in from file fort 72 Parameter Type only for MD simulations if tdyn is true nseed integer seed used to generate initial velocities Parameter Type nsp integer number of atomic species 2 2 The input files inp mod and start inp 16 Parameter Type Value only for MD simulations if tdyn is true nthm integer simulation ensemble for the ions only ac tive if tdyn is set to true 0 micro canonical ensemble 1 micro canonical ensemble but the veloci ties are rescaled 2 canonical ensemble Nose Hoover connected to nfi_rescale and Q Parameter Type Value pgind integer point group index 0 automatic symmetries and center 1 no symmetries assumed see also parameter ibrav and latgen f connected to ibrav see also How to set up atomic geometries Parameter Type only for MD simulations if tdyn is true Q real gt 0 if nthm 2 mass of thermostat in a u Parameter Type t_init_basis 3 x logical if set to true include s p and d LCAO orbitals in the initialization only active if init_basis is 2 or 3 connected to init_basis P
67. tep Dvolt Wa 0000 0000 3 6 How to set up a band structure calculation 56 Energy eV X WwW K Tr Figure 3 8 GaAs bulk band structure calculated with fhi98md within the local density approximation LDA The Kohn Sham eigenvalues extracted from the above output file fort 6 are denoted by o and the Brillouin zone of the fcc lattice is shown in the inset The energy is measured with respect to the Erfermi dashed line obtained from the self consistent run used to produce the fort 72 file see Pre settings section Note that in case of systems with gap in the energy spectrum Efermi in a band structure run is only assured to be between the valence band maximum and the conduction band minimum thus having no particular physical meaning Bibliography 1 P Hohenberg and W Kohn Phys Rev 136B 864 1964 2 W Andreoni F Gygi and M Parrinello Phys Rev Lett 68 823 1992 3 N Troullier and J L Martins Phys Rev B 46 1754 1992 4 G Ortiz Phys Rev B 45 11328 1992 5 A Garcia et al Phys Rev B 46 9829 1992 6 S B Zhang and J E Northrup Phys Rev Lett 67 2339 1991 7 M Fuchs M Bockstedte E Pehlke and M Scheffler Phys Rev B 57 2134 1998 8 J Neugebauer and M Scheffler Phys Rev B 46 16067 1992 9 R Stumpf and M Scheffler Phys Rev B 53 4958 1996 10 C Stampfl and M Scheffler Phys Rev Lett 78 15
68. this routine for nsp gt 6 2 4 Pseudopotential file s 24 1 yes ionic pseudopotential in l hartree Goo aal field description gen number of valence electrons number of pseudopot com Imax 1 ponents number of radial mesh block for 1 0 Mmax points 13 i Qmesh mesh increment rm41 Tm l l l l m mesh index l l ene Tm radial coordinate in bohr gt a ae ete N f l up radial pseudo wavefunction l l l l l l Mmax 11 I i appears only if core valence l exchange correlation applies l see Ref 26 l I l Figure 2 2 Format of the pseudopotential file name cpi as generated by the psgen tool When setting up the start inp file one should ensure that the values of the zv and lmaz parameters match those of the Zion and lmax 1 fields respectively in the pseudopotential files for each species The fhi98md program stops if zv Zion or if not enough angular momentum components have been provided i e lmax 1 lt Imax Parameter Type Imax integer highest angular momentum of the pseu dopotential 1 gt s 2 gt p 3 gt d lloc integer angular momentum of the local pseudo potential lloc lt lmaz Special care should be paid also to the lloc parameter s in file start inp it must be set to the same value used in generating the pseudopotential If for any Ni reason the user would like to change lloc then a new pseudopotential has to be constructed according to the new lloc value Some
69. times an explicit account of the core valence nonlinearity of the exchange correlation functional is required for instance in studies involving alkali metal atoms 26 In this case the psgen tool appends at the end of the name cpi file a data block containing the partial core density Fig 2 2 The fhi98md program automatically recognizes the use of such a pseudopotential and the proper informa tion record is written to fort 6 during the initialization Note also that pseudopotentials should be generated and used within the same exchange correlation scheme 2 5 Input files for advanced users 25 2 5 Input files for advanced users 2 5 1 The input file constraints ini This file allows the user to specify constraints for the atomic motion when a struc tural relaxation run is performed For using this feature consult the section 3 4 If no constraints are required this files contains two lines both with a O zero in it 2 5 2 The input file inp occ This input file allows for a calculation where the occupation of the eigenstates of the system is user specified This feature is mostly used in conjunction with calculations of atoms and molecules To activate this feature set tdegen to true in the input file start inp The file inp occ should consist of one line with the occupation numbers real numbers between 0 0 and 2 0 listed in the order of ascending energy eigenvalues A typical file that would impose spherical symmetry on a three
70. ulk band structure Use the same bulk unit cell as the one for which the band structure will be calculated Usually this is the minimum cell e g two atoms for Si or GaAs e Projected bulk band structure To make the bulk unit cell first take the same lateral unit cell as for the surface then take the minimum periodicity in the third direction This is required in order to have correspondence between the kj coordinates for the bulk and the surface respectively e Surface band structure Use the same surface unit cell as the one for which the band structure will be calculated 1 Band structure run e Switch on the band structure run mode of fhi98md by setting the in dexemtband parameter in start inp to true Parameter Value tband true the k point set is not reduced by fhi98start utility the electron density is not recalculated after the initialization e Set the value of nrho parameter in start inp to 3 Parameter Value nrho 3 the initial electron density is read in from fort 72 3 6 How to set up a band structure calculation 50 e Sample the symmetry lines in the Brillouin zone along which you want to calculate the energy bands the resulting k point set should be declared in start inp in the standard format see also Choice of the k point mesh Be sure to set the proper value for the tkpoint_rel parameter in the same file depending on the reference frame used to determine the k point coordinates ak 1 3 parameters
71. valent single atom e g Ga In etc regardless of the symmetry of the unit cell employed in the calculation is shown below inp occ 2 0 0 333333 0 333333 0 333333 0 0 Thus the s subshell is completely occupied by 2 electrons and to each of the three p orbitals pz is assigned 1 3 occupancy Note however that in this particular case the nempty variable in start inp should be set to 2 and not to 1 as one could deduce from inp occ This is easy to realize having in mind that the maximum number nz of electronic states per k point is defined as nz nempty nel 1 0 2 0 2 6 Runtime control files The files stopprogram and stopfile allow to control the program execution during runtime They both consist of one line with one integer number The file stop program can be used to stop the program execution deliberately If the number of electronic iterations already performed exceeds the number given in stopprogram the fhi98md program terminates The file stopfile allows to reset the variable nomore while the program is run ning Please notice that the meaning of nomore depends on the mode in which the program is run electronic structure only or relaxation molecular dynamics run In both cases the program only ends after all output files are written thus enabling a continuation of the run see section How to set up a continuation run 2 7 The output files 26 2 7 The output files 1 Main output files filename ener
72. vided in celldm 3 In certain cases ibrav 4 10 12 the input format allows also for ny X Ng X ng scaling of the super cell The three scaling factors ni i 1 2 3 are specified in celldm 4 6 respectively celldm i ibrav 1 2 3 4 5 6 iL Ao a 4 a c a 2 n n2 n3 8 a b c S 10 a c a u n no n3 12 a a ni n2 n3 Specify coordinates of the nuclei parameter tau0 1 3 in start inp How to set up atomic geometries 30 The value of ibrav determines the units in which tau0 1 3 are to be given ibrav units taw0 1 3 1 2 aiat celldm 1 3 4 8 10 12 atomic units ag 0 529177 A e For ibrav 4 10 and 12 positions of the nuclei are solely determined S by celldm therefore the supplied values of tau0 1 3 are not significant Thus the following fragment from start inp is an allowed input 10 0 ibrav pgind 7 1 2 67 0 227 1 1 1 celldm 2 5 Arsenic number of atoms zv name 1 0 74 92 0 7 3 3 gauss radius mass damping 1_max 1_loc eGo Can of t_init_basis 0 0 0 0 0 0 f taud tford 0 0 0 0 0 0 f taud tford _J Other features e Cluster calculations To set up a cluster calculation with fhi98md you need to follow the same steps described above It is important however that you take a sufficiently large super cell in order to avoid coupling between the periodic images of the finite system The common practice is to place the cluster in a cubic ibrav 1 or orthorhombic ibrav

FHI98MD Computer code for density-functional theory calculations

Contents

Download Pdf Manuals

Related Search

Related Contents