diff --git a/pyscf/df/grad/casdm2_util.py b/pyscf/df/grad/casdm2_util.py
index c0cdf0af16..b1ee96f008 100644
--- a/pyscf/df/grad/casdm2_util.py
+++ b/pyscf/df/grad/casdm2_util.py
@@ -24,7 +24,7 @@
 from pyscf import df
 from pyscf import mcscf
 from pyscf.ao2mo import _ao2mo
-from pyscf.grad.rhf import GradientsMixin
+from pyscf.grad.rhf import GradientsBase
 from pyscf.df.grad.rhf import _int3c_wrapper
 from pyscf.ao2mo.outcore import balance_partition
 from pyscf.ao2mo.incore import _conc_mos
@@ -107,7 +107,7 @@ def solve_df_rdm2 (mc_or_mc_grad, mo_cas=None, ci=None, casdm2=None):
 
     # Initialize mol and auxmol
     mol = mc_or_mc_grad.mol
-    if isinstance (mc_or_mc_grad, GradientsMixin):
+    if isinstance (mc_or_mc_grad, GradientsBase):
         mc = mc_or_mc_grad.base
     else:
         mc = mc_or_mc_grad
@@ -157,7 +157,7 @@ def solve_df_eri (mc_or_mc_grad, mo_cas=None, compact=True):
 
     # Initialize mol and auxmol
     mol = mc_or_mc_grad.mol
-    if isinstance (mc_or_mc_grad, GradientsMixin):
+    if isinstance (mc_or_mc_grad, GradientsBase):
         mc = mc_or_mc_grad.base
     else:
         mc = mc_or_mc_grad
@@ -213,7 +213,7 @@ def energy_elec_dferi (mc, mo_cas=None, ci=None, dfcasdm2=None, casdm2=None):
             List of energies corresponding to the dfcasdm2s,
             E = (P|ij) d_Pij / 2 = (P|ij) (P|Q)^-1 (Q|kl) d_ijkl / 2
     '''
-    if isinstance (mc, GradientsMixin): mc = mc.base
+    if isinstance (mc, GradientsBase): mc = mc.base
     if mo_cas is None:
         ncore = mc.ncore
         nocc = ncore + mc.ncas
@@ -273,7 +273,7 @@ def gfock_dferi (mc, mo_cas=None, ci=None, dfcasdm2=None, casdm2=None, max_memor
         gfock: ndarray of shape (nset, nmo[0], nmo[1]) or (nset, nao, nao)
 
     '''
-    if isinstance (mc, GradientsMixin): mc = mc.base
+    if isinstance (mc, GradientsBase): mc = mc.base
     if mo_cas is None:
         ncore = mc.ncore
         nocc = ncore + mc.ncas
@@ -328,7 +328,7 @@ def grad_elec_auxresponse_dferi (mc_grad, mo_cas=None, ci=None, dfcasdm2=None, c
     Returns:
         dE: list of ndarray of shape (len (atmlst), 3) '''
 
-    if isinstance (mc_grad, GradientsMixin):
+    if isinstance (mc_grad, GradientsBase):
         mc = mc_grad.base
     else:
         mc = mc_grad
@@ -427,7 +427,7 @@ def grad_elec_dferi (mc_grad, mo_cas=None, ci=None, dfcasdm2=None, casdm2=None,
 
     Returns:
         dE: ndarray of shape (len (dfcasdm2), len (atmlst), 3) '''
-    if isinstance (mc_grad, GradientsMixin):
+    if isinstance (mc_grad, GradientsBase):
         mc = mc_grad.base
     else:
         mc = mc_grad
diff --git a/pyscf/geomopt/addons.py b/pyscf/geomopt/addons.py
index d30f1c5c26..e4ad7b5a41 100644
--- a/pyscf/geomopt/addons.py
+++ b/pyscf/geomopt/addons.py
@@ -33,7 +33,7 @@ def as_pyscf_method(mol, scan_function):
     >>> m = as_pyscf_method(mol, scan_fn)
     >>> pyscf.geomopt.berny_solver.kernel(m)
     '''
-    from pyscf.grad.rhf import GradientsMixin
+    from pyscf.grad.rhf import GradientsBase
     class OmniGrad(lib.GradScanner):
         def __init__(self, g):
             self.__dict__.update(g.__dict__)
@@ -45,7 +45,7 @@ def __call__(self, mol):
         def converged(self):
             return True
 
-    class Gradients(GradientsMixin):
+    class Gradients(GradientsBase):
         def as_scanner(self):
             return OmniGrad(self)
 
diff --git a/pyscf/geomopt/berny_solver.py b/pyscf/geomopt/berny_solver.py
index 4e039e12f8..6159a310bc 100644
--- a/pyscf/geomopt/berny_solver.py
+++ b/pyscf/geomopt/berny_solver.py
@@ -35,7 +35,7 @@
 from pyscf.geomopt.addons import (as_pyscf_method, dump_mol_geometry,
                                   symmetrize)  # noqa
 from pyscf import __config__
-from pyscf.grad.rhf import GradientsMixin
+from pyscf.grad.rhf import GradientsBase
 
 from berny import Berny, geomlib, coords
 
@@ -119,7 +119,7 @@ def kernel(method, assert_convergence=ASSERT_CONV,
 
     if isinstance(method, lib.GradScanner):
         g_scanner = method
-    elif isinstance(method, GradientsMixin):
+    elif isinstance(method, GradientsBase):
         g_scanner = method.as_scanner()
     elif getattr(method, 'nuc_grad_method', None):
         g_scanner = method.nuc_grad_method().as_scanner()
diff --git a/pyscf/geomopt/geometric_solver.py b/pyscf/geomopt/geometric_solver.py
index f11c46f9c6..bedb0f42da 100644
--- a/pyscf/geomopt/geometric_solver.py
+++ b/pyscf/geomopt/geometric_solver.py
@@ -29,7 +29,7 @@
 from pyscf.geomopt.addons import (as_pyscf_method, dump_mol_geometry,
                                   symmetrize)  # noqa
 from pyscf import __config__
-from pyscf.grad.rhf import GradientsMixin
+from pyscf.grad.rhf import GradientsBase
 
 try:
     from geometric import internal, optimize, nifty, engine, molecule
@@ -120,7 +120,7 @@ def kernel(method, assert_convergence=ASSERT_CONV,
     '''
     if isinstance(method, lib.GradScanner):
         g_scanner = method
-    elif isinstance(method, GradientsMixin):
+    elif isinstance(method, GradientsBase):
         g_scanner = method.as_scanner()
     elif getattr(method, 'nuc_grad_method', None):
         g_scanner = method.nuc_grad_method().as_scanner()
diff --git a/pyscf/grad/casci.py b/pyscf/grad/casci.py
index 7baf057d76..144f202436 100644
--- a/pyscf/grad/casci.py
+++ b/pyscf/grad/casci.py
@@ -253,14 +253,14 @@ def __call__(self, mol_or_geom, state=None, **kwargs):
         return e_tot, de
 
 
-class Gradients(rhf_grad.GradientsMixin):
+class Gradients(rhf_grad.GradientsBase):
     '''Non-relativistic restricted Hartree-Fock gradients'''
     def __init__(self, mc):
         if isinstance(mc, StateAverageMCSCFSolver):
             self.state = None  # not a specific state
         else:
             self.state = 0  # of which the gradients to be computed.
-        rhf_grad.GradientsMixin.__init__(self, mc)
+        rhf_grad.GradientsBase.__init__(self, mc)
 
     def dump_flags(self, verbose=None):
         log = logger.new_logger(self, verbose)
diff --git a/pyscf/grad/ccsd.py b/pyscf/grad/ccsd.py
index 297d964def..2d8841c0dc 100644
--- a/pyscf/grad/ccsd.py
+++ b/pyscf/grad/ccsd.py
@@ -417,7 +417,7 @@ def _load_block_tril(h5dat, row0, row1, nao, out=None):
 def _cp(a):
     return numpy.array(a, copy=False, order='C')
 
-class Gradients(rhf_grad.GradientsMixin):
+class Gradients(rhf_grad.GradientsBase):
 
     grad_elec = grad_elec
 
diff --git a/pyscf/grad/cisd.py b/pyscf/grad/cisd.py
index 1bac19ba17..692db23ba3 100644
--- a/pyscf/grad/cisd.py
+++ b/pyscf/grad/cisd.py
@@ -134,10 +134,10 @@ def converged(self):
             ci_conv = ci_scanner.converged
         return all((ci_scanner._scf.converged, ci_conv))
 
-class Gradients(rhf_grad.GradientsMixin):
+class Gradients(rhf_grad.GradientsBase):
     def __init__(self, myci):
         self.state = 0  # of which the gradients to be computed.
-        rhf_grad.GradientsMixin.__init__(self, myci)
+        rhf_grad.GradientsBase.__init__(self, myci)
 
     def dump_flags(self, verbose=None):
         log = logger.new_logger(self, verbose)
diff --git a/pyscf/grad/dhf.py b/pyscf/grad/dhf.py
index 3f44e003a3..86d3fba309 100644
--- a/pyscf/grad/dhf.py
+++ b/pyscf/grad/dhf.py
@@ -124,7 +124,7 @@ def get_coulomb_hf(mol, dm, level='SSSS'):
 get_veff = get_coulomb_hf
 
 
-class GradientsMixin(rhf_grad.GradientsMixin):
+class GradientsBase(rhf_grad.GradientsBase):
     '''
     Basic nuclear gradient functions for 4C relativistic methods
     '''
@@ -158,13 +158,13 @@ def get_ovlp(self, mol=None):
         return get_ovlp(mol)
 
 
-class Gradients(GradientsMixin):
+class Gradients(GradientsBase):
     '''Unrestricted Dirac-Hartree-Fock gradients'''
 
     _keys = set(['level'])
 
     def __init__(self, scf_method):
-        GradientsMixin.__init__(self, scf_method)
+        GradientsBase.__init__(self, scf_method)
         if scf_method.with_ssss:
             self.level = 'SSSS'
         else:
diff --git a/pyscf/grad/lagrange.py b/pyscf/grad/lagrange.py
index eb66360323..4b31ed9bc7 100644
--- a/pyscf/grad/lagrange.py
+++ b/pyscf/grad/lagrange.py
@@ -26,7 +26,7 @@
 default_conv_rtol = getattr (__config__, 'grad_lagrange_Gradients_conv_rtol', 1e-7)
 default_max_cycle = getattr (__config__, 'grad_lagrange_Gradients_max_cycle', 50)
 
-class Gradients (rhf_grad.GradientsMixin):
+class Gradients (rhf_grad.GradientsBase):
     r''' Dummy parent class for calculating analytical nuclear gradients using the technique of
     Lagrange multipliers:
     L = E + \sum_i z_i L_i
@@ -73,7 +73,7 @@ def __init__(self, method, nlag):
         self.conv_atol = default_conv_atol
         self.conv_rtol = default_conv_rtol
         self.max_cycle = default_max_cycle
-        rhf_grad.GradientsMixin.__init__(self, method)
+        rhf_grad.GradientsBase.__init__(self, method)
 
     def debug_lagrange (self, Lvec, bvec, Aop, Adiag, **kwargs):
         logger.debug (self, "{} gradient Lagrange factor debugging not enabled".format (
diff --git a/pyscf/grad/mp2.py b/pyscf/grad/mp2.py
index 62c8c709c7..2e72892926 100644
--- a/pyscf/grad/mp2.py
+++ b/pyscf/grad/mp2.py
@@ -276,7 +276,7 @@ def _index_frozen_active(frozen_mask, mo_occ):
     VF = numpy.where((~frozen_mask) & (mo_occ==0))[0] # virtual frozen orbitals
     return OA, VA, OF, VF
 
-class Gradients(rhf_grad.GradientsMixin):
+class Gradients(rhf_grad.GradientsBase):
 
     grad_elec = grad_elec
 
diff --git a/pyscf/grad/rhf.py b/pyscf/grad/rhf.py
index a20753465c..2ec791aeee 100644
--- a/pyscf/grad/rhf.py
+++ b/pyscf/grad/rhf.py
@@ -273,7 +273,7 @@ def __call__(self, mol_or_geom, **kwargs):
         return e_tot, de
 
 
-class GradientsMixin(lib.StreamObject):
+class GradientsBase(lib.StreamObject):
     '''
     Basic nuclear gradient functions for non-relativistic methods
     '''
@@ -441,7 +441,11 @@ def _tag_rdm1 (self, dm, mo_coeff, mo_occ):
         to be split into alpha,beta in DF-ROHF subclass'''
         return lib.tag_array (dm, mo_coeff=mo_coeff, mo_occ=mo_occ)
 
-class Gradients(GradientsMixin):
+# export the symbol GradientsMixin for backward compatibility.
+# GradientsMixin should be dropped in the future.
+GradientsMixin = GradientsBase
+
+class Gradients(GradientsBase):
     '''Non-relativistic restricted Hartree-Fock gradients'''
 
     def get_veff(self, mol=None, dm=None):
diff --git a/pyscf/grad/tdrhf.py b/pyscf/grad/tdrhf.py
index 76f5b289d1..e4c73735d3 100644
--- a/pyscf/grad/tdrhf.py
+++ b/pyscf/grad/tdrhf.py
@@ -233,7 +233,7 @@ def converged(self):
                     td_scanner.converged[self.state]))
 
 
-class Gradients(rhf_grad.GradientsMixin):
+class Gradients(rhf_grad.GradientsBase):
 
     cphf_max_cycle = getattr(__config__, 'grad_tdrhf_Gradients_cphf_max_cycle', 20)
     cphf_conv_tol = getattr(__config__, 'grad_tdrhf_Gradients_cphf_conv_tol', 1e-8)
diff --git a/pyscf/grad/uhf.py b/pyscf/grad/uhf.py
index 7fed70c1c9..0f46458975 100644
--- a/pyscf/grad/uhf.py
+++ b/pyscf/grad/uhf.py
@@ -90,7 +90,7 @@ def make_rdm1e(mo_energy, mo_coeff, mo_occ):
                           rhf_grad.make_rdm1e(mo_energy[1], mo_coeff[1], mo_occ[1])))
 
 
-class Gradients(rhf_grad.GradientsMixin):
+class Gradients(rhf_grad.GradientsBase):
     '''Non-relativistic unrestricted Hartree-Fock gradients
     '''
     def get_veff(self, mol=None, dm=None):
diff --git a/pyscf/md/integrators.py b/pyscf/md/integrators.py
index e9f5ceabf2..f295256dbc 100644
--- a/pyscf/md/integrators.py
+++ b/pyscf/md/integrators.py
@@ -21,7 +21,7 @@
 from pyscf import data
 from pyscf import lib
 from pyscf.lib import logger
-from pyscf.grad.rhf import GradientsMixin
+from pyscf.grad.rhf import GradientsBase
 
 
 class Frame:
@@ -124,7 +124,7 @@ class _Integrator(lib.StreamObject):
     but inherited from.
 
     Args:
-        method : lib.GradScanner, rhf.GradientsMixin instance, or
+        method : lib.GradScanner, rhf.GradientsBase instance, or
         has nuc_grad_method method.
             Method by which to compute the energy gradients and energies
             in order to propogate the equations of motion. Realistically,
@@ -187,7 +187,7 @@ def __init__(self, method, **kwargs):
 
         if isinstance(method, lib.GradScanner):
             self.scanner = method
-        elif isinstance(method, GradientsMixin):
+        elif isinstance(method, GradientsBase):
             self.scanner = method.as_scanner()
         elif getattr(method, 'nuc_grad_method', None):
             self.scanner = method.nuc_grad_method().as_scanner()
@@ -423,7 +423,7 @@ class VelocityVerlet(_Integrator):
     '''Velocity Verlet algorithm
 
     Args:
-        method : lib.GradScanner or rhf.GradientsMixin instance, or
+        method : lib.GradScanner or rhf.GradientsBase instance, or
         has nuc_grad_method method.
             Method by which to compute the energy gradients and energies
             in order to propagate the equations of motion. Realistically,
diff --git a/pyscf/pbc/geomopt/geometric_solver.py b/pyscf/pbc/geomopt/geometric_solver.py
index 6523320258..d6f468d84a 100644
--- a/pyscf/pbc/geomopt/geometric_solver.py
+++ b/pyscf/pbc/geomopt/geometric_solver.py
@@ -25,7 +25,7 @@
 from pyscf import lib
 from pyscf.geomopt.addons import dump_mol_geometry
 from pyscf import __config__
-from pyscf.pbc.grad.krhf import GradientsMixin
+from pyscf.pbc.grad.krhf import GradientsBase
 
 try:
     from geometric import internal, optimize, nifty, engine, molecule
@@ -114,7 +114,7 @@ def kernel(method, assert_convergence=ASSERT_CONV,
     '''
     if isinstance(method, lib.GradScanner):
         g_scanner = method
-    elif isinstance(method, GradientsMixin):
+    elif isinstance(method, GradientsBase):
         g_scanner = method.as_scanner()
     elif getattr(method, 'nuc_grad_method', None):
         g_scanner = method.nuc_grad_method().as_scanner()
diff --git a/pyscf/pbc/grad/krhf.py b/pyscf/pbc/grad/krhf.py
index a7ba41a3b8..9fd628882f 100644
--- a/pyscf/pbc/grad/krhf.py
+++ b/pyscf/pbc/grad/krhf.py
@@ -283,14 +283,14 @@ def make_rdm1e(mo_energy, mo_coeff, mo_occ):
     dm1e = [molgrad.make_rdm1e(mo_energy[k], mo_coeff[k], mo_occ[k]) for k in range(nkpts)]
     return np.asarray(dm1e)
 
-class GradientsMixin(molgrad.GradientsMixin):
+class GradientsBase(molgrad.GradientsBase):
     '''
     Basic nuclear gradient functions for non-relativistic methods
     '''
     def __init__(self, method):
         self.cell = method.cell
         self.kpts = method.kpts
-        molgrad.GradientsMixin.__init__(self, method)
+        molgrad.GradientsBase.__init__(self, method)
 
     def get_hcore(self, cell=None, kpts=None):
         if cell is None: cell = self.cell
@@ -379,7 +379,7 @@ def __call__(self, cell_or_geom, **kwargs):
         return e_tot, de
 
 
-class Gradients(GradientsMixin):
+class Gradients(GradientsBase):
     '''Non-relativistic restricted Hartree-Fock gradients'''
 
     def get_veff(self, dm=None, kpts=None):
diff --git a/pyscf/x2c/sfx2c1e.py b/pyscf/x2c/sfx2c1e.py
index 61c003ce44..c4d6cf8dc8 100644
--- a/pyscf/x2c/sfx2c1e.py
+++ b/pyscf/x2c/sfx2c1e.py
@@ -155,7 +155,7 @@ def _transfer_attrs_(self, dst):
         return hf.SCF._transfer_attrs_(self, dst)
 
 
-class SpinFreeX2CHelper(x2c.X2CHelperMixin):
+class SpinFreeX2CHelper(x2c.X2CHelperBase):
     '''1-component X2c (spin-free part only)
     '''
     def get_hcore(self, mol=None):
diff --git a/pyscf/x2c/x2c.py b/pyscf/x2c/x2c.py
index 2bf017ea56..b0b8f7b8f5 100644
--- a/pyscf/x2c/x2c.py
+++ b/pyscf/x2c/x2c.py
@@ -31,7 +31,7 @@
 
 LINEAR_DEP_THRESHOLD = 1e-9
 
-class X2CHelperMixin(lib.StreamObject):
+class X2CHelperBase(lib.StreamObject):
     '''2-component X2c (including spin-free and spin-dependent terms) in
     the j-adapted spinor basis.
     '''
@@ -270,7 +270,7 @@ def reset(self, mol):
         self.mol = mol
         return self
 
-class SpinorX2CHelper(X2CHelperMixin):
+class SpinorX2CHelper(X2CHelperBase):
     '''2-component X2c (including spin-free and spin-dependent terms) in
     the j-adapted spinor basis.
     '''
@@ -278,7 +278,7 @@ class SpinorX2CHelper(X2CHelperMixin):
 
 X2C = SpinorX2CHelper
 
-class SpinOrbitalX2CHelper(X2CHelperMixin):
+class SpinOrbitalX2CHelper(X2CHelperBase):
     '''2-component X2c (including spin-free and spin-dependent terms) in
     the Gaussian type spin-orbital basis (as the spin-orbital basis in GHF)
     '''
@@ -332,7 +332,7 @@ def get_hcore(self, mol=None):
             h1 = reduce(lib.dot, (contr_coeff.T, h1, contr_coeff))
         return h1
 
-    @lib.with_doc(X2CHelperMixin.picture_change.__doc__)
+    @lib.with_doc(X2CHelperBase.picture_change.__doc__)
     def picture_change(self, even_operator=(None, None), odd_operator=None):
         mol = self.mol
         xmol, c = self.get_xmol(mol)