Basic tutorial to CPMD calculations
Contents
1. ECLASSIC and EHAM during a CPMD calculation The effect of the modification of u and or Atmax is shown on the evolution of ECLASSIC and EKINC 24 Parameters of the Nos thermostat For details see references 10 13 4 1 Nos thermostats in the CPMD method When introducing Nos thermostats on ions and or fictitious electrons the CPMD equations of motion Eq 3 2 and 3 3 are modified and friction terms are introduced to couple atoms and or wave function motions to the thermostats ER OR uv ATEO xvi i PHP HM constraints M R xp 4 1 MR RO OR constraints u Xf Xe 4 2 where the last term of each equation MjRjig and LV jXe are friction terms that couple respectively atoms and wave functions dynamics to the thermostats These frictions terms are governed by the variable xg and x which obey the following equations of motion 1 1 Orig 2 Y MR 5 8 hy T 4 3 QX 2 Y n Wilh Erne 4 4 where 1 gKkpT is the average kinetic energy of the ionic sub system g 3N is the number of degrees of freedom for the atomic motion in a system with N atoms k is the Boltzmann constant and T the physical temperature of the system And Ekin e is the average kinetic energy of the fictitious electronic sub system The masses Or and O determine the time scale of the thermal fluctuations of the thermostats 25 Chapter 4 Parameters of the Nos thermostat 26 42. Negative numbers all atoms to be used 1 2 7 2 Impose constraints on atoms Blue moon or fix atom positions amp ATOMS CONSTRAINTS FIX STRUCTURE il DIST 1 2 7 200 FIX ATOMS 2 10 11 amp END With FIX STRUCTURE we fix the distance IR 1 R 2 I to 7 2 a u With FIX ATOMS we fix some atoms here 2 atoms number 10 and 11 and no force will act on these atoms Remarks The section DFT describe the XC functional other sections like VDW and ATOMS require keyword s parameters sensitive to the XC functional chosen by the user in the DFT section Consistency in the choice of the XC functional in each of these sections is an absolute pre requisite to any CPMD calculation CPMD calculations 2 1 Wavefunction optimization Before any molecular dynamics calculation it is mandatory to optimize the wavefunction of the system Since in the CPMD method the time independent Schrodinger equation is not solved at each step it is mandatory to start the calculation with the best possible electronic structure and thus the first step of any CPMD calculation is to optimize the wavefunction In the CPMD section add the keywords amp CPMD OPTIMIZE WAVEFUNCTION INITIALIZE WAVEFUNCTION RANDOM TIMESTEP 340 MAXITER 50000 PCG CONVERGENCE ORBITALS 8 0H 06 EMASS 2000 amp END e PCG Use the preconditioned conjugate gradients for optimization other are available DIIS SD e CONVE
3. 13 12210 13 12210 0 322E 08 0 96 10 0 00000 0 0 13 12210 13 12210 13 12210 0 488E 08 0 96 e NFI step number e EKINC kinetic energy of the fictitious electronic sub system a u e TEMPP temperature of the system K e EKS Kohn Sham energy a u like potential energy in classical MD e ECLASSIC EKS Ionic Kinetic energy a u e EHAM EHAM ECLASSIC EKINC the conserved parameter in CPMD e DIS Mean square displacement of atoms since the first MD step a u e TCPU CPU time for the corresponding MD step Table 2 2 Typical output of a Car Parrinello molecular dynamics run in the CPMD code 12 13 2 4 Temperature control CPMD dynamics Remarks 1 When changing the calculation from OPTIMIZE GEOMETRY to MOLECULAR DYNAMICS one must not use the velocities from the restart file 2 In the output file labelled ENERGIES the first column presents the kinetic energy of the fictitious electronic part If the values in this column are too important and or if they are subject to too much variation then there is a problem in the calculation 2 4 Temperature control CPMD dynamics In the CPMD section add the keywords amp CPMD MOLECULAR DYNAMICS TIMESTEP 3 MAXSTEP 50000 TEMPCONTROL IONS 800 100 amp END The first number is the target temperature and the second is the tolerance around this value Both are given in K In this example the temperature of the system is to be oscillated around 80
4. ES and the ACCUMULATORS ex time step number derivatives of many variables The keyword ACCUMULATORS is interesting during a molecular dynamics if the details of the calculation remain identical between each run If the temperature the pressure or else changes then the keyword ACCUMULATORS must not be used In this case the RESTART files will be written alternatively every 50 steps in this example 2 restart files will be written RESTART 1 and RESTART 2 it can be more than 2 files RESTART 1 will be written at 50 then RESTART 2 will be written at 100 then RESTART 1 will be re written at 150 and so on 19 Controlling adiabaticity in CPMD For details see references 1 2 9 At this point a short introduction on the Car and Parrinello method is required 3 1 Car Parrinello equations of motion Car and Parrinello postulated the following class of Lagrangians 1 1 Lcp IR UL UIT V 71 constraints 3 1 2 2 L A Sie en ec dt Potential energy o ensure orthonorma ity Kinetic energies The generic Car Parrinello equations of motion are of the form MiR ft am VIA ag constraints 3 2 d d w t ra a 3 3 Equations of motion Eq 3 2 and 3 3 refer respectively to the ionic and the fictitious electronic sub systems u u are the fictitious masses assigned to the electronic orbitals degrees of freedom According to the Car Parrinello equations of motio
5. 0 100 K if the temperature goes bellow 700 K or higher than 900 K then it is readjusted Optional keywords for the fictitious electronic sub system are amp CPMD TEMPCONTROL ELECTRONS 0 04 0 01 amp END The first number is the target fictitious kinetic energy and the second is the tolerance around this value Both are given in atomic units 13 Chapter 2 CPMD calculations 14 2 5 Nos thermostat In the CPMD section add the keywords amp CPMD MOLECULAR DYNAMICS TIMESTEP Jol MAXSTEP 50000 NOSE IONS 1000 200 NOSE ELECTRONS 0 2 600 amp END In the first case IONS the first number represents the target thermodynamic temperature in K the second the thermostat frequency in cm In the second case ELECTRONS the first number represents the target fictitious kinetic energy the second the thermostat frequency in cm The value of the target fictitious kinetic energy has to be of the order of the value of the fictitious kinetic energy observed in a free dynamics simulation little bit higher is good When the Nos thermostat is used the frequency must be much higher for the ELECTRONS than for the IONS The two Nos thermostats IONS and ELECTRONS must played on distinct frequencies yet leading to a similar thermalization for the two sub systems see chapter Chap 4 for details 2 6 Electronic density of states Note that several methods exist in order to achieve this calculation two of these methods
6. 2 Choosing the parameters In a CPMD run controlled using Nos thermostat s the physical temperature T is the most basic requirement for the thermostat that controls the ionic degrees of freedom Then it is necessary to input the mass Or associated to this thermostat Similarly when using a thermostat to control the fictitious electronic degrees of freedom it is required to input Ej and Qe 421 Exine A possibility to determine a value for Egin is to estimate the kinetic energy Exin ad related to the adiabatic motion of a model system of well known separated atoms The adiabatic atomic wave functions of this model follow the corresponding atom rigidly thus it is possible to relate the kinetic energy to the velocities of the atoms U 1 Exinad 2kyT R lt V O 4 5 i M is the mass of one atom The value of Exin aa gives a guideline for setting the average kinetic energy of the electronic wave functions since appropriate values of Ekin have to be larger than Ein a4 For Exin e the recommended value is about twice Exin ad 13 42 20 Qpr and 0 The frequencies Org of the thermal fluctuations associated to the thermostat that controls the ionic degrees of freedom are determined through 2gk T OR OTR 4 6 And for the frequencies 7e of the thermal fluctuations associated to the thermostat that con trols fictitious electronic degrees of freedom through 4 Exin Qe OTe 4 7 Basically the i
7. ANNIER amp END This section is used to describe how to treat Van Der Waals interactions amp VDW WANNIER CORRECTION VERSION 2 FRAGMENT BOND 15 96 TOLERANCE WANNIER do TOLERANCE REFERENCE 1 0 PRINT INFO FRAGMENT END WANNIER CORRECTION amp END 1 2 5 The PROP section If the calculation described in the CPMD section of the input file is PROPERTIES amp CPMD PROPERTIES WANNIER amp END This section is used to provide details about the physical properties to be calculated Chapter 1 CPMD basics 6 amp PROP LOCALIZE PROJECT WAVEFUNCTION POPULATION ANALYSIS MULLIKEN amp END 1 2 6 The SYSTEM section This section is used to describe the geometry of the system to be studied SYSTEM ANGSTROM SYMMETRY il CELL ABSOLUTE DEGREE 10 00 10 00 10 00 90 00 390 00 20 00 CUTOFF 100 0 KPOINTS MONKHORST PACK 2 2 i amp END In this example e The ANGSTROM keyword specify that cell parameters A B and C as well as atomic coor dinates see section Sec 1 2 7 will be provided in Angstrom In the CPMD code the default unit for length and positions is atomic unit e The SYMMETRY keyword number on the following line between 0 and 14 0 isolated system molecule 1 cubic FCC describes the symmetry e The CELL keyword specify that cells parameters and angles will be read on the next line 6 values e The ABSOLUTE keyword specify that cell parameters will be read in the
8. BE o sG 4 a EROR we a erde Electron localization function ELF 22 2 CC mn 2 8 1 CPMDkeywords lt se sa saose em ema rs 28 2 CBMD2CUBE e seti eck ae ed un a en RUE oe AE A eg Maximally localized Wannier functions a 2 9 1 CPMD keywords iii ii 292 EBMD2CUBE 3 aid ct e ES ee ud a a ege 18 2 10 Restarts 5 coco u ua a e aa ew a BG E A ee Ve 19 Controlling adiabaticity in CPMD 21 3 1 Car Parrinello equations of motion ee 21 3 2 Howto control adiabaticity 2 2 22 2 Como 22 Parameters of the Nos thermostat 25 4 1 Nos thermostats in the CPMD method 2 02 200004 25 4 2 Choosing the parameters 2 rs 26 ALL Feme wen Sw ual a Bow hah DEG a ho ar qw Rege RICE ase 26 4 2 2 Or and Qo nce ox excog Ree m eee UROR aa E Oe ee 26 List of Figures 3 1 Schematic representation of the vibration frequencies of the ionic and the electronic sub system in CPMD Electronic frequencies have to be higher than the ionic frequencies and the two must not OVerldpi s hose Baer Oe So REN UE ea ee Se as Swe we ooh ae 22 3 2 Schematic representation of the evolution of EKINC EKS ECLASSIC and EHAM during a CPMD calculation The effect of the modification of u and or Atmax is shown on the evolution Of ECEASSIC and EK TN Gyere EA A e ee wk Oe dp doe ecw ER ERR de 24 4 1 Schematic representation of the thermal frequencies of the ionic and the electronic sub system as well as the Nos
9. Car and Parinnello Molecular Dynamics http www cpmd org Basic tutorial to CPMD calculations S bastien LE ROUX sebastien leroux ipcms unistra fr INSTITUT DE PHYSIQUE ET DE CHIMIE DES MATERIAUX DE STRASBOURG DEPARTEMENT DES MATERIAUX ORGANIQUES 23 RUE DU LOESS BP43 F 67034 STRASBOURG CEDEX 2 FRANCE Contents Contents List of figures List of tables 1 CPMD basics General ideas on the CPMD method 2 KK m m m nn 1 2 Structure of the CPMD input file oo Con 1 1 1 2 1 The INFO section 1 2 2 The CPMD section 1 2 3 The DFT section 1 2 4 The VDW section 1 2 5 The PROP section 12 6 The SYSTEM sections a s i R a eso a e di eee RO doc RW amp N 1 2 7 The ATOMS section 2 CPMD calculations 2 1 2 2 2 3 2 4 2 5 2 6 2 7 2 8 2 9 Wavefunction optimization 2 2 lees Geometry optimization 2 4 42e rs 2 2 1 Using the OPTIMIZE GEOMETRY keywords the automatic way 2 2 2 Using Molecular Dynamics Friction forces the manual way Free CPMD dynamics 2 rs Temperature control CPMD dynamics 0 00002 ee eee Nos thermostat duck wa u a an a a a e a ROSE ee el ge eo Electronic density of states ee 2 6 1 Kohn Sham energies calculation 2 lens 2 6 2 Free energy calculation 22e 3D visualization of orbitals electronic densities 22 2 2 CC mo mn 2 7 1 CPMDkeywords ae scs a cosac ea rs 2 13 CPMD2CU
10. PMD calculations 18 2 9 Maximally localized Wannier functions 2 9 1 CPMD keywords In the CPMD section the calculation to be run must be changed to PROPERTIES followed by the WANNIER keyword and some additional keywords to describe the details of the analysis amp CPMD PROPERTIES WANNIER WANNIER WFNOUT ALL amp END Then create a new properties PROP section that contains the following amp PROP LOCALIZE amp END This calculation will output in particular two files e WANNIER_CENTER that contains the positions of the Wannier centers e IONS CENTERS xyz that contains the positions of both the Wannier centers and the atoms 2 9 2 CPMD2CUBE The WANNIER_ file is a density file similar to the DENSITY files and has to be converted using the CPMD2CUBE code however the dens keyword is not to be used cpmd2cube x centre WANNIER_1 1 18 19 2 10 Restarts In the CPMD section add the keywords amp CPMD 2 10 Restarts RESTART WAVEFUNCTION COORDINATES VELOCITIES ACCUMULATORS LATEST STORE 50 RESTFILE 2 amp END Optional keywords if Nos thermostats are used amp CPMD RESTART NOSEE NOSEP LATEST amp END With these keywords the calculation will be resumed with the information contained in a CPMD restart file which names is specified in the file LATEST using the data regarding the WAVEFUNCTION the atomic COORDINATES the V ELOCITI
11. RG ENCE ORBITALS Convergence criteria maximum value for the biggest element of the gradient of the wavefunction default value 10 a u Chapter 2 CPMD calculations 10 e MAXITER Maximum number of optimization steps Default 10000 in the case of wave function optimization the keyword MAXSTEP which is a parameter required for geometry optimization or molecular dynamics can also be used e WAVEFUNCTION can also be initialized from the atom positions then instead of RANDOM use ATOMS as keyword e EMASS Fiction electron mass in a u default value 400 a u 2 2 Geometry optimization 2 2 1 Using the OPTIMIZE GEOMETRY keywords the automatic way In the CPMD section add the keywords amp CPMD OPTIMIZE GEOMETRY HESSIAN UNIT MAXSTEP 50000 CONVERGENCE GEOMETRY 1 0704 EMASS 600 amp END The HESSIAN UNITS keyword is required for the first run if more than one are required this is the intial approximate Hessian for geometry optimization The optimize geometry routines in CPMD are automatic the calculation will end when the atomic forces get lower than 1074 a u default is 5x107 a u in other words when the energy will reach a minimum value and its derivative will become almost equal to zero It can work well nevertheless if we start the calculation with a geometry too far from the real exact minimum then this minimum may be impossible to reach and the optimization will fail To
12. ac distribution in order to fill the electronic states as a function of the temperature The first part of the calculation is a Born Oppenheimer wave function optimization to calculate the electronic structure of the system The parameters provided above likely depend on the system under study 15 Chapter 2 CPMD calculations 16 2 7 3D visualization of orbitals electronic densities The calculation and the extraction of the orbitals and densities is done using the CPMD code afterwards the creation of the input files for 3D visualization is done using the CPMD2CUBE code 2 7 1 CPMD keywords In the CPMD section add the keywords amp CPMD RHOOUT BANDS 4 805 806 807 808 amp END Using the RHOOUT BANDS keywords a certain number of bands or orbitals will be plotted The number of orbitals is given in the second line and their index starting from 1 is given in the third line e positive index ex 806 the electronic density is written output files are stored in DENSITY index ex DENSITY 806 e negative index ex 806 the wavefunction is written output files are stored in WAVEFUNCTION index ex WAVEFUNCTION 806 2 7 2 CPMD2CUBE The exact command depends on the nature of the object to be visualized e Fora DENSITY cpmd2cube x dens centre DENSITY 806 e For a WAVEFUNCTION cpmd2cube x wave centre WAVEFUNCTION 806 It is also possible to use the inbox keyword to force atoms
13. are introduced hereafter 2 6 1 Kohn Sham energies calculation In the CPMD section add the keywords amp CPMD RESTART WAVEFUNCTION COORDINATES LATEST KOHN SHAM ENERGIES 100 amp END 14 15 2 6 Electronic density of states The number that appears after the keyword KOHN SHAM ENERGIES represents the number of unoccupied electronic states this number depends on the system and has to be large enough so that the energy of the electronic states that surround the HOMO LUMO gap are perfectly converged This calculation can only be performed using a converged electronic structure that is why one has to use a restart file see section Sec 2 10 for detail from a wave function opti mization see section Sec 2 1 for detail Also worth to mention is the fact that this calculation requires significantly more memory than a standard CPMD calculation both wavefunctions and Hamiltonian required to be stored in memory the latest being diagonalized 2 6 2 Free energy calculation In the CPMD section add the keywords amp CPMD OPTIMIZE WAVEFUNCTIONS FREE ENERGY FUNCTIONAL INITIALIZE WAVEFUNCTIONS RANDOM TROTTER FACTOR 0 001 BOGOLIUBOV CORRECTION LANCZOS DIAGONALISATION LANCZOS PARAMETERS n 5 1000 22 1 D 8 0505 1 DL 0 01 1 D 12 0 0025 1 D 15 0 001 AD ANDERSON MIXING 3 Ike 0 adl 0 01 0 1 0 005 0 11 TEMPERATURE ELECTRON 300 0 amp END The free functional calculation uses the Fermi Dir
14. be run amp CPMD OP INITIALIZE WAVEFUNCTION ATOMS TIMIZE WAVEFUNCTION TIMESTEP MAXITER and or MAXRUNTIME Jol 50000 PCG CONVERGENCE ORBITALS 8 0E 06 EMASS amp END 2000 1 2 Structure of the CPMD input file Many different type of calculation can be performed using the CPMD code among others OPTIMIZE WAVEFUNCTION OPTIMIZE GEOMETRY MOLECULAR DYNAMICS MOLECULAR DYNAMICS BO KOHN SHAM ENERGIES PROPERTIES LINEAR RESPONSE EHERENFEST DYNAMICS TDDFT FREE ENERGY MD many of them are briefly introduced in the next chapter see Chap 2 Chapter 1 CPMD basics 4 1 2 3 The DFT section This section is used to describe the parameters of the Density Functional Theory DFT amp DFT FUNCTIONAL BLYP GC CUTOFF 107 amp END And in particular the exchange correlation functional XC to be used in the calculation among LDA BP BLYP PBE revPBE XLYP HCTH OLYP And Hybrid XC functional B3LYP PBEO HSE06 1 2 4 The VDW section If required in the CPMD section of the input file 1 2 4 1 Empirical VDW corrections amp CPMD VDW CORRECTION SEND This section is used to describe how to treat Van Der Waals interactions amp VDW EMPIRICAL CORRECTION VDW PARAMETERS ALL DFT 2 S6GRIMME BLYP END EMPIRICAL CORRECTION amp END 1 2 Structure of the CPMD input file 1 2 4 2 Wannier VDW corrections amp CPMD VDW W
15. ergy minimization and the resulting energy will be too far from the Born Oppenheimer minimum energy surface This will result in more important vibrations around the total energy for the Hamiltonian EHAM which by construction has to be perfectly conserved when running free dynamics 11 Chapter 2 CPMD calculations 12 2 3 Free CPMD dynamics In the CPMD section add the keywords amp CPMD MOLECULAR DYNAMICS TIMESTEP 2 0 MAXSTEP 50000 amp END In free dynamics the term EHAM which describes the total energy from the Hamiltonian has to be perfectly conserved therefore the initial electronic structure must be determined by optimizing the wave functions as shown in section Sec 2 1 otherwise EHAM will not be conserved and the dynamics will be bad See table Tab 2 2 for the typical output of a Car Parrinello molecular dynamics run in the CPMD code NFI EKINC TEMPP EKS ECLASSIC EHAM DIS TCPU 1 0 00000 0 0 13 12210 13 12210 13 12210 0 515E 12 0 99 2 0 00000 0 0 13 12210 13 12210 13 12210 0 822E 11 0 96 3 0 00000 0 0 13 12210 13 12210 13 12210 0 414E 10 0 96 4 0 00000 0 0 13 12210 13 12210 13 12210 0 130E 09 0 96 5 0 00000 0 0 13 12210 13 12210 13 12210 0 316E 09 0 97 6 0 00000 0 0 13 12210 13 12210 13 12210 0 651E 09 0 96 7 0 00000 0 0 13 12210 13 12210 13 12210 0 120E 08 0 96 8 0 00000 0 0 13 12210 13 12210 13 12210 0 203E 08 0 95 9 0 00000 0 0 13 12210
16. format A B and C The default format in CPMD is A B A and C A e The DEGREE keyword specify that cell angles will be provided in degrees for a B and y The default format in CPMD is cos a cos D and cos y e The CUTOFF keyword value on the following line specify the cutoff for the plane waves in Rydberg e The K points if needed only in Monkhorst Pack representation Here k 2 x 1 x 1 7 1 2 Structure of the CPMD input file 1 2 7 The ATOMS section This section is used to give the atomic positions it follows the structure amp ATOMS PP_FILE_NAME pps PP_GENERATION_METHOD LMAX P LOC P N_ATOMS R yay z1 xe z2 X N ATOMS y N ATOMS z N ATOMS amp END The inner part being repeated as many times as the number of chemical species ex amp ATOMS GE MT BLYP pps KLEINMAN BYLANDER LMAX P LOC P 1 4 7367371096 6 6107492731 4 6811927823 SE MT BLYP pps KLEINMAN BYLANDER LMAX P LOC P 4 5139 8999 SAA 8 7467531007 112930278730 1 8904372056 10 8283093176 4 1633837562 5 4988512174 4 9899517884 10 6568153345 0 00539935793 0 1450182640 4 3887246936 amp END This required to sort the atomic coordinates by chemical species Chapter 1 CPMD basics 1 2 7 1 Generate dummy atoms e g the center of mass of the system amp ATOMS DUMMY ATOMS il Wee 2 i 2 amp END Only 1 dummy atom will be used TYPE2 center of mass TYPE1 to TYPE4 available see manual and two atoms number 1 and 2 are used
17. g is the KS eigenvalue of an empty or occupied state u is the fictitious mass parameter of the electronic sub system It is possible to estimate the lowest frequency of the electronic sub system 07 e E a prin Se gap 3 5 u Egap is the energy difference between the lowest unoccupied LUMO and the highest occupied HOMO orbital c increases like the square root of Egap and decrease similarly with u Remark The relations Eq 3 4 and 3 5 have the important consequence that it is not possible to study metals within the CPMD framework since in that case no gap exists between the HOMO and LUMO 22 23 3 2 How to control adiabaticity To ensure the adiabaticity between the two sub systems the difference omin 07 must be important represents the highest ionic phonon frequency Since both Egap and are quantities whose values are dictated by the physics of the system the only parameter to control adiabaticity is the fictitious mass u therefore also called adia baticity parameter Furthermore the highest frequency of the electronic sub system 07 pa also depends on u Eeu Ur o Z 3 6 U where Ecu is the largest kinetic energy of the wave functions Therefore decreasing u not only shifts the fictitious electronic frequencies upwards but also stretches the entire frequency spectrum Finally the molecular dynamics technique itself introduces a limitation in the decrease of u due to
18. n the nuclei evolve in time at a certain instantaneous temperature proportional to Y MR Following the same idea the fictitious electronic sub system evolves in a fictitious temperature proportional to Y u V V The fast fictitious electronic sub system must have a low fictitious electronic temperature cold electrons where as simultaneously the slow ionic sub system nuclei are kept at much higher temperature hot nuclei The fast electronic sub system stays cold for long times but still follows the slow nuclear motion adiabatically Adiabaticity means that the two sub systems are decoupled and that there is no 21 Chapter 3 Controlling adiabaticity in CPMD 22 M f a miu Figure 3 1 Schematic representation of the vibration frequencies of the ionic and the elec tronic sub system in CPMD Electronic frequencies have to be higher than the ionic frequencies and the two must not overlap energy transfer between the two sub systems in time This is possible if the vibration spectra from both sub systems do not overlap in their frequency domain so that energy transfers from the hot nuclei to the cold electrons become impossible on the relevant time scales this is illustrated by a schematic in figure Fig 3 1 3 2 How to control adiabaticity The dynamics of Kohn Sham orbitals can be described as a superposition of oscillations whose frequencies are given by 28 i ey 2 64 U where
19. onic and fake electronic frequencies will decrease when increasing respectively Or and Qe The schematic in figure Fig 4 1 illustrates how the frequencies of the two thermostats should be set up for a CPMD calculation e Frequencies of the thermostat coupled to the electronic sub system must be lower than the frequencies of the electronic sub system e Frequencies of the thermostat coupled to the ionic sub system must be lower than the frequencies of the ionic sub system 26 21 4 2 Choosing the parameters Bi ionic sub system thermostat electronic sub system thermostat Figure 4 1 Schematic representation of the thermal frequencies of the ionic and the elec tronic sub system as well as the Nos thermostats for both sub systems in CPMD Electronic frequencies fictitious electrons thermostat have to be higher than the ionic frequencies ions thermostat and the two sets of frequencies must not overlap e Electronic frequencies fictitious electrons thermostat have to be higher than the ionic frequencies ions thermostat e Two sets of frequencies electronic and ionic must not overlap 27 Bibliography 1 R Car and M Parrinello Phys Rev Lett 55 22 2471 2474 1985 2 D Marx and J Hutter Mod Met Algo O Chem 1 301 449 2000 3 http www cpmd org 4 A D Becke and K E Edgecombe J Chem Phys 92 5397 5403 1990 5 A Savin A D Becke J Flad R Nesper H Preus
20. reach it or at least to control the results obtained with this method it is possible to use the friction dynamics 10 11 2 2 Geometry optimization 2 2 2 Using Molecular Dynamics Friction forces the manual way In the CPMD section add the keywords amp CPMD MOLECULAR DYNAMICS ANNEALING IONS 008 TIMESTEP Jol MAXSTEP 50000 amp END Remarks 1 For both sections Sec 2 2 1 and 2 2 2 we can consider that the optimization is complete if the forces gradients on X Y Z are lt 1074 An example is presented in table Tab 2 1 KKKKKKKKKKKKKKKKKKKKKKKKKKKKK KKK KKK KKK KKKKKKKKKKKKKKKKKKKKKKKKKK FINAL RESULTS i kkxkxkxkxkxkxkxkxkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkxkXx k ATOM COORDINATES GRADIENTS FORCES 8 4600 6 6095 11 2683 1 619E 04 1 945E 04 4 204E 04 9 7962 14 9558 3596 13 8968 1 1 1899 10 5870 6 7304 6 0976 1 8 5220 8 7119 1 0558 6 261E 05 3 067E 04 6 960E 05 0554 1 658E 04 6 944E 05 8 040E 06 50163 2 748E 04 5 825E 05 3 162E 05 1723 9 913E 05 3 407E 04 3 046E 04 0357 2 546E 05 2 293E 04 3 923E 04 Nos np TEXT e C5 O CO Oo KKKKKKKKKKKKKKKKKKKKKKKKK ko kck kk ckckckck ck ck kck kck ck ck ck ck kc k ck k ck kk k kk kk kkkkkkk Table 2 1 Typical output of a geometry optimization run in the CPMD code 2 If the fictitious mass of the electronic sub system is too small or too big then the frequen cies arising from this mass will interfere with the en
21. rt from the ions and the one from the electrons orbitals in other words there should be absolutely no exchange between the two kinetic parts In the program the word ELECTRONS appears often for the keywords concerning the parametrization of the dynamics of the electronic sub system this terminology is not ap propriate since the dynamical objects are actually the orbitals and not the electrons Chapter 1 CPMD basics 2 1 2 Structure of the CPMD input file The CPMD input file is divided in sections each of them contains specific keywords amp SECTION KEYWORD comment s SEND The section named SECTION start by amp SECTION and end by amp END The CPMD code read keywords in capital letters otherwise the text is considered as comments Many sections can be used to create a CPMD input file see the CPMD user manual for details but only 4 sections are mandatory to run a CPMD calculation CPMD DFT SYSTEM and ATOMS Because of their significant interest few additional sections are presented thereafter 1 2 1 The INFO section This section is used to present details about the calculation and the system under study No keyword is required here the information will simply be appended to the output file at runtime amp INFO GeSe9 480 atoms DFT GGA BLYP FPMD or whatever information I want to use to describ what I am doing amp END 1 2 2 The CPMD section This section is used to describe the calculation to
22. s and H G von Schnering Angew Chem Int Ed in English 30 4 409 412 1991 6 A Savin O Jepsen J Flad O K Andersen H Preuss and H G von Schnering Angew Chem Int Ed in English 31 2 187 188 1992 7 A Savin R Nesper S Wengert and T F F ssler Angew Chem Int Ed in English 36 17 1808 1832 1997 8 D Marx and A Savin Angew Chem Int Ed in English 36 19 2077 2080 1997 9 G Pastore E Smargiassi and F Buda Phys Rev A 44 10 6334 6347 1991 10 S Nos J Chem Phys 81 1 130511 1984 11 S Nos Mol Phys 52 2 255 258 1984 12 W G Hoover Phys Rev A 31 3 1695 1697 1985 13 P E Bl chl and M Parrinello Phys Rev B 45 16 9413 9416 1992 This document has been prepared using the Linux operating system and free softwares The text editor gVim The GNU image manipulation program The Gimp The WYSIWYG plotting tool Grace And the document preparation system BIEX 2e
23. the maximum time step Atmax that can be used The time step Atmax is indeed inversely proportional to the highest electronic frequency in the system which happens to be 97 U x max Oe E cut Atmax 9 3 7 Equation Eq 3 7 thus governs the largest possible time step in a CPMD calculation The choice of the adiabaticity parameter is a compromise between equations Eq 3 5 and 3 7 As illustrated in table Tab 2 2 in a CPMD calculation the quantities to follow and con trol are the kinetic energy of the electronic sub system EKINC the potential energy from the Kohn Sham equations EKS the classical first principles energy ECLASSIC and the conserved energetic quantity of the CPMD Hamiltonian EHAM Figure Fig 3 2 illustrates the evolution of these quantities during the dynamics and highlights the effect of the modification of u and or Atmax On the evolution of ECLASSIC and EKINC EHAM ECLASSIC EKINC then if the variation on ECLASSIC are so important that EKINC interferes with EKS then there is no adiabaticity anymore Remark EKINC should be smaller than 20 of the difference EHAM EKS 23 Chapter 3 Controlling adiabaticity in CPMD 24 EKINC soana KP DR ie AE QUA A uam c erase Dn wat d a 3 rao lo eu EHAM ECLASSIC EKS EKIN EKS Time No adiabaticity Figure 3 2 Schematic representation of the evolution of EKINC EKS
24. thermostats for both sub systems in CPMD Electronic frequencies fictitious electrons thermostat have to be higher than the ionic frequencies ions thermostat and the two sets of frequencies must not overlap 0 aaa 27 ill List of Tables 2 1 Typical output of a geometry optimization run in the CPMD code 0 2 2 Typical output of a Car Parrinello molecular dynamics run in the CPMD code CPMD basics 1 1 General ideas on the CPMD method For details see references 1 3 In the CPMD method the time independent Schrodinger equation is not solved at each step the main idea of the method is to transform the quantum adiabatic time scale separation of fast electronic and slow nuclear motion into classical adiabatic energy scale separation for the dynamics To achieve this goal the two component quantum classical problem of First Principles Ab Initio Molecular Dynamics FPMD AIMD is transformed into a two component purely classical problem with 2 separate energy scales The explicit time dependence of the quantum sub system dynamics is also lost in the process Furthermore since the kinetic evolution of the electronic structure is going to be considered apart from the one of the ions it is necessary to give a fictitious mass to the electronic degrees of freedom E Km Also a term is added to the Hamiltonian of the system to ensure the adiabaticity of the relationship between the kinetic pa
25. to be put inside the simulation box when periodic boundary conditions are used 16 17 2 8 Electron localization function ELF 2 8 Electron localization function ELF The formalism of the electron localization function ELF 4 7 provides a deeper insight on the bonding localization More precisely the ELF gives information on the degree of localization of the electronic density When the ELF values are close to the maximum it is possible to identify for a given value centers in space so called attractors which are surrounded by disjoint ELF isosurfaces so called basins defining different charge localization domains Merging of dif ferent basins for slightly lower values of ELF may reveal the existence of localization domains having in common more than one attractor This gives indications on the centers involved in the bonding and on the spatial extension of this latter 8 The calculation of the ELF requires to add a single keyword to the CPMD section of the input file afterwards the creation of the input files for 3D visualization is done using the CPMD2CUBE code 2 8 1 CPMD keywords In the CPMD section add the keywords amp CPMD ELF SEND At CPMD runtime this keyword will generate the creation of a file named ELF in the active directory 2 8 2 CPMD2CUBE The ELF file is a density file similar to the DENSITY files and has to be converted using the CPMD2CUBE code cpmd2cube x dens centre ELF I7 Chapter 2 C
Download Pdf Manuals
Related Search