1#!/usr/bin/env python
2
3'''
4User defined XC functional
5
6See also
7* The parser parse_xc function implemented in pyscf.dft.libxc
8* Example 24-define_xc_functional.py to input a functional which is not
9  provided in Libxc or XcFun library.
10'''
11
12from pyscf import gto
13from pyscf import dft
14
15mol = gto.M(
16    atom = '''
17    O  0.   0.       0.
18    H  0.   -0.757   0.587
19    H  0.   0.757    0.587 ''',
20    basis = 'ccpvdz')
21
22#
23# DFT can parse the custom XC functional, following the rules:
24# * The given functional description must be a one-line string.
25# * The functional description is case-insensitive.
26# * The functional description string has two parts, separated by ",".  The
27#   first part describes the exchange functional, the second is the correlation
28#   functional.
29#   - If "," was not presented in string, the entire string is treated as a
30#     compound XC functional (including both X and C functionals, such as b3lyp).
31#   - To input only X functional (without C functional), leave the second part
32#     blank. E.g. description='slater,' for pure LDA functional.
33#   - To neglect the X functional (only input C functional), leave the first
34#     part blank. E.g. description=',vwn' means pure VWN functional
35#   - If compound XC functional is specified, no matter whehter it is in the X
36#     part (the string in front of comma) or the C part (the string behind
37#     comma), both X and C functionals of the compound XC functional will be
38#     used.
39# * The functional name can be placed in arbitrary order.  Two names needs to
40#   be separated by operators + or -.  Blank spaces are ignored.
41#   NOTE the parser only reads operators + - *.  / is not supported.
42# * A functional name can have at most one factor.  If the factor is not
43#   given, it is set to 1.  Compound functional can be scaled as a unit. For
44#   example '0.5*b3lyp' is equivalent to
45#   'HF*0.1 + .04*LDA + .36*B88, .405*LYP + .095*VWN'
46# * String "HF" stands for exact exchange (HF K matrix).  It is allowed to
47#   put "HF" in C (correlation) functional part.
48# * String "RSH" means range-separated operator. Its format is
49#   RSH(omega, alpha, beta).  Another way to input RSH is to use keywords
50#   SR_HF and LR_HF: "SR_HF(0.1) * alpha_plus_beta" and "LR_HF(0.1) *
51#   alpha" where the number in parenthesis is the value of omega.
52# * Be careful with the libxc convention on GGA functional, in which the LDA
53#   contribution has been included.
54
55mf = dft.RKS(mol)
56# B3LYP can be constructed
57mf.xc = 'HF*0.2 + .08*LDA + .72*B88, .81*LYP + .19*VWN'
58e1 = mf.kernel()
59print('E = %.15g  ref = -76.3832244350081' % e1)
60
61#
62# No correlation functional
63#
64mf.xc = '.2*HF + .08*LDA + .72*B88'
65e1 = mf.kernel()
66print('E = %.15g  ref = -75.9807850596666' % e1)
67
68#
69# If not given, the factor for each functional equals 1 by default.
70#
71mf = dft.RKS(mol)
72mf.xc = 'b88,lyp'
73e1 = mf.kernel()
74
75mf = dft.RKS(mol)
76mf.xc = 'b88*1,lyp*1'
77eref = mf.kernel()
78print('%.15g == %.15g' % (e1, eref))
79
80#
81# Compound functionals can be used as part of the definition
82#
83mf = dft.RKS(mol)
84mf.xc = '0.5*b3lyp, 0.5*lyp'
85e1 = mf.kernel()
86print('E = %.15g  ref = -71.9508340443282' % e1)
87
88# Compound XC functional can be presented in the C part (the string behind
89# comma). Both X and C functionals of the compound XC functional will be used.
90# Compound XC functional can be scaled as a unit.
91#
92mf = dft.RKS(mol)
93mf.xc = '0.5*b3lyp, 0.5*b3p86'
94e1 = mf.kernel()
95
96mf = dft.RKS(mol)
97mf.xc = '0.5*b3lyp + 0.5*b3p86'
98e2 = mf.kernel()
99print('E1 = %.15g  E2 = %.15g  ref = -76.3923625924023' % (e1, e2))
100
101#
102# More examples of customized functionals. NOTE These customized functionals
103# are presented for the purpose of demonstrating the feature of the XC input.
104# They are not reported in any literature. DO NOT use them in the actual
105# calculations.
106#
107# Half HF exchange plus half B3LYP plus half VWN functional
108mf.xc = '.5*HF+.5*B3LYP,VWN*.5'
109
110# "-" to subtract one functional from another
111mf.xc = 'B88 - SLATER*.5'
112
113# The functional below gives omega = 0.33, alpha = 0.6 * 0.65 = 0.39
114# beta = -0.46 * 0.6 + 0.4 * 0.2(from HF of B3P86) = -0.196
115mf.xc = '0.6*CAM_B3LYP+0.4*B3P86'
116
117# The input XC description does not depend on the order of functionals. The
118# functional below is the same to the functional above
119mf.xc = '0.4*B3P86+0.6*CAM_B3LYP'
120
121# Use SR_HF/LR_HF keywords to input range-separated functionals.
122# When exact HF exchange is presented, it is split into SR and LR parts
123# alpha = 0.8(from HF) + 0.22
124# beta = 0.5 + 0.8(from HF)
125mf.xc = '0.5*SR-HF(0.3) + .8*HF + .22*LR_HF'
126
127# RSH is another keyword to input range-separated functionals
128mf.xc = '0.5*RSH(0.3,2.04,0.56) + 0.5*BP86'
129
130# A shorthand to input 'PBE,PBE', which is a compound functional. Note the
131# shorthand input is different to the two examples 'PBE,' and ',PBE' below.
132mf.xc = 'PBE'
133
134# Exchange part only
135mf.xc = 'PBE,'
136
137# Correlation part only
138mf.xc = ',PBE'
139
140# When defining custom functionals, compound functional will affect both the
141# exchange and correlation parts
142mf.xc = 'PBE + SLATER*.5'
143
144# The above input is equivalent to
145mf.xc = 'PBE + SLATER*.5, PBE'
146
147