Mod:Hunt Research Group/OTRSH Funct

Links to files

Zip Archive: Media:OTRSH_hpc_ver_V2_Mar2018.zip

What's an OTRSH Functional?

It is an optimally tuned range separated hybrid functional (OTRSH). This is a class of functionals that are useful for ccalculating electronic properties such as UV-Vis spectra (especially if charge-transfer excitations are important) or photoelectron spectra. Before talking about OTRSH functionals, it is useful to first talk about how standard range separated hybrid functionals work, and why they were created.

The Problem

Standard DFT has problems describing long range interactions, due to the form of the exchange energy in standard functionals. Actual (exact) Exchange energy scales with distance as $\frac{1}{r_{12}}$ where $r_{12}$ is the distance between two electrons. However, standard pure DFT exchange functionals depend directly on the electron density, which decays exponentially with distance from an atomic center (i.e. DFT exchange decays faster than actual exchange). This means that standard DFT functionals are poor at describing long-range interactions (charge transfer excitations in DFT is an example). Range seperated hybrid functionals try to fix this by using Hartree-Fock exchange energy (which is exact) at large distances (long-range) while keeping DFT exchange energy at short range.

One Solution - Normal Range Separated Hybrids

A range separated hydrbrid functional divides the exchange energy into two terms - short range exchange ( $E_{X}^{S R}$ ) and long range exchange ( $E_{X}^{L R}$ ). The exchange energy at any one distance is a weighted combination of short-range and long range exchange, and can essentially be described with three variables; $α$ , $β$ and $μ$ (See picture below).

The equations for total, short range (DFT) and long-range (HF) exchange are:

$E_{X}^{T o t a l} = E_{X}^{S R} + E_{X}^{L R}$

$E_{X}^{S R} = [1 - (α + β . erf (μ r_{12})] E_{X}^{D F T} (r_{12})$

$E_{X}^{L R} = [α + β . erf (μ r_{12})] E_{X}^{H F} (r_{12})$

where erf is the error function. $α$ represents the amount of HF exchange at short range. For example, $α = 0.2$ would mean 20% HF exchange was used at $r_{12} = 0$ . $α + β$ represents the amount of HF exchange at long range (i.e as $r_{12}$ approaches infinity). A value of $α + β = 1$ means that 100% HF exchange is employed at long range, which should give the correct long-range behavior for the exchange energy. $μ$ is the range seperation parameter, its value determines how quickly we change from short-range to long-range exchange expresssions. The higher $μ$ is, the faster (shorter distances) we start using the long-range exchange energy expression.

A problem with range separated hybrids is in picking the correct value for $μ$ (i.e when should we stop using DFT exchange?). This is often determined empirically by varying $μ$ for a test data set, and seeing which $μ$ value leads to a functional that best reproduces an experimental property (e.g. atomic ionisation energies). This has the same problem as any kind of empirical fitting - its not clear how well the functional will work outside the test set of molecules/properties. Also this empirical tuning of parameters can lead to unphysical values. For example, CAM-B3LYP uses $α + β = 0.65$ , meaning it definitely has the wrong long-range behavior.

OTRSH Functionals

OTRSH functionals provide a non-empirical method for choosing $μ$ ; $μ$ is "tuned" on a per-calculation basis. The $μ$ value is chosen so that the calculated ionisation energy matches the negative of the HOMO energy, i.e.:

$I P_{C a l c} = - E_{H O M O}$

$I P_{C a l c} = E_{N - 1} - E_{N}$

Where $E_{N - 1}$ and $E_{N}$ represent the energy of the N-1 and N electron system (both in the N electron system geometry). Note that all values in the equations above can be obtained from two single-point energy calculations. Also note that $I P_{C a l c} = - E_{H O M O}$ is true for the exact (unknown) KS-DFT functional, hence it is a sensible thing to try to reproduce in a functional.

To summarise an OTRSH functional is one whereby $μ$ is chosen to minimise $I P_{C a l c} + E_{H O M O}$ , hence the $μ$ value will be different for every calculation you carry out. Sometimes the IP and HOMO energy of the N-1 electron system are also used to tune $μ$ , but unfortunately thats not possible the code on this page.

Running The Code

Quick Start Guide (Using the test file)

Once you download the .zip archive, copy the files to somewhere on your hpc path [Note: make sure you dont miss the .rng_OTRSH_V2 file]. Personally i use "/home/rf614/bin" as the location. To use the gf_OTRSH_V2 function straight away you'll need to put them in $HOME/bin. "echo $HOME" while on the hpc to find where your home directory is if your not sure.

Move the test file ('H2O_SPE_6311plgdp_GP.com') to whereever you like in your work directory, and navigate to this folder. Then type either of the following commands (assuming you put the scripts in ~/bin):

qsub -v in=H2O_SPE_6311plgdp_GP -N test ~/bin/.rng_OTRSH_V2
gf_OTRSH_V2 H2O_SPE_6311plgdp_GP.com

If you saved the files somewhere else then simply modify "~/bin/.rng_OTRSH" to whatevers appropriate for the top command. For details on how the gf command works see this wiki link.

Quick Start Guide - The output files

Assuming you ran the job above you will now have a loada output files in your current directory, i.e.:

H2O_SPE_6311plgdp_GP_gsstate_mu_opt.*

H2O_SPE_6311plgdp_GP_ionstate_mu_opt.*

H2O_SPE_6311plgdp_GP_OTRSH.out

The H2O_SPE_6311plgdp_GP_OTRSH.out file is essentially a progress file for how the tuning process is going. Once the calculation is over its main purpose is really to keep track of the options you used to create the OTRSH functional. Its a comma delimited file, and the headings should be pretty self-explanatory (at least they will be if you read the above section on "What is an OTRSH")

The H2O_SPE_6311plgdp_GP_gsstate_mu_opt* files are for single energy calculations using the OTRSH functional tuned for this particular system. If your goal was to calculate ground state properties with an OTRSH then these files might be all you need. If you want to do something other than a single-point energy calculation, then just add the relevant keywords (e.g TD) to the H2O_SPE_6311plgdp_GP_gsstate_mu_opt.com file - this file contains all the keywords needed to define the OTRSH functional.

The H2O_SPE_6311plgdp_GP_ionstate_mu_opt* files are the single point energy calculations (using the OTRSH) for the ionised system. You probably have no need for these.

You will also have a new folder named temp_files_H2O_SPE_6311plgdp_GP. At every tested value for $μ$ , the .log files for the ground and ionised state are written to this folder. The purpose of this folder is to check nothings going wrong during the calculation. If something stops the code extracting scf energies for each $μ$ value, then it might just loop infinitely (or at least for a loooooong time). The code SHOULD catch these sorts of errors and kill the job, but if it runs a long time without any change in the *.out file then make sure the *.log files in this folder are ok. The secondary purpose of this folder is therefore to tell you WHY your job got killed, if it got killed (an error message should be given in the .out file as well but....)

Quick Note on Creating Input Files

If you want to run the code, with default options, on your own systems then theres not much to say. Simply prepare a normal SINGLE POINT ENERGY *.com file for your system, but remove the functional from the input line. Also dont use the "#p" command for excess printing - this will break the OTRSH code (sorry, the code essentially uses the # for a search + replace to insert a functional at each step of the tuning process). For example, this input line is fine(and used in the example file):

# 6-311+g(d,p) geom=connectivity int=ultrafine scf=conver=9

Make sure you use a single point energy calculation input file, using a keyword like "freq" wont throw any errors, but you'll be waiting a really long time for your job to finish.

Understanding (And Modifying) The Default Options

The easily changeable variables in this code are stored in the .rng_OTRSH_V2 file. They are (in order of "likelyhood of needing altering"):

es_mult=2
alpha_val=0
beta_val=0
mu_lower=0.1
mu_upper=0.2
max_step=0.03
min_step=0.001
xc_funct="LC-wPBE"

Its likely you'll only ever need/want to change the es_mult here. This is the multiplicity of the ionised state, and is obviously equal to 2 if you have a closed shell system. If you use this code on an open shell system you need to set es_mult to the appropriate value for the ionised state. Personally i create seperate .rng files for each value of multiplicity i use, to save me editing the .rng_OTRSH file too often.

alpha_val corresponds to the fraction HF exchange at short range. Setting this to zero means use the default value for the functional.Its plausible that raising this to say alpha_val=0.2 might improve your results in some cases, so it might be worth trying the functional with two different alpha values if your trying to match an experimental spectrum or have something else you can validate with.

beta_val - as stated above "alpha_val + beta_val" corresponds to the fraction HF exchange at long range. Setting this to zero means use the default value for the functional. If alpha_val and beta_val are both zero, then only the $μ$ value should differ between the base functional and the OTRSH

mu_lower and mu_upper control the initial range of $μ$ values to try out. If you have a good idea what your final $μ$ will be, then settings mu_lower and mu_upper to nearby values should speed up the optimisation process. However, there is NO NEED for the final mu_val to be between these initial search values, it will just take longer if its far away from this range.

max_step and min_step correspond to the maximum and minimum change in $μ$ between optimisation steps. Note that i think the minimum possible min_step=0.0001, so dont set it lower than that or it could cause serious problems (it'll lead to a gradient of zero between two mu values, i have no idea how the optimiser will react to that but it wont be good). Theres probably NEVER a need to mess with these values. If the optimiser isnt efficient enough, changing these values wont help much. Writing a better optimiser (the current ones pretty terrible by the way,but it works so.....) would probably be easier than working out the best values for these numbers.

LC-wPBE corresponds to the underlying functional to use. The functional this code generates is essentially LC-wPBE with a modified $μ$ value. The underlying DFT exchange/ DFT correlation functionals are the same as those used in LC-wPBE. The underlying functional needs to be a range separated one (i think). If you change it, make sure you can reproduce the base functional with the code (e.g. run the code, change mu to that of the base functional and see if the results are the same).

Finally, the convergence criteria for the OTRSH (the IP_HOMO - E_HOMO value) is set to 0.01 eV by default, and controlled by the IP_conver_thres variable in optimiser.py. I guess its possible you may want to relax that criteria in certain circumstances (obviously it would speed up the optimisation and lower the accuracy).