Tuning the range separation

In this example we optimize the value of gamma parameter for long-range corrected XC functional (in our case: LCY-PBE) in ADF. Long-range corrected XC functionals can be used in ADF with XCfun (see ADF manual).

The optimal range separation parameter gamma yields the HOMO energy equal to the ionization potential (IP). Given a molecular system, we simultaneously minimize the difference between HOMO and IP for that system (N) and its anion (A) (system with one more electron). We define the J function as:

\[J = \sqrt{N^2+A^2}\]

and find the value of gamma (within a certain range) which minimizes J. See also this article by Kronik and coworkers.

We first define a new job type GammaJob by extending MultiJob. The goal of GammaJob is to calculate the J function for one fixed value of gamma To do that we need to perform 3 different single point calculations: 1 for the given system (let’s call it S0), 1 for the system with one more electron (S-) and 1 for the system with one less electron (S+). S+ calculation is needed to find the ionization potential of S0.

The constructor (__init__) of GammaJob accepts several new arguments and simply stores them. These new arguments define: the value of gamma, the Molecule together with its initial charge, and the values of spin for S-, S0 and S+ (as a tuple of length 3). Then the prerun() method is used to create three children jobs with different values of total charge and spin multiplicity. A dedicated Results subclass features a simple method for extracting the value of J based on results on three children jobs.

We can then treat our newly defined GammaJob as a blackbox with simple interface: input gamma -> run -> extract J. The next step is to create multiple instances of GammaJob for a range of different gammas. That task can be conveniently wrapped in a simple function gamma_scan.

class GammaResults(Results):

    def get_difference(job, jobplus):
        """Calculate the difference between HOMO and IP.
        *jobplus* should be the counterpart of *job* with one less electron."""
        homo = job.results.readrkf('Properties','HOMO', file='engine')
        IP = jobplus.results.get_energy() - job.results.get_energy()
        return IP + homo

    def get_J(self):
        N = GammaResults.get_difference(self.job.children[1], self.job.children[2])
        A = GammaResults.get_difference(self.job.children[0], self.job.children[1])
        return (N**2 + A**2)**0.5

class GammaJob(MultiJob):
    _result_type = GammaResults

    def __init__(self, molecule, gamma, charge, spins, **kwargs):
        MultiJob.__init__(self, **kwargs)
        self.molecule = molecule
        self.charge = charge
        self.spins = spins
        self.gamma = gamma

    def prerun(self):
        charges = [self.charge-1, self.charge, self.charge+1]
        for charge, spin in zip(charges, self.spins):
            name = '{}_charge_{}'.format(self.name, charge)
            newjob = AMSJob(name=name, molecule=self.molecule, settings=self.settings)
            newjob.settings.input.ams.System.charge = charge
            newjob.settings.input.adf.xc.rangesep = "gamma={:f}".format(self.gamma)
            if spin != 0:
                newjob.settings.input.adf.unrestricted = True
                newjob.settings.input.adf.SpinPolarization = spin


def gamma_scan(gammas, settings, molecule, name='scan', charge=0, spins=(1,0,1)):
    """Calculate values of J function for given range of gammas.

    gammas   - list of gamma values to calculate the J function for
    settings - Settings object for an ADF calculation
    molecule - Molecule object with the system of interest
    name     - base name of all the jobs
    charge   - base charge of the system of interest. The J function is going to be
               calculated based on two systems: with charge, and charge-1
    spins    - values of spin polarization for jobs with, respectively, charge-1, 
               charge and charge +1

    In other words, if charge=X and spins=(a,b,c) the three resulting jobs
    are going to have the following values for charge and spin:

    Charge=X-1  SpinPolarization=a
    Charge=X    SpinPolarization=b
    Charge=X+1  SpinPolarization=c

    Returns a list of pairs (gamma, J) of the same length as the parameter *gammas*
    jobs = [GammaJob(molecule=molecule, settings=settings, gamma=g,
            charge=charge, spins=spins, name=name+'_gamma_'+str(g)) for g in gammas]
    results = [j.run() for j in jobs]
    js = [r.get_J() for r in results]
    return list(zip(gammas, js))

# =============================================================
# Now we simply use the gamma_scan function to find the optimal 
# gamma value for a toy system (H2)
# =============================================================

import numpy as np
import multiprocessing

# Run as many jobs in parallel as there are cores:
config.default_jobrunner = JobRunner(parallel=True, maxjobs=multiprocessing.cpu_count())

# Settings of the ADF calculations
# ================================

s = Settings()
s.input.ams.task = 'SinglePoint'
s.input.adf.basis.type = 'DZP'
s.input.adf.basis.core = 'None'
s.input.adf.xc.gga = 'PBE'
s.input.adf.xc.xcfun = True
s.runscript.nproc = 1

# The molecule (here we just use H2)
# ==================================

mol = Molecule()
mol.add_atom(Atom(symbol='H', coords=(0,0,-0.3540)))
mol.add_atom(Atom(symbol='H', coords=(0,0, 0.3540)))

# The list of gamma values
# ========================

# Here we scan just a few values for gamma. 
# In practice, you want to scan a wider range and smaller step.
gammas = np.around(np.arange(1.2, 1.9, 0.2), decimals=3)

results = gamma_scan(gammas, s, mol)

print('== Results ==')
print('gamma \t J')
for g,j in results:
    print('{:.4f} \t {:.8f}'.format(g,j))
print('Optimal gamma value: {:.4f}'.format(min(results,key=lambda x:x[1])[0]))


To execute this PLAMS script:


[17:24:51] JOB scan_gamma_1.2 STARTED
[17:24:51] JOB scan_gamma_1.4 STARTED
[17:24:51] JOB scan_gamma_1.6 STARTED
[17:24:51] JOB scan_gamma_1.2 RUNNING
[17:24:51] JOB scan_gamma_1.4 RUNNING
[17:26:21] JOB scan_gamma_1.8 FINISHED
[17:26:21] JOB scan_gamma_1.8 SUCCESSFUL
== Results ==
gamma 	 J
1.2000 	 0.01138875
1.4000 	 0.00755822
1.6000 	 0.00858682
1.8000 	 0.01111191
Optimal gamma value: 1.4000
[17:26:21] PLAMS run finished. Goodbye
Test duration in seconds: 90