Improve XC parser: if the XC description does not contain ",", the de…

…scription is treated as XC

examples/mpi/00-pbc_df.py

@@ -39,7 +39,7 @@
 kmf = dft.KRKS(cell, kpts)
 # Turn to the atomic grids if you like
 kmf.grids = dft.gen_grid.BeckeGrids(cell)
-kmf.xc = 'm06'
+kmf.xc = 'm06,m06'
 kmf.with_df = mpidf.FFTDF(cell, kpts)
 kmf.kernel()
 

examples/pbc/10-gamma_point_scf.py

@@ -35,7 +35,7 @@
 print("HF energy (per unit cell) = %.17g" % ehf)
 
 mf = dft.RKS(cell)
-mf.xc = 'm06'
+mf.xc = 'm06,m06'
 edft = mf.kernel()
 print("DFT energy (per unit cell) = %.17g" % edft)
 

examples/pbc/20-k_points_scf.py

@@ -34,7 +34,7 @@
 kmf = dft.KRKS(cell, kpts)
 # Turn to the atomic grids if you like
 kmf.grids = dft.gen_grid.BeckeGrids(cell)
-kmf.xc = 'm06'
+kmf.xc = 'm06,m06'
 kmf.kernel()
 
 

pyscf/dft/libxc.py

@@ -292,7 +292,7 @@
 'XC_HYB_GGA_XC_HJS_B88'        :  431, # HJS hybrid screened exchange B88 version
 'XC_HYB_GGA_XC_HJS_B97X'       :  432, # HJS hybrid screened exchange B97x version
 'XC_HYB_GGA_XC_CAM_B3LYP'      :  433, # CAM version of B3LYP
-'XC_HYB_GGA_XC_TUNED_CAM_B3LYP' :  434, # CAM version of B3LYP tuned for excitations
+'XC_HYB_GGA_XC_TUNED_CAM_B3LYP':  434, # CAM version of B3LYP tuned for excitations
 'XC_HYB_GGA_XC_BHANDH'         :  435, # Becke half-and-half
 'XC_HYB_GGA_XC_BHANDHLYP'      :  436, # Becke half-and-half with B88 exchange
 'XC_HYB_GGA_XC_MB3LYP_RC04'    :  437, # B3LYP with RC04 LDA
@@ -442,6 +442,18 @@
 'CAMYB3LYP'     : 'XC_HYB_GGA_XC_CAMY_B3LYP',
 }
 
+# Some XC functionals have conventional name, like M06-L means M06-L for X
+# functional and M06-L for C functional, PBE mean PBE-X plus PBE-C. If the
+# conventional name was placed in the XC_CODES, it may lead to recursive
+# reference when parsing the xc description.  These names (as exceptions of
+# XC_CODES) are listed in XC_ALIAS below and they should be treated as a
+# shortcut for XC functional.
+XC_ALIAS = {
+    # Conventional name : name in XC_CODES
+    'M06-L'             : 'M06_L,M06_L',
+    'PBE'               : 'PBE,PBE'
+}
+
 XC_KEYS = set(XC_CODES.keys())
 
 VV10_XC = set(('B97M_V', 'WB97M_V', 'WB97X_V', 'VV10', 'LC_VV10'))
@@ -578,10 +590,12 @@ def parse_xc(description):
       first part describes the exchange functional, the second is the correlation
       functional.
 
-      - If "," was not appeared in string, the entire string is considered as
-        X functional.
+      - If "," was not in string, the entire string is considered as a XC
+        functional (if applicable).
+      - To apply only X functional (without C functional), leave blank in the
+        second part, e.g. description='lda,' for pure LDA functional
       - To neglect X functional (just apply C functional), leave blank in the
-        first part, eg description=',vwn' for pure VWN functional
+        first part, e.g. description=',vwn' for pure VWN functional
       - If compound XC functional (including both X and C functionals, such as
         b3lyp) is specified, no matter whehter it is in the X part (the string
         in front of comma) or the C part (the string behind comma), both X and C
@@ -768,25 +782,24 @@ def remove_dup(fn_facs):
                 n += 1
         return list(zip(fn_ids, facs))
 
+    description = description.replace(' ','').upper()
+    if ',' not in description and description in XC_ALIAS:
+        description = XC_ALIAS[description]
+
     if '-' in description:  # To handle e.g. M06-L
         for key in _NAME_WITH_DASH:
             if key in description:
                 description = description.replace(key, _NAME_WITH_DASH[key])
 
     if ',' in description:
-        x_code, c_code = description.replace(' ','').upper().split(',')
+        x_code, c_code = description.split(',')
         for token in x_code.replace('-', '+-').split('+'):
             parse_token(token, possible_x_k_for)
         for token in c_code.replace('-', '+-').split('+'):
             parse_token(token, possible_c_for)
     else:
-        x_code = description.replace(' ','').upper()
-        try:
-            for token in x_code.replace('-', '+-').split('+'):
-                parse_token(token, possible_xc_for)
-        except KeyError:
-            for token in x_code.replace('-', '+-').split('+'):
-                parse_token(token, possible_x_k_for)
+        for token in description.replace('-', '+-').split('+'):
+            parse_token(token, possible_xc_for)
     return hyb, remove_dup(fn_facs)
 
 _NAME_WITH_DASH = {'SR-HF'  : 'SR_HF',

pyscf/dft/test/test_dks.py

@@ -28,7 +28,7 @@ def test_dks_lda(self):
         mol.output = '/dev/null'
         mol.build()
         mf = dks.DKS(mol)
-        mf.xc = 'lda'
+        mf.xc = 'lda,'
         eks4 = mf.kernel()
         self.assertAlmostEqual(eks4, -126.041808355268, 9)
 

pyscf/dft/test/test_libxc.py

@@ -66,7 +66,7 @@ def test_parse_xc(self):
         hyb, fn_facs = dft.libxc.parse_xc('APBE,')
         self.assertEqual(fn_facs[0][0], 184)
 
-        hyb, fn_facs = dft.libxc.parse_xc('TF')
+        hyb, fn_facs = dft.libxc.parse_xc('TF,')
         ref = [(1, 1), (7, 1)]
         self.assertEqual(dft.libxc.parse_xc_name('LDA,VWN'), (1,7))
         self.assertEqual(dft.libxc.parse_xc(('LDA','VWN'))[1], ref)
@@ -76,12 +76,12 @@ def test_parse_xc(self):
 
         self.assertTrue (dft.libxc.is_meta_gga('m05'))
         self.assertFalse(dft.libxc.is_meta_gga('pbe0'))
-        self.assertFalse(dft.libxc.is_meta_gga('tf'))
+        self.assertFalse(dft.libxc.is_meta_gga('tf,'))
         self.assertFalse(dft.libxc.is_meta_gga('vv10'))
         self.assertTrue (dft.libxc.is_gga('PBE0'))
         self.assertFalse(dft.libxc.is_gga('m05'))
-        self.assertFalse(dft.libxc.is_gga('tf'))
-        self.assertTrue (dft.libxc.is_lda('tf'))
+        self.assertFalse(dft.libxc.is_gga('tf,'))
+        self.assertTrue (dft.libxc.is_lda('tf,'))
         self.assertFalse(dft.libxc.is_lda('vv10'))
         self.assertTrue (dft.libxc.is_hybrid_xc('m05'))
         self.assertTrue (dft.libxc.is_hybrid_xc('pbe0,'))

pyscf/dft/test/test_numint.py

@@ -193,7 +193,7 @@ def test_rks_vxc(self):
         numpy.random.seed(10)
         nao = mol.nao_nr()
         dms = numpy.random.random((2,nao,nao))
-        v = mf._numint.nr_vxc(mol, mf.grids, 'B88', dms, spin=0, hermi=0)[2]
+        v = mf._numint.nr_vxc(mol, mf.grids, 'B88,', dms, spin=0, hermi=0)[2]
         self.assertAlmostEqual(finger(v), -0.70124686853021512, 8)
 
         v = mf._numint.nr_vxc(mol, mf.grids, 'HF', dms, spin=0, hermi=0)[2]
@@ -205,7 +205,7 @@ def test_uks_vxc(self):
         numpy.random.seed(10)
         nao = mol.nao_nr()
         dms = numpy.random.random((2,nao,nao))
-        v = mf._numint.nr_vxc(mol, mf.grids, 'B88', dms, spin=1)[2]
+        v = mf._numint.nr_vxc(mol, mf.grids, 'B88,', dms, spin=1)[2]
         self.assertAlmostEqual(finger(v), -0.73803886056633594, 8)
 
         v = mf._numint.nr_vxc(mol, mf.grids, 'HF', dms, spin=1)[2]
@@ -223,20 +223,20 @@ def test_rks_fxc(self):
         dm0 = numpy.einsum('pi,i,qi->pq', mo_coeff, mo_occ, mo_coeff)
         dms = numpy.random.random((2,nao,nao))
         ni = dft.numint._NumInt()
-        v = ni.nr_fxc(mol1, mf.grids, 'B88', dm0, dms, spin=0, hermi=0)
+        v = ni.nr_fxc(mol1, mf.grids, 'B88,', dm0, dms, spin=0, hermi=0)
         self.assertAlmostEqual(finger(v), -7.5671368618070343, 8)
 
         # test cache_kernel
-        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'B88', mo_coeff, mo_occ, spin=0)
-        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'B88', dm0, dms, spin=0, hermi=0,
+        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'B88,', mo_coeff, mo_occ, spin=0)
+        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'B88,', dm0, dms, spin=0, hermi=0,
                                rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
-        v = ni.nr_fxc(mol1, mf.grids, 'LDA', dm0, dms[0], spin=0, hermi=0)
+        v = ni.nr_fxc(mol1, mf.grids, 'LDA,', dm0, dms[0], spin=0, hermi=0)
         self.assertAlmostEqual(finger(v), -3.0019207112626876, 8)
         # test cache_kernel
-        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'LDA', mo_coeff, mo_occ, spin=0)
-        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'LDA', dm0, dms[0], spin=0, hermi=0,
+        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'LDA,', mo_coeff, mo_occ, spin=0)
+        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'LDA,', dm0, dms[0], spin=0, hermi=0,
                                rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
@@ -256,31 +256,31 @@ def test_rks_fxc_st(self):
         dms = numpy.random.random((2,nao,nao))
         ni = dft.numint._NumInt()
 
-        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'B88', [mo_coeff,mo_coeff],
+        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'B88,', [mo_coeff,mo_coeff],
                                  [mo_occ*.5]*2, spin=1)
-        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88', dm0, dms, singlet=True)
+        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88,', dm0, dms, singlet=True)
         self.assertAlmostEqual(finger(v), -7.5671368618070343*2, 8)
-        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88', dm0, dms, singlet=True,
+        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88,', dm0, dms, singlet=True,
                                       rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
-        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88', dm0, dms, singlet=False)
+        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88,', dm0, dms, singlet=False)
         self.assertAlmostEqual(finger(v), -7.5671368618070343*2, 8)
-        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88', dm0, dms, singlet=False,
+        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'B88,', dm0, dms, singlet=False,
                                       rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
-        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'LDA', [mo_coeff,mo_coeff],
+        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'LDA,', [mo_coeff,mo_coeff],
                                  [mo_occ*.5]*2, spin=1)
-        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA', dm0, dms[0], singlet=True)
+        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA,', dm0, dms[0], singlet=True)
         self.assertAlmostEqual(finger(v), -3.0019207112626876*2, 8)
-        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA', dm0, dms[0], singlet=True,
+        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA,', dm0, dms[0], singlet=True,
                                       rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
-        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA', dm0, dms[0], singlet=False)
+        v = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA,', dm0, dms[0], singlet=False)
         self.assertAlmostEqual(finger(v), -3.0019207112626876*2, 8)
-        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA', dm0, dms[0], singlet=False,
+        v1 = dft.numint.nr_rks_fxc_st(ni, mol1, mf.grids, 'LDA,', dm0, dms[0], singlet=False,
                                       rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
@@ -294,20 +294,20 @@ def test_uks_fxc(self):
         dm0 = numpy.einsum('xpi,xi,xqi->xpq', mo_coeff, mo_occ, mo_coeff)
         dms = numpy.random.random((2,nao,nao))
         ni = dft.numint._NumInt()
-        v = ni.nr_fxc(mol1, mf.grids, 'B88', dm0, dms, spin=1)
+        v = ni.nr_fxc(mol1, mf.grids, 'B88,', dm0, dms, spin=1)
         self.assertAlmostEqual(finger(v), -10.316443204083185, 8)
 
         # test cache_kernel
-        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'B88', mo_coeff, mo_occ, spin=1)
-        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'B88', dm0, dms, hermi=0, spin=1,
+        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'B88,', mo_coeff, mo_occ, spin=1)
+        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'B88,', dm0, dms, hermi=0, spin=1,
                                rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
-        v = ni.nr_fxc(mol1, mf.grids, 'LDA', dm0, dms[0], spin=1)
+        v = ni.nr_fxc(mol1, mf.grids, 'LDA,', dm0, dms[0], spin=1)
         self.assertAlmostEqual(finger(v), -5.6474405864697967, 8)
         # test cache_kernel
-        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'LDA', mo_coeff, mo_occ, spin=1)
-        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'LDA', dm0, dms[0], hermi=0, spin=1,
+        rvf = ni.cache_xc_kernel(mol1, mf.grids, 'LDA,', mo_coeff, mo_occ, spin=1)
+        v1 = dft.numint.nr_fxc(mol1, mf.grids, 'LDA,', dm0, dms[0], hermi=0, spin=1,
                                rho0=rvf[0], vxc=rvf[1], fxc=rvf[2])
         self.assertAlmostEqual(abs(v-v1).max(), 0, 8)
 
@@ -338,9 +338,9 @@ def test_uks_gga_wv1(self):
         rho1 = [numpy.random.random((4,5))]*2
         weight = numpy.ones(5)
 
-        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88', rho0, 1, 0, 3)
+        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88,', rho0, 1, 0, 3)
         wva, wvb = dft.numint._uks_gga_wv1(rho0, rho1, vxc, fxc, weight)
-        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88', rho0[0]+rho0[1], 0, 0, 3)
+        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88,', rho0[0]+rho0[1], 0, 0, 3)
         wv = dft.numint._rks_gga_wv1(rho0[0]+rho0[1], rho1[0]+rho1[1], vxc, fxc, weight)
         self.assertAlmostEqual(abs(wv - wva).max(), 0, 10)
         self.assertAlmostEqual(abs(wv - wvb).max(), 0, 10)
@@ -351,9 +351,9 @@ def test_uks_gga_wv2(self):
         rho1 = [numpy.random.random((4,5))]*2
         weight = numpy.ones(5)
 
-        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88', rho0, 1, 0, 3)
+        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88,', rho0, 1, 0, 3)
         wva, wvb = dft.numint._uks_gga_wv2(rho0, rho1, fxc, kxc, weight)
-        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88', rho0[0]+rho0[1], 0, 0, 3)
+        exc, vxc, fxc, kxc = dft.libxc.eval_xc('b88,', rho0[0]+rho0[1], 0, 0, 3)
         wv = dft.numint._rks_gga_wv2(rho0[0]+rho0[1], rho1[0]+rho1[1], fxc, kxc, weight)
         self.assertAlmostEqual(abs(wv - wva).max(), 0, 10)
         self.assertAlmostEqual(abs(wv - wvb).max(), 0, 10)

pyscf/dft/xcfun.py

@@ -149,6 +149,18 @@
 'TF'            : 'TFK',
 }
 
+# Some XC functionals have conventional name, like M06-L means M06-L for X
+# functional and M06-L for C functional, PBE mean PBE-X plus PBE-C. If the
+# conventional name was placed in the XC_CODES, it may lead to recursive
+# reference when parsing the xc description.  These names (as exceptions of
+# XC_CODES) are listed in XC_ALIAS below and they should be treated as a
+# shortcut for XC functional.
+XC_ALIAS = {
+    # Conventional name : name in XC_CODES
+    'M06-L'             : 'M06L,M06L',
+    'PBE'               : 'PBE,PBE'
+}
+
 LDA_IDS = set([0, 2, 3, 13, 14, 15, 24, 28, 45])
 GGA_IDS = set([1, 4, 5, 6, 7, 8, 9, 16, 17, 18, 19, 20, 21, 22, 23, 25, 26,
                27, 44, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 57, 58, 59, 61,
@@ -257,10 +269,12 @@ def parse_xc(description):
       first part describes the exchange functional, the second is the correlation
       functional.
 
-      - If "," was not appeared in string, the entire string is considered as
-        X functional.
+      - If "," was not in string, the entire string is considered as a XC
+        functional (if applicable).
+      - To apply only X functional (without C functional), leave blank in the
+        second part, e.g. description='lda,' for pure LDA functional
       - To neglect X functional (just apply C functional), leave blank in the
-        first part, eg description=',vwn' for pure VWN functional
+        first part, e.g. description=',vwn' for pure VWN functional
       - If compound XC functional (including both X and C functionals, such as
         b3lyp) is specified, no matter whehter it is in the X part (the string
         in front of comma) or the C part (the string behind comma), both X and C
@@ -288,10 +302,14 @@ def parse_xc(description):
     elif not isinstance(description, str): #isinstance(description, (tuple,list)):
         return parse_xc('%s,%s' % tuple(description))
 
+    description = description.replace(' ','').upper()
+    if ',' not in description and description in XC_ALIAS:
+        description = XC_ALIAS[description]
+
     if ',' in description:
-        x_code, c_code = description.replace(' ','').upper().split(',')
+        x_code, c_code = description.split(',')
     else:
-        x_code, c_code = description.replace(' ','').upper(), ''
+        x_code, c_code = description, ''
 
     def assign_omega(omega):
         if hyb[2] == 0:

pyscf/geomopt/berny_solver.py

@@ -141,7 +141,7 @@ def optimize(method, assert_convergence=ASSERT_CONV,
     print(scf.RHF(mol1).kernel() - -153.222680852335)
 
     mf = dft.RKS(mol)
-    mf.xc = 'pbe'
+    mf.xc = 'pbe,'
     mf.conv_tol = 1e-7
     mol1 = optimize(mf)
 

pyscf/grad/tdrks.py

@@ -357,7 +357,7 @@ def grad_elec(self, xy, singlet, atmlst=None):
     mol.build()
 
     mf = dft.RKS(mol)
-    mf.xc = 'LDA'
+    mf.xc = 'LDA,'
     mf.grids.prune = False
     mf.conv_tol = 1e-14
 #    mf.grids.level = 6
@@ -399,7 +399,7 @@ def grad_elec(self, xy, singlet, atmlst=None):
 
     mol.set_geom_('H 0 0 1.804; F 0 0 0', unit='B')
     mf = dft.RKS(mol)
-    mf.xc = 'lda'
+    mf.xc = 'lda,'
     mf.conv_tol = 1e-14
     mf.kernel()
     td = tddft.TDA(mf)

pyscf/grad/tduks.py

@@ -446,7 +446,7 @@ def grad_elec(self, xy, singlet, atmlst=None):
     mol.build()
 
     mf = dft.UKS(mol).set(conv_tol=1e-14)
-    mf.xc = 'LDA'
+    mf.xc = 'LDA,'
     mf.grids.prune = False
     mf.kernel()
 

pyscf/grad/test/test_rks.py

@@ -275,7 +275,7 @@ def test_get_vxc(self):
         grids1 = dft.gen_grid.Grids(mol1)
         grids1.atom_grid = (20,86)
         grids1.build(with_non0tab=False)
-        xc = 'lda'
+        xc = 'lda,'
         exc0 = dft.numint.nr_rks(mf0._numint, mol0, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_rks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
@@ -299,7 +299,7 @@ def test_get_vxc(self):
         self.assertAlmostEqual(dexc_t, exc1_approx[2], 2)
         self.assertAlmostEqual(dexc_t, exc1_full[2], 5)
 
-        xc = 'pbe'
+        xc = 'pbe,'
         exc0 = dft.numint.nr_rks(mf0._numint, mol0, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_rks(mf1._numint, mol1, grids1, xc, dm0)[1]