======== Examples ======== or how-to's =========== .. contents:: .. section-numbering:: ------------- CO molecule ------------- Setting up a CO molecule in a box. This is the minimum script (or close to the mimimum) one can define: CO_in_a_box.py -------------------------- Relaxing the CO molecule -------------------------- This is the second minimum script one could define. This relaxes the CO geometry until the magnitude of the force on each atom is less than 0.05. CO_relaxed_in_a_box_ .. _CO_relaxed_in_a_box: attachment: CO_relaxed_in_a_box.py -------------------- Al Equation of State -------------------- This example uses the `EquationOfState Module`_ and `GeometricTransforms Module`_ to find the ground state lattice constant and bulk modulus of bulk Al in the fcc crystal structure. Al_equation_of_state_ .. _Al_equation_of_state: inline :Al_equation_of_state.py Once the jobs are done, you could run the command line utility ``eos`` to view the results again and save the figure as a png file.:: eos -m -s Al_murn.png bulkAl_*.nc ``eos`` takes many options:: eos [-m -b --bm -a -t -v -p --all -q -s [filename]] [ncfiles]' these flags specify the equation of state to use. Several can be selected at a time -m 'Murnaghan', -b 'Birch', --bm 'BirchMurnaghan', -t 'Taylor', -p 'PourierTarantola', -v 'Vinet', -a 'AntonSchmidt' --all runs all equations of state -s filename saves figure in file, format by extension. png and eps supported. -q quiet, no gnuplot window if no ncfiles are provided, it uses all ncfiles found in the current directory. shell patterns like bulk*.nc are allowed to use only files matching that pattern. Typical results look like:: Murnaghan: PRB 28, 5480 (1983) V0(A^3) B(Gpa) E0 (eV) pressure(GPa) 16.60 78.27 -56.4688 -0.00 chisq = 0.0000 lattice constant = 4.05 angstroms |Almurn| .. |Almurn| image:: Al_murn.png ---------- H/Co(1000) ---------- This example calculates the total energy of a Hydrogen atom on a Co(1000) surface: H_Co_ontop_ .. _H_Co_ontop: attachment: H_Co_ontop.py This is how the surface looks like, using: >>> from ASE.Visualization.RasMol import RasMol >>> RasMol(slab, (2,2,1)) |HCo| .. |HCo| image:: HCo.jpg :width: 300 -------------------------- Using filters and dynamics -------------------------- The CO molecule is relaxed using a constrain filter `Subset`_ and the `MDMin`_ dynamics: filter_ .. _filter: attachment: filter.py A `NetCDF trajectory`_ is added to the dynamics, writing to a file ``CO-trajectory.nc``. The minimization path can be viewed after the run using the tool `plottrajectory`_ from the ASE/Tools directory:: plottrajectory Co-trajectory.py ---------------- h2o dissociation ---------------- This example use the `Subset`_ filter together with the `FixBondLenght`_ filter to find the barrier for dissociating the h2o molecule: h2o-quasi-Newton_ .. _h2o-quasi-Newton:attachment:h2o-quasi-Newton.py The script use the `quasi-Newton`_ algorithm, reusing the Hessian matrix from one bond length to the next. This can speed up the convergence significantly. The plot below shows the energy versus the total number of iterations, ie. for all bond lengths, comparing different minimization algorithms: |h2o-hessian| .. |h2o-hessian| image:: h2o-hessian.png :width: 600 The total number of iteration is reduced from 95 for MDMin, to 45 for the quasi-Newton and finally to 20 for the quasi-Newton, re-using the Hessian from the previous bond length. Using the `plottrajectory`_ tool one can plot the constrained bond length against the total energy: plottrajectory.py -p [[BondLength(0,1),TotalEnergy]] h2o-dissociation.nc This will look like: |plottrajectory| .. |plottrajectory| image:: plottrajectory.gif :width: 600 ------------ VTK plotting ------------ This example illustrates how one can use VTK to make a combined plot of the wavefunction of a CO molecule, together with the molecule itself: plotwave1_ .. _plotwave1: {{attachment:plotwavefunction.py}} Run the example like:: python -i plotwavefunction.py This example also illustrates how a calculation can be restarted from a existing dacapo NetCDF file. Using:: python -i plotwavefunction.py 3: the following plot is produced: |vtkplot| .. |vtkplot| image:: vtk.gif :width: 400 -------------------------- Calculating a band diagram -------------------------- The example shows how the band diagram between the Gamma point and the X point in the Brillouin zone can be calculated for Cu in a fcc lattice. The plot produced by this example is just plotting the individual eigenvalues, not making a real bandstructure plot. In order to do a band structure calculation one can proceed as follows: 1. First do a calculation with a usual k point sampling. 2. From this calculation get the self-consistent electron density and use this density as input to the band diagram calculation. 3. Do the band structure calculation non-self-consistently, i.e. a Harris calculation where the electron density is not updated. Note that the Fermi levels for the two calculations may differ (by a significant amount !!). In general the Fermi energy from the original calculation should be used. harris_ .. _harris: inline : harris.py The resulting plot looks like: |harrisplot| .. |harrisplot| image:: harris.gif :width: 400 ------------------- Nudged Elastic Band ------------------- The Nudged Elastic Band method is a technique for finding transition paths (and corresponding energy barriers) between given initial and final states. The method constructs a chain of "images" of the system. See the ASE `NudgedElasticBand filter`_ for more details. To illustrate this functionality, we will consider O diffusion on Al(111) between an fcc and hcp hollow site. This example use the `quasi-Newton`_ algorithm for minimizing the band. See also `Mimimizing the Nudged Elastic Band`_ for a general description of this method: neb_ .. _neb: {{attachment:neb.py}} The initial and final state looks like |initial| |final| .. |initial| image:: initial.jpg :width: 200 .. |final| image:: final.jpg :width: 200 After the run a set of dacapo files out.*.nc are produced and a set of trajectory files image*.nc, one for each image. Each of these files will contain the number of timesteps done by the quasi-Newton algorithm. You can check the overall convergence from the quasi-Newton log file:: > grep QN qn.log QN: energy maxforce -8013.409959 2.705231 QN: energy maxforce -8014.081961 0.704026 QN: energy maxforce -8014.149499 0.152324 QN: energy maxforce -8014.151499 0.042264 In this case the minimization is done in 4 steps, each steps involves calculating the energy and force for each image. A restart script can look like this: restart-neb_ .. _restart-neb: attachment: restart-neb.py This script reuse the Hessian matrix written (by default) to the file ``hessian.pickle``. ``restart-neb.py`` write to the log file `qn-restart.log`:: > grep QN restart.log QN: energy maxforce -8014.151534 0.042264 Running the restart-neb.py script should not involve doing any new calculations, because the band has allready been converged in the neb.py script. At the end of ``neb.py`` and ``neb-restart.py`` the last timestep for each image are put in a new trajectory file ``neb-path.nc``. Instead of a timeseries this trajectory contains the M images for the converged path. (This last step will probably be added to ASE.Utilities) The result of: plottrajectory neb-path.nc should look like this, showing the point close to the transition state: |nebpath| .. |nebpath| image:: nebpath.gif :width: 400 -------------------------- Vibrational Calculation -------------------------- You can calculate the vibrational modes of a ListOfAtoms in the harmonic approximation using the ASE.Utilities.VibrationalCalculation module. The ListOfAtoms should either be in a ground state (i.e. fully relaxed) or at a saddle point (i.e. a transition state, as calculated by NEB). You can select a subset of the ListOfAtoms for which to find the vibrational modes, specify the displacements for each atom, and whether forward, backward or centered finite differences are used to calculate the Hessian matrix. The module calculates vibrational frequencies, zero-point energy and the vibrational eigenmodes of the system you specify. The module should be calculator independent, but it does rely on the calculator having a ReadAtoms method, and it uses netcdf files to get data from. **Note**: problems with running vibrational calculations in dacaop were reported ``_, and running dacapo in StayAlive (``_) mode turned out to be a workaround. You must run the `CO_relaxed_in_a_box.py `__ script above before you can calculate the vibrational frequencies of CO. CO_vibrations_ .. _CO_vibrations: {{attachment:CO_vibrations.py}} The output of this script looks like this:: 2124 1/cm Imaginary frequency Imaginary frequency 4 1/cm 398 1/cm 401 1/cm Zero-point energy = 0.18 eV The main experimental vibrational frequency_ for CO is 2180 1/cm. Note the calculation produces 6 frequencies, while for CO we expect only one (3N-5 because its linear). The module calculates 6 modes, but 3 of them are translational; these modes should have frequencies of 0, and correspond to the 3 smallest (including the imaginary ones, which likely had small, but negative eigenvalues. There are also 3 rotational modes for molecules (only 2 for linear molecules). The other 2 modes correspond to the rotational modes, and the magnitude of them suggests that the box is not big enough (in other words, the CO molecules are close enough to interact). Try this set of calculations in a bigger box, you will see it makes a difference. On slabs it can be a little different. If you only look at the modes of an adsorbate, then depending on the potential energy surface, the three translational modes can be converted to three frustrated translations, and likewise with the three rotational modes. So an adsorbate could have up to 3N vibrational modes. You should read the source file in detail to see all the options and things that can be done with it. You can find it in ASE/Utilities/VibrationalCalculation.py .. _frequency: "http://webbook.nist.gov/cgi/cbook.cgi?Formula=CO&NoIon=on&Units=SI&cIR=on&cES=on#IR-Spec" --------------------------------------------- Dipole correction and electrostatic potential --------------------------------------------- Using the method calc.DipoleCorrection() or the keyword ``dipole`` you can activates the interslab dipole electrostatic decoupling. Presently, this feature is only active along the third unit cell vector so you need to choose unit cell accordingly. Corrections for higher electrostatic multipoles are not implemented presently. This correction is calculated for a O atom on a Al(111) surface in the example: electrostatic_ .. _electrostatic: {{attachment:electrostatic.py}} You can get information on the dipole correction by using: grep DIP Al111-O.out This will give the output:: DIP: u_damp u Field_1 Field_2 Work_f1 Work_f2 V_jump E_dip DIP: [eA] [eA] [V/A] [V/A] [eV] [eV] [V] [eV] DIP: 0.000 0.142 - - - - 0.000 0.000 .. DIP: 0.338 0.339 0.000 0.000 4.178 6.331 2.155 3.256 DIP: 0.338 0.338 0.000 0.000 4.178 6.332 2.155 3.256 `electrostatic.py` plots the electrostatic potential: |dipole| .. |dipole| image:: dipole.gif :width: 500 For more details see the `workfunction document`_ .. _workfunction document: attachment: workfunction.pdf -------------------- STM images - Al(100) -------------------- This example use the ASE `STMTool`_ class to make a STM images of a Al(100) surface: STM.py_ .. _STM.py: attachment: STM.py For more details see the `STM tutorial`_. The resulting image is shown below: |stm| .. |stm| image:: stm.jpg :width: 300 .. _STM Tutorial: ASE:ASE_Tutorial_3 .. _STMTool: ASE:Analysis#stmtool -------------------------- Localized Wannier orbitals -------------------------- Below are examples of using the ASE `Wannier`_ class. This examples shows how to construct Wannier functions for an ethylene molecule (C2H4): Wannier-ethylene_ .. _Wannier-ethylene: {{attachment:Wannier-ethylene.py}} This example shows how to construct Wannier functions for a linear chain of 4 Pt with periodic boundary conditions. The gamma-point approximation is applied: Wannier-Pt4_ .. _Wannier-Pt4: [[attachment:Wannier-Pt4.py]] This example shows how to construct Wannier functions for an infinite,linear chain of Pt atoms. The Wannier functions are constructed using k-points: Wannier-Ptwire_ .. _Wannier-Ptwire: [[attachment:Wannier-Ptwire.py]] This example shows how to construct Wannier functions for iron-bcc. Wannier functions are constructed for spin up and spin down states separately: Wannier-Fe-bcc_ .. _Wannier-Fe-bcc: [[attachment:Wannier-Fe-bcc.py]] For more details see the `Wannier tutorial`_. .. _Wannier Tutorial: ASE:ASE_Tutorial_4 .. _Wannier: ASE:Utilities#wannier ---------------------- Transport calculations ---------------------- This example illustrates the use of the transport module by application to a simple 1d tight-binding model: transport_1dmodel_ .. _transport_1dmodel: {{attachment:transport_1dmodel.py}} For more details see the `Transport tutorial`_. This example shows how to construct a Wannier function Hamiltonian matrix to be used for a subsequent transport calculation: setupham_ .. _setupham: {{attachment:setupham.py}} For more details see the `Transport tutorial`_. ------------------------------------------------------ Bayesian Error estimation in Density Functional Theory ------------------------------------------------------ Two examples of error estimation using ensembles are given below. For more details see the `GPAW Ensemble Tutorial`_. Error estimation of the binding energy of the H2 molecule: bee_ .. _bee: {{attachment:bee.py}} Error estimation of the bond distance and frequency of the H2 molecule: bee2_ .. _bee2: {{attachment:bee2.py}} .. _Transport tutorial: ASE:ASE_Tutorial_5 .. _MDmin: ASE:Dynamics#mdmin .. _NudgedElasticBand filter: https://wiki.fysik.dtu.dk/ase/Filters#nudged-elastic-band-method .. _Subset: https://wiki.fysik.dtu.dk/ase/Filters#subset .. _plottrajectory: https://wiki.fysik.dtu.dk/ase/Trajectories#the-plottrajectory-tool .. _NetCDF trajectory: https://wiki.fysik.dtu.dk/ase/Trajectories#trajectories .. _FixBondLenght: https://wiki.fysik.dtu.dk/ase/Filters#fixbondlength .. _VibrationalCalculation module: https://wiki.fysik.dtu.dk/ase/Analysis#vibration-analysis .. _quasi-Newton: https://wiki.fysik.dtu.dk/ase/Dynamics#quasi-newton .. _Mimimizing the Nudged Elastic Band: https://wiki.fysik.dtu.dk/ase/Dynamics#mimimizing-the-nudged-elastic-band .. _EquationOfState Module: https://wiki.fysik.dtu.dk/ase/Utilities#the-equationofstate-module .. _GeometricTransforms Module: https://wiki.fysik.dtu.dk/ase/Utilities#the-geometrictransforms-module .. _GPAW Ensemble Tutorial: GPAW:Ensembles