diff --git a/pyscf/agf2/ragf2.py b/pyscf/agf2/ragf2.py
index 49a1939577..99c7c07102 100644
--- a/pyscf/agf2/ragf2.py
+++ b/pyscf/agf2/ragf2.py
@@ -869,8 +869,7 @@ def density_fit(self, auxbasis=None, with_df=None):
             myagf2.with_df = with_df
 
         if auxbasis is not None and myagf2.with_df.auxbasis != auxbasis:
-            import copy
-            myagf2.with_df = copy.copy(myagf2.with_df)
+            myagf2.with_df = myagf2.with_df.copy()
             myagf2.with_df.auxbasis = auxbasis
 
         return myagf2
diff --git a/pyscf/agf2/uagf2.py b/pyscf/agf2/uagf2.py
index 6c6778c290..e08cab67c8 100644
--- a/pyscf/agf2/uagf2.py
+++ b/pyscf/agf2/uagf2.py
@@ -633,8 +633,7 @@ def density_fit(self, auxbasis=None, with_df=None):
             myagf2.with_df = with_df
 
         if auxbasis is not None and myagf2.with_df.auxbasis != auxbasis:
-            import copy
-            myagf2.with_df = copy.copy(myagf2.with_df)
+            myagf2.with_df = myagf2.with_df.copy()
             myagf2.with_df.auxbasis = auxbasis
 
         return myagf2
diff --git a/pyscf/cc/ccsd.py b/pyscf/cc/ccsd.py
index 30284b51e6..91e790b9c6 100644
--- a/pyscf/cc/ccsd.py
+++ b/pyscf/cc/ccsd.py
@@ -1355,8 +1355,7 @@ def density_fit(self, auxbasis=None, with_df=None):
         if with_df is not None:
             mycc.with_df = with_df
         if mycc.with_df.auxbasis != auxbasis:
-            import copy
-            mycc.with_df = copy.copy(mycc.with_df)
+            mycc.with_df = mycc.with_df.copy()
             mycc.with_df.auxbasis = auxbasis
         return mycc
 
diff --git a/pyscf/cc/test/test_ccsd_lambda.py b/pyscf/cc/test/test_ccsd_lambda.py
index 05a57fc40b..3bbc7d0ebd 100644
--- a/pyscf/cc/test/test_ccsd_lambda.py
+++ b/pyscf/cc/test/test_ccsd_lambda.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import tempfile
 import unittest
 import numpy
@@ -114,7 +113,7 @@ def test_ccsd(self):
 
     def test_restart(self):
         ftmp = tempfile.NamedTemporaryFile()
-        cc1 = copy.copy(mycc)
+        cc1 = mycc.copy()
         cc1.max_cycle = 5
         cc1.solve_lambda()
         l1ref = cc1.l1
diff --git a/pyscf/cc/test/test_dfccsd.py b/pyscf/cc/test/test_dfccsd.py
index f1aa92486d..a887020ca4 100644
--- a/pyscf/cc/test/test_dfccsd.py
+++ b/pyscf/cc/test/test_dfccsd.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import unittest
-import copy
 import numpy
 from functools import reduce
 
@@ -24,7 +23,7 @@
 from pyscf.cc import dfccsd, eom_rccsd
 
 def make_mycc1():
-    mf1 = copy.copy(mf)
+    mf1 = mf.copy()
     no = mol.nelectron // 2
     n = mol.nao_nr()
     nv = n - no
diff --git a/pyscf/cc/test/test_eom_gccsd.py b/pyscf/cc/test/test_eom_gccsd.py
index 2937f71f88..05fa173db2 100644
--- a/pyscf/cc/test/test_eom_gccsd.py
+++ b/pyscf/cc/test/test_eom_gccsd.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import unittest
-import copy
 import numpy
 from functools import reduce
 
diff --git a/pyscf/cc/test/test_eom_rccsd.py b/pyscf/cc/test/test_eom_rccsd.py
index b1210ac7b3..8b8f867c11 100644
--- a/pyscf/cc/test/test_eom_rccsd.py
+++ b/pyscf/cc/test/test_eom_rccsd.py
@@ -26,7 +26,7 @@
 from pyscf.cc import ccsd, rccsd, eom_rccsd, rintermediates, gintermediates
 
 def make_mycc1():
-    mf1 = copy.copy(mf)
+    mf1 = mf.copy()
     no = mol.nelectron // 2
     n = mol.nao_nr()
     nv = n - no
@@ -67,7 +67,7 @@ def setUpModule():
     mf1, mycc1, eris1 = make_mycc1()
     no, nv = mycc1.t1.shape
 
-    mycci = copy.copy(mycc1)
+    mycci = mycc1.copy()
     erisi = copy.copy(eris1)
     erisi.oooo = eris1.oooo + numpy.sin(eris1.oooo)*1j
     erisi.oooo = erisi.oooo + erisi.oooo.conj().transpose(1,0,3,2)
diff --git a/pyscf/cc/test/test_eom_uccsd.py b/pyscf/cc/test/test_eom_uccsd.py
index 754eeb873e..d4b551df20 100644
--- a/pyscf/cc/test/test_eom_uccsd.py
+++ b/pyscf/cc/test/test_eom_uccsd.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 import numpy
 import numpy as np
@@ -47,7 +46,7 @@ def setUpModule():
     mol1.spin = 2
     mol1.build()
     mf0 = scf.UHF(mol1).run(conv_tol=1e-12)
-    mf1 = copy.copy(mf0)
+    mf1 = mf0.copy()
 
     nocca, noccb = mol1.nelec
     nmo = mol1.nao_nr()
@@ -218,7 +217,7 @@ def test_ucc_eomee_ccsd_matvec(self):
         self.assertAlmostEqual(float(abs(gr2-uee1.spatial2spin(r2, orbspin)).max()), 0, 9)
         self.assertAlmostEqual(lib.fp(vec1), 49.499911123484523, 9)
 
-        ucc2 = copy.copy(ucc1)
+        ucc2 = ucc1.copy()
         ucc2.direct = True
         uee1 = eom_uccsd.EOMEESpinKeep(ucc2)
         vec1 = uee1.matvec(vec)
@@ -293,7 +292,7 @@ def test_ucc_eomsf_ccsd_matvec(self):
         self.assertAlmostEqual(float(abs(gr1-myeom.spatial2spin(v1, orbspin)).max()), 0, 9)
         self.assertAlmostEqual(float(abs(gr2-myeom.spatial2spin(v2, orbspin)).max()), 0, 9)
 
-        ucc2 = copy.copy(ucc1)
+        ucc2 = ucc1.copy()
         ucc2.direct = True
         myeom = eom_uccsd.EOMEESpinFlip(ucc2)
         vec1 = myeom.matvec(vec)
@@ -322,7 +321,7 @@ def test_ucc_eomip_matvec(self):
         self.assertAlmostEqual(float(abs(gr1-myeom.spatial2spin(v1, orbspin)).max()), 0, 9)
         self.assertAlmostEqual(float(abs(gr2-myeom.spatial2spin(v2, orbspin)).max()), 0, 9)
 
-        ucc2 = copy.copy(ucc1)
+        ucc2 = ucc1.copy()
         ucc2.direct = True
         myeom = eom_uccsd.EOMIP(ucc2)
         vec1 = myeom.matvec(vec)
@@ -351,7 +350,7 @@ def test_ucc_eomea_matvec(self):
         self.assertAlmostEqual(float(abs(gr1-myeom.spatial2spin(v1, orbspin)).max()), 0, 9)
         self.assertAlmostEqual(float(abs(gr2-myeom.spatial2spin(v2, orbspin)).max()), 0, 9)
 
-        ucc2 = copy.copy(ucc1)
+        ucc2 = ucc1.copy()
         ucc2.direct = True
         myeom = eom_uccsd.EOMEA(ucc2)
         vec1 = myeom.matvec(vec)
diff --git a/pyscf/cc/test/test_gccsd_t.py b/pyscf/cc/test/test_gccsd_t.py
index 86efcb1220..e448cb3d3e 100644
--- a/pyscf/cc/test/test_gccsd_t.py
+++ b/pyscf/cc/test/test_gccsd_t.py
@@ -15,7 +15,6 @@
 
 import unittest
 import numpy
-import copy
 from functools import reduce
 
 from pyscf import gto, scf, lib, symm
@@ -36,7 +35,7 @@ def setUpModule():
     mol.basis = '3-21g'
     mol.symmetry = 'C2v'
     mol.build()
-    mol1 = copy.copy(mol)
+    mol1 = mol.copy()
     mol1.symmetry = False
 
     mf = scf.UHF(mol1).run(conv_tol=1e-14)
@@ -53,7 +52,7 @@ def test_gccsd_t_compare_uccsd_t(self):
         self.assertAlmostEqual(myucc.ccsd_t(), mygcc.ccsd_t(t1=None), 7)
 
     def test_gccsd_t(self):
-        mf1 = copy.copy(mf)
+        mf1 = mf.copy()
         nao, nmo = mf.mo_coeff[0].shape
         numpy.random.seed(10)
         mf1.mo_coeff = numpy.random.random((2,nao,nmo))
diff --git a/pyscf/cc/test/test_rccsd.py b/pyscf/cc/test/test_rccsd.py
index 6c1ebe3037..8f5687efe7 100644
--- a/pyscf/cc/test/test_rccsd.py
+++ b/pyscf/cc/test/test_rccsd.py
@@ -112,7 +112,7 @@ def ao2mofn(mos):
         self.assertAlmostEqual(abs(eris1.vvvv.imag-eris.vvvv).max(), 0, 11)
 
     def test_dump_chk(self):
-        cc1 = copy.copy(mycc)
+        cc1 = mycc.copy()
         cc1.nmo = mf.mo_energy.size
         cc1.nocc = mol.nelectron // 2
         cc1.dump_chk()
@@ -221,7 +221,7 @@ def test_vector_size(self):
         self.assertEqual(mycc.vector_size(), 860)
 
     def test_rccsd_frozen(self):
-        cc1 = copy.copy(mycc)
+        cc1 = mycc.copy()
         cc1.frozen = 1
         self.assertEqual(cc1.nmo, 12)
         self.assertEqual(cc1.nocc, 4)
@@ -383,7 +383,7 @@ def test_add_vvvv(self):
         nocc, nvir = t1.shape
         tau = t2 + numpy.einsum('ia,jb->ijab', t1, t1)
         eris1 = copy.copy(eris)
-        mycc1 = copy.copy(mycc)
+        mycc1 = mycc.copy()
         ovvv = eris1.get_ovvv()
         tmp = -numpy.einsum('ijcd,ka,kdcb->ijba', tau, t1, ovvv)
         t2a = tmp + tmp.transpose(1,0,3,2)
diff --git a/pyscf/cc/test/test_rccsd_lambda.py b/pyscf/cc/test/test_rccsd_lambda.py
index 17824baa16..ae46eb507e 100644
--- a/pyscf/cc/test/test_rccsd_lambda.py
+++ b/pyscf/cc/test/test_rccsd_lambda.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 import numpy
 import numpy as np
@@ -170,7 +169,7 @@ def test_rdm(self):
         self.assertAlmostEqual(e1, mycc.e_tot, 6)
 
         d1 = ccsd_rdm._gamma1_intermediates(mycc, mycc.t1, mycc.t2, mycc.l1, mycc.l2)
-        mycc1 = copy.copy(mycc)
+        mycc1 = mycc.copy()
         mycc1.max_memory = 0
         d2 = ccsd_rdm._gamma2_intermediates(mycc1, mycc.t1, mycc.t2, mycc.l1, mycc.l2, True)
         dm2 = ccsd_rdm._make_rdm2(mycc, d1, d2, with_dm1=True, with_frozen=True)
diff --git a/pyscf/cc/test/test_uccsd.py b/pyscf/cc/test/test_uccsd.py
index 52a26c2c45..a05f3bb4de 100644
--- a/pyscf/cc/test/test_uccsd.py
+++ b/pyscf/cc/test/test_uccsd.py
@@ -187,7 +187,7 @@ def test_amplitudes_from_rccsd(self):
         self.assertAlmostEqual(abs(t2[2]-myucc.t2[2]).max(), 0, 5)
 
     def test_uccsd_frozen(self):
-        ucc1 = copy.copy(myucc)
+        ucc1 = myucc.copy()
         ucc1.frozen = 1
         self.assertEqual(ucc1.nmo, (12,12))
         self.assertEqual(ucc1.nocc, (4,4))
diff --git a/pyscf/cc/test/test_uccsd_t.py b/pyscf/cc/test/test_uccsd_t.py
index b7b7452c1c..2c39758dd0 100644
--- a/pyscf/cc/test/test_uccsd_t.py
+++ b/pyscf/cc/test/test_uccsd_t.py
@@ -15,7 +15,6 @@
 
 import unittest
 import numpy
-import copy
 from functools import reduce
 
 from pyscf import gto, scf, lib, symm
@@ -43,7 +42,7 @@ def setUpModule():
     mol.basis = '3-21g'
     mol.symmetry = 'C2v'
     mol.build()
-    mol1 = copy.copy(mol)
+    mol1 = mol.copy()
     mol1.symmetry = False
 
     mf = scf.UHF(mol1).run(conv_tol=1e-14)
@@ -58,7 +57,7 @@ def tearDownModule():
 
 class KnownValues(unittest.TestCase):
     def test_uccsd_t(self):
-        mf1 = copy.copy(mf)
+        mf1 = mf.copy()
         nao, nmo = mf.mo_coeff[0].shape
         numpy.random.seed(10)
         mf1.mo_coeff = numpy.random.random((2,nao,nmo)) - .5
diff --git a/pyscf/df/addons.py b/pyscf/df/addons.py
index b042963b3a..7166677203 100644
--- a/pyscf/df/addons.py
+++ b/pyscf/df/addons.py
@@ -17,7 +17,6 @@
 #
 
 import sys
-import copy
 import numpy
 from pyscf.lib import logger
 from pyscf import gto
@@ -204,7 +203,7 @@ def make_auxmol(mol, auxbasis=None):
 
     See also the paper JCTC, 13, 554 about generating auxiliary fitting basis.
     '''
-    pmol = copy.copy(mol)  # just need shallow copy
+    pmol = mol.copy(deep=False)
 
     if auxbasis is None:
         auxbasis = make_auxbasis(mol)
diff --git a/pyscf/df/df.py b/pyscf/df/df.py
index 7261283ed8..42928454bd 100644
--- a/pyscf/df/df.py
+++ b/pyscf/df/df.py
@@ -21,7 +21,6 @@
 '''
 
 
-import copy
 import tempfile
 import contextlib
 import numpy
@@ -277,7 +276,7 @@ def range_coulomb(self, omega):
         if key in self._rsh_df:
             rsh_df = self._rsh_df[key]
         else:
-            rsh_df = self._rsh_df[key] = copy.copy(self).reset()
+            rsh_df = self._rsh_df[key] = self.copy().reset()
             if hasattr(self, '_dataname'):
                 rsh_df._dataname = f'{self._dataname}/lr/{key}'
             logger.info(self, 'Create RSH-DF object %s for omega=%s', rsh_df, omega)
diff --git a/pyscf/df/df_jk.py b/pyscf/df/df_jk.py
index 9a6d03dbf2..ee1e27a211 100644
--- a/pyscf/df/df_jk.py
+++ b/pyscf/df/df_jk.py
@@ -16,8 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
-import copy
-
 import ctypes
 import numpy
 import scipy.linalg
@@ -83,7 +81,7 @@ def density_fit(mf, auxbasis=None, with_df=None, only_dfj=False):
             mf.with_df = with_df
         elif getattr(mf.with_df, 'auxbasis', None) != auxbasis:
             #logger.warn(mf, 'DF might have been initialized twice.')
-            mf = copy.copy(mf)
+            mf = mf.copy()
             mf.with_df = with_df
             mf.only_dfj = only_dfj
         return mf
diff --git a/pyscf/df/grad/sacasscf.py b/pyscf/df/grad/sacasscf.py
index 247a5de055..cef6373d2c 100644
--- a/pyscf/df/grad/sacasscf.py
+++ b/pyscf/df/grad/sacasscf.py
@@ -17,7 +17,6 @@
 #
 
 import numpy as np
-import copy
 import time
 import gc
 from functools import reduce
diff --git a/pyscf/dft/dks.py b/pyscf/dft/dks.py
index 59e170e0a3..783c3c546b 100644
--- a/pyscf/dft/dks.py
+++ b/pyscf/dft/dks.py
@@ -21,7 +21,6 @@
 '''
 
 
-import copy
 import numpy
 from pyscf import lib
 from pyscf.lib import logger
@@ -118,7 +117,7 @@ def to_dhf(self):
 
     def to_dks(self, xc=None):
         if xc is not None and xc != self.xc:
-            mf = copy.copy(self)
+            mf = self.copy()
             mf.xc = xc
             mf.converged = False
         return self
diff --git a/pyscf/dft/libxc.py b/pyscf/dft/libxc.py
index de025b16bf..4ac047f4c8 100644
--- a/pyscf/dft/libxc.py
+++ b/pyscf/dft/libxc.py
@@ -24,7 +24,6 @@
 
 import sys
 import warnings
-import copy
 import ctypes
 import math
 import numpy
@@ -1661,5 +1660,5 @@ def define_xc_(ni, description, xctype='LDA', hyb=0, rsh=(0,0,0)):
     return ni
 
 def define_xc(ni, description, xctype='LDA', hyb=0, rsh=(0,0,0)):
-    return define_xc_(copy.copy(ni), description, xctype, hyb, rsh)
+    return define_xc_(ni.copy(), description, xctype, hyb, rsh)
 define_xc.__doc__ = define_xc_.__doc__
diff --git a/pyscf/dft/xcfun.py b/pyscf/dft/xcfun.py
index 3001df721d..a5d88e7e26 100644
--- a/pyscf/dft/xcfun.py
+++ b/pyscf/dft/xcfun.py
@@ -22,7 +22,6 @@
 U. Ekstrom et al, J. Chem. Theory Comput., 6, 1971
 '''
 
-import copy
 import ctypes
 from functools import lru_cache
 import math
@@ -1080,5 +1079,5 @@ def define_xc_(ni, description, xctype='LDA', hyb=0, rsh=(0,0,0)):
     return ni
 
 def define_xc(ni, description, xctype='LDA', hyb=0, rsh=(0,0,0)):
-    return define_xc_(copy.copy(ni), description, xctype, hyb, rsh)
+    return define_xc_(ni.copy(), description, xctype, hyb, rsh)
 define_xc.__doc__ = define_xc_.__doc__
diff --git a/pyscf/eph/eph_fd.py b/pyscf/eph/eph_fd.py
index aa89e51b66..f0fa4bb421 100644
--- a/pyscf/eph/eph_fd.py
+++ b/pyscf/eph/eph_fd.py
@@ -20,7 +20,6 @@
 electron-phonon matrix from finite difference
 '''
 
-import copy
 import numpy as np
 from pyscf import scf, dft, gto, hessian
 from pyscf.eph import rhf as rhf_eph
@@ -35,9 +34,9 @@ def run_mfs(mf, mols_a, mols_b):
     dm0 = mf.make_rdm1()
     mflist = []
     for i in range(nconfigs):
-        mf1 = copy.copy(mf)
+        mf1 = mf.copy()
         mf1.reset(mols_a[i])
-        mf2 = copy.copy(mf)
+        mf2 = mf.copy()
         mf2.reset(mols_b[i])
         mf1.kernel(dm0=dm0)
         mf2.kernel(dm0=dm0)
diff --git a/pyscf/fci/addons.py b/pyscf/fci/addons.py
index f48143a035..caa549b6a0 100644
--- a/pyscf/fci/addons.py
+++ b/pyscf/fci/addons.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import sys
-import copy
 import warnings
 import numpy
 from pyscf import lib
@@ -604,7 +603,7 @@ class FCISolver (fciobj_class, SpinPenaltyFCISolver):
         _keys = set(('ss_value', 'ss_penalty', 'base'))
 
         def __init__(self, fcibase):
-            self.base = copy.copy (fcibase)
+            self.base = fcibase.copy()
             self.__dict__.update (fcibase.__dict__)
             self.ss_value = ss_value
             self.ss_penalty = shift
@@ -619,7 +618,7 @@ def base_contract_2e (self, *args, **kwargs):
     fciobj.__dict__.update (new_fciobj.__dict__)
     return fciobj
 def fix_spin(fciobj, shift=.1, ss=None):
-    return fix_spin_(copy.copy(fciobj), shift, ss)
+    return fix_spin_(fciobj.copy(), shift, ss)
 
 def transform_ci_for_orbital_rotation(ci, norb, nelec, u):
     '''
diff --git a/pyscf/grad/sacasscf.py b/pyscf/grad/sacasscf.py
index 0cb0d758f0..13f086713b 100644
--- a/pyscf/grad/sacasscf.py
+++ b/pyscf/grad/sacasscf.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 #
 
-import copy
 import gc
 import numpy as np
 from functools import reduce
@@ -505,7 +504,7 @@ def make_fcasscf (self, state=None, casscf_attr={}, fcisolver_attr={}):
         # Spin penalty method is inapplicable to response calc'ns
         # It must be deactivated for Lagrange multipliers to converge
         if isinstance (fcasscf.fcisolver, SpinPenaltyFCISolver):
-            fcasscf.fcisolver = copy.copy (fcasscf.fcisolver)
+            fcasscf.fcisolver = fcasscf.fcisolver.copy()
             fcasscf.fcisolver.ss_penalty = 0
         fcasscf.__dict__.update (casscf_attr)
         fcasscf.nelecas = nelecas
@@ -536,8 +535,7 @@ def make_fcasscf_sa (self, casscf_attr={}, fcisolver_attr={}):
         fcasscf.__dict__.update (self.base.__dict__)
         if isinstance (self.base, StateAverageMCSCFSolver):
             if isinstance (self.base.fcisolver, StateAverageMixFCISolver):
-                fcisolvers = [copy.copy (f) for f in
-                              self.base.fcisolver.fcisolvers]
+                fcisolvers = [f.copy() for f in self.base.fcisolver.fcisolvers]
                 # Spin penalty method is inapplicable to response calc'ns
                 # It must be deactivated for Lagrange multipliers to converge
                 for i in range (len (fcisolvers)):
@@ -550,7 +548,7 @@ def make_fcasscf_sa (self, casscf_attr={}, fcisolver_attr={}):
         # Spin penalty method is inapplicable to response calc'ns
         # It must be deactivated for Lagrange multipliers to converge
         if isinstance (fcasscf.fcisolver, SpinPenaltyFCISolver):
-            fcasscf.fcisolver = copy.copy (fcasscf.fcisolver)
+            fcasscf.fcisolver = fcasscf.fcisolver.copy()
             fcasscf.fcisolver.ss_penalty = 0
         fcasscf.__dict__.update (casscf_attr)
         fcasscf.fcisolver.__dict__.update (fcisolver_attr)
diff --git a/pyscf/grad/test/test_rks.py b/pyscf/grad/test/test_rks.py
index 777626a8c4..fe760651de 100644
--- a/pyscf/grad/test/test_rks.py
+++ b/pyscf/grad/test/test_rks.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import unittest
-import copy
 import numpy
 from pyscf import gto, dft, lib
 from pyscf.dft import radi
@@ -411,9 +410,9 @@ def test_get_vxc(self):
         exc0 = dft.numint.nr_rks(mf0._numint, mol0, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_rks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
-        grids0_w = copy.copy(grids0)
+        grids0_w = grids0.copy()
         grids0_w.weights = grids1.weights
-        grids0_c = copy.copy(grids0)
+        grids0_c = grids0.copy()
         grids0_c.coords = grids1.coords
         exc0_w = dft.numint.nr_rks(mf0._numint, mol0, grids0_w, xc, dm0)[1]
         exc0_c = dft.numint.nr_rks(mf1._numint, mol1, grids0_c, xc, dm0)[1]
@@ -435,9 +434,9 @@ def test_get_vxc(self):
         exc0 = dft.numint.nr_rks(mf0._numint, mol0, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_rks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
-        grids0_w = copy.copy(grids0)
+        grids0_w = grids0.copy()
         grids0_w.weights = grids1.weights
-        grids0_c = copy.copy(grids0)
+        grids0_c = grids0.copy()
         grids0_c.coords = grids1.coords
         exc0_w = dft.numint.nr_rks(mf0._numint, mol0, grids0_w, xc, dm0)[1]
         exc0_c = dft.numint.nr_rks(mf1._numint, mol1, grids0_c, xc, dm0)[1]
@@ -462,9 +461,9 @@ def test_get_vxc(self):
         exc0 = dft.numint.nr_rks(mf0._numint, mol0, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_rks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
-        grids0_w = copy.copy(grids0)
+        grids0_w = grids0.copy()
         grids0_w.weights = grids1.weights
-        grids0_c = copy.copy(grids0)
+        grids0_c = grids0.copy()
         grids0_c.coords = grids1.coords
         exc0_w = dft.numint.nr_rks(mf0._numint, mol0, grids0_w, xc, dm0)[1]
         exc0_c = dft.numint.nr_rks(mf1._numint, mol1, grids0_c, xc, dm0)[1]
@@ -486,9 +485,9 @@ def test_get_vxc(self):
         exc0 = dft.numint.nr_rks(mf0._numint, mol0, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_rks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
-        grids0_w = copy.copy(grids0)
+        grids0_w = grids0.copy()
         grids0_w.weights = grids1.weights
-        grids0_c = copy.copy(grids0)
+        grids0_c = grids0.copy()
         grids0_c.coords = grids1.coords
         exc0_w = dft.numint.nr_rks(mf0._numint, mol0, grids0_w, xc, dm0)[1]
         exc0_c = dft.numint.nr_rks(mf1._numint, mol1, grids0_c, xc, dm0)[1]
diff --git a/pyscf/grad/test/test_uks.py b/pyscf/grad/test/test_uks.py
index 653ae22088..effd7218a7 100644
--- a/pyscf/grad/test/test_uks.py
+++ b/pyscf/grad/test/test_uks.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import unittest
-import copy
 import numpy
 from pyscf import gto, dft, lib
 from pyscf.dft import radi
@@ -208,9 +207,9 @@ def test_get_vxc(self):
         exc0 = dft.numint.nr_uks(mf._numint, mol, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_uks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
-        grids0_w = copy.copy(grids0)
+        grids0_w = grids0.copy()
         grids0_w.weights = grids1.weights
-        grids0_c = copy.copy(grids0)
+        grids0_c = grids0.copy()
         grids0_c.coords = grids1.coords
         exc0_w = dft.numint.nr_uks(mf._numint, mol, grids0_w, xc, dm0)[1]
         exc0_c = dft.numint.nr_uks(mf._numint, mol1, grids0_c, xc, dm0)[1]
@@ -232,9 +231,9 @@ def test_get_vxc(self):
         exc0 = dft.numint.nr_uks(mf._numint, mol, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_uks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
-        grids0_w = copy.copy(grids0)
+        grids0_w = grids0.copy()
         grids0_w.weights = grids1.weights
-        grids0_c = copy.copy(grids0)
+        grids0_c = grids0.copy()
         grids0_c.coords = grids1.coords
         exc0_w = dft.numint.nr_uks(mf._numint, mol, grids0_w, xc, dm0)[1]
         exc0_c = dft.numint.nr_uks(mf._numint, mol1, grids0_c, xc, dm0)[1]
@@ -256,9 +255,9 @@ def test_get_vxc(self):
         exc0 = dft.numint.nr_uks(mf._numint, mol, grids0, xc, dm0)[1]
         exc1 = dft.numint.nr_uks(mf1._numint, mol1, grids1, xc, dm0)[1]
 
-        grids0_w = copy.copy(grids0)
+        grids0_w = grids0.copy()
         grids0_w.weights = grids1.weights
-        grids0_c = copy.copy(grids0)
+        grids0_c = grids0.copy()
         grids0_c.coords = grids1.coords
         exc0_w = dft.numint.nr_uks(mf._numint, mol, grids0_w, xc, dm0)[1]
         exc0_c = dft.numint.nr_uks(mf._numint, mol1, grids0_c, xc, dm0)[1]
diff --git a/pyscf/gto/mole.py b/pyscf/gto/mole.py
index 67a645cf41..092f8214ac 100644
--- a/pyscf/gto/mole.py
+++ b/pyscf/gto/mole.py
@@ -602,8 +602,7 @@ def decontract_basis(mol, atoms=None, to_cart=False):
     >>> abs(s-mol.intor('int1e_ovlp')).max()
     0.0
     '''
-    import copy
-    pmol = copy.copy(mol)
+    pmol = mol.copy(deep=False)
 
     # Some input basis may be segmented basis from a general contracted set.
     # This may lead to duplicated pGTOs. First contract all basis to remove
@@ -1187,7 +1186,7 @@ def tot_electrons(mol):
                     nelectron, nelectron_int)
     return nelectron_int
 
-def copy(mol):
+def copy(mol, deep=True):
     '''Deepcopy of the given :class:`Mole` object
 
     Some attributes are shared between the original and copied objects.
@@ -1195,7 +1194,12 @@ def copy(mol):
     not affect the original object.
     '''
     import copy
-    newmol = copy.copy(mol)
+    # Avoid copy.copy(mol) for shallow copy because copy.copy automatically
+    # calls __copy__, __reduce__, __getstate__, __setstate__ methods
+    newmol = mol.view(mol.__class__)
+    if not deep:
+        return newmol
+
     newmol._atm    = numpy.copy(mol._atm)
     newmol._bas    = numpy.copy(mol._bas)
     newmol._env    = numpy.copy(mol._env)
@@ -2541,11 +2545,10 @@ def set_geom_(self, atoms_or_coords, unit=None, symmetry=None,
                   inplace=True):
         '''Update geometry
         '''
-        import copy
         if inplace:
             mol = self
         else:
-            mol = copy.copy(self)
+            mol = self.copy(deep=False)
             mol._env = mol._env.copy()
         if unit is None:
             unit = mol.unit
@@ -3170,7 +3173,6 @@ def __init__(self):
         # _TemporaryMoleContext which is protected by the _ctx_lock.
         self._ctx_lock = None
 
-    # NOTE: do not overwrite __copy__. It causes recursive calls in copy.copy(mol)
     copy = copy
 
     pack = pack
diff --git a/pyscf/gto/test/test_mole.py b/pyscf/gto/test/test_mole.py
index 294233e168..4adaa7e32a 100644
--- a/pyscf/gto/test/test_mole.py
+++ b/pyscf/gto/test/test_mole.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 import tempfile
 from functools import reduce
@@ -555,13 +554,13 @@ def test_nelectron(self):
         self.assertRaises(RuntimeError, lambda *args: mol0.nelec)
         mol0.spin = 1
 
-        mol1 = copy.copy(mol0)
+        mol1 = mol0.copy()
         self.assertEqual(mol1.nelec, (5, 4))
         mol1.nelec = (3, 6)
         self.assertEqual(mol1.nelec, (3, 6))
 
     def test_multiplicity(self):
-        mol1 = copy.copy(mol0)
+        mol1 = mol0.copy()
         self.assertEqual(mol1.multiplicity, 2)
         mol1.multiplicity = 5
         self.assertEqual(mol1.multiplicity, 5)
@@ -569,7 +568,7 @@ def test_multiplicity(self):
         self.assertRaises(RuntimeError, lambda:mol1.nelec)
 
     def test_ms(self):
-        mol1 = copy.copy(mol0)
+        mol1 = mol0.copy()
         self.assertEqual(mol1.ms, 0.5)
         mol1.ms = 1
         self.assertEqual(mol1.multiplicity, 3)
diff --git a/pyscf/lib/misc.py b/pyscf/lib/misc.py
index 9faf46f176..bcf5fa353c 100644
--- a/pyscf/lib/misc.py
+++ b/pyscf/lib/misc.py
@@ -666,6 +666,10 @@ def dir_attributes(cls):
         return {key for key, val in vars(cls).items()
                 if not (key.startswith('__') or callable(val))}
 
+    def copy(self):
+        '''Returns a shallow copy'''
+        return self.view(self.__class__)
+
 _warn_once_registry = {}
 def check_sanity(obj, keysref, stdout=sys.stdout):
     '''Check misinput of class attributes, check whether a class method is
diff --git a/pyscf/mcscf/df.py b/pyscf/mcscf/df.py
index 6522ac02e9..3451cce414 100644
--- a/pyscf/mcscf/df.py
+++ b/pyscf/mcscf/df.py
@@ -17,7 +17,6 @@
 #
 
 
-import copy
 import ctypes
 from functools import reduce
 import numpy
@@ -72,7 +71,7 @@ def density_fit(casscf, auxbasis=None, with_df=None):
             casscf.with_df = with_df
         elif getattr(casscf.with_df, 'auxbasis', None) != auxbasis:
             #logger.warn(casscf, 'DF might have been initialized twice.')
-            casscf = copy.copy(casscf)
+            casscf = casscf.copy()
             casscf.with_df = with_df
         return casscf
 
diff --git a/pyscf/mcscf/mc1step.py b/pyscf/mcscf/mc1step.py
index 5959c974d1..9ea4e58f70 100644
--- a/pyscf/mcscf/mc1step.py
+++ b/pyscf/mcscf/mc1step.py
@@ -18,7 +18,6 @@
 
 import sys
 
-import copy
 from functools import reduce
 import numpy
 import scipy.linalg
@@ -429,8 +428,6 @@ def kernel(casscf, mo_coeff, tol=1e-7, conv_tol_grad=None,
 
         eris = None
         # keep u, g_orb in locals() so that they can be accessed by callback
-        u = u.copy()
-        g_orb = g_orb.copy()
         mo = casscf.rotate_mo(mo, u, log)
         eris = casscf.ao2mo(mo)
         t2m = log.timer('update eri', *t3m)
diff --git a/pyscf/mcscf/newton_casscf.py b/pyscf/mcscf/newton_casscf.py
index a2aecb007b..d61a3e2ea8 100644
--- a/pyscf/mcscf/newton_casscf.py
+++ b/pyscf/mcscf/newton_casscf.py
@@ -21,7 +21,6 @@
 '''
 
 
-import copy
 from functools import reduce
 from itertools import product
 import numpy
diff --git a/pyscf/mcscf/test/test_addons.py b/pyscf/mcscf/test/test_addons.py
index 145019ecf8..06773e245a 100644
--- a/pyscf/mcscf/test/test_addons.py
+++ b/pyscf/mcscf/test/test_addons.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 from functools import reduce
 import numpy
@@ -118,7 +117,7 @@ def test_canonicalize1(self):
         self.assertAlmostEqual(e1, 44.2658681077, 7)
         self.assertAlmostEqual(lib.fp(mo_e), 4.1206025804989173, 4)
 
-        mcr1 = copy.copy(mcr)
+        mcr1 = mcr.copy()
         mcr1.frozen = 2
         mo, ci, mo_e = mcr1.canonicalize(mo1)
         self.assertAlmostEqual(lib.fp(mo_e), 6.6030999409178577, 5)
diff --git a/pyscf/mcscf/test/test_casci.py b/pyscf/mcscf/test/test_casci.py
index 7f80a11d97..121c5d3056 100644
--- a/pyscf/mcscf/test/test_casci.py
+++ b/pyscf/mcscf/test/test_casci.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 import numpy
 from pyscf import lib
@@ -165,7 +164,7 @@ def test_with_ci_init_guess(self):
         self.assertAlmostEqual(mc2.e_tot, -108.62009625745821, 7)
 
     def test_slight_symmetry_broken(self):
-        mf = copy.copy(msym)
+        mf = msym.copy()
         mf.mo_coeff = numpy.array(msym.mo_coeff)
         u = numpy.linalg.svd(numpy.random.random((4,4)))[0]
         mf.mo_coeff[:,5:9] = mf.mo_coeff[:,5:9].dot(u)
@@ -203,7 +202,7 @@ def test_state_average(self):
 
     def test_state_average_mix(self):
         mc = mcscf.CASCI(m, 4, 4)
-        cis1 = copy.copy(mc.fcisolver)
+        cis1 = mc.fcisolver.copy()
         cis1.spin = 2
         cis1.nroots = 3
         mc = mcscf.addons.state_average_mix(mc, [cis1, mc.fcisolver], [.25, .25, .25, .25])
diff --git a/pyscf/mcscf/test/test_mc1step.py b/pyscf/mcscf/test/test_mc1step.py
index 838a9ac31f..f94bb57f4d 100644
--- a/pyscf/mcscf/test/test_mc1step.py
+++ b/pyscf/mcscf/test/test_mc1step.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 import tempfile
 import numpy
@@ -307,7 +306,7 @@ def test_state_average1(self):
 
     def test_state_average_mix(self):
         mc = mcscf.CASSCF(m, 4, 4)
-        cis1 = copy.copy(mc.fcisolver)
+        cis1 = mc.fcisolver.copy()
         cis1.spin = 2
         mc = mcscf.addons.state_average_mix(mc, [cis1, mc.fcisolver], [.5, .5])
         mc.run()
diff --git a/pyscf/mcscf/test/test_ucasci.py b/pyscf/mcscf/test/test_ucasci.py
index 117a6a974f..2e67576194 100644
--- a/pyscf/mcscf/test/test_ucasci.py
+++ b/pyscf/mcscf/test/test_ucasci.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 import numpy
 from pyscf import lib
diff --git a/pyscf/mcscf/test/test_umc1step.py b/pyscf/mcscf/test/test_umc1step.py
index b4aeca34b2..ab47b23e0d 100644
--- a/pyscf/mcscf/test/test_umc1step.py
+++ b/pyscf/mcscf/test/test_umc1step.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import copy
 import unittest
 import numpy
 from pyscf import lib
@@ -96,7 +95,7 @@ def test_state_average(self):
 
     def test_state_average_mix(self):
         mc = mcscf.UCASSCF(m, 5, (4, 2))
-        cis1 = copy.copy(mc.fcisolver)
+        cis1 = mc.fcisolver.copy()
         cis1.spin = 0
         mc = mcscf.addons.state_average_mix(mc, [cis1, mc.fcisolver], [0.5, 0.5])
         mc.kernel()
diff --git a/pyscf/mcscf/umc1step.py b/pyscf/mcscf/umc1step.py
index e409f4661d..7af8dc1b9f 100644
--- a/pyscf/mcscf/umc1step.py
+++ b/pyscf/mcscf/umc1step.py
@@ -23,7 +23,6 @@
 
 import sys
 
-import copy
 from functools import reduce
 import numpy
 import pyscf.gto
@@ -310,8 +309,6 @@ def kernel(casscf, mo_coeff, tol=1e-7, conv_tol_grad=None,
         totinner += njk
 
         eris = None
-        u = copy.copy(u)
-        g_orb = copy.copy(g_orb)
         mo = casscf.rotate_mo(mo, u, log)
         eris = casscf.ao2mo(mo)
         t2m = log.timer('update eri', *t3m)
diff --git a/pyscf/mcscf/umc2step.py b/pyscf/mcscf/umc2step.py
index f490f80baa..3e946464ef 100644
--- a/pyscf/mcscf/umc2step.py
+++ b/pyscf/mcscf/umc2step.py
@@ -23,7 +23,6 @@
 
 
 import numpy
-import copy
 import pyscf.lib.logger as logger
 
 def kernel(casscf, mo_coeff, tol=1e-7, conv_tol_grad=None,
@@ -78,8 +77,6 @@ def kernel(casscf, mo_coeff, tol=1e-7, conv_tol_grad=None,
             t3m = log.timer('orbital rotation', *t3m)
 
             eris = None
-            u = copy.copy(u)
-            g_orb = copy.copy(g_orb)
             mo = casscf.rotate_mo(mo, u, log)
             eris = casscf.ao2mo(mo)
             t3m = log.timer('update eri', *t3m)
diff --git a/pyscf/mp/gmp2.py b/pyscf/mp/gmp2.py
index d1ad8dc09c..f91bcc3863 100644
--- a/pyscf/mp/gmp2.py
+++ b/pyscf/mp/gmp2.py
@@ -17,7 +17,6 @@
 E(MP2) = 1/4 <ij||ab><ab||ij>/(ei+ej-ea-eb)
 '''
 
-import copy
 import numpy
 from pyscf import lib
 from pyscf import ao2mo
@@ -203,7 +202,7 @@ def density_fit(self, auxbasis=None, with_df=None):
         if with_df is not None:
             mymp.with_df = with_df
         if mymp.with_df.auxbasis != auxbasis:
-            mymp.with_df = copy.copy(mymp.with_df)
+            mymp.with_df = mymp.with_df.copy()
             mymp.with_df.auxbasis = auxbasis
         return mymp
 
diff --git a/pyscf/mp/mp2.py b/pyscf/mp/mp2.py
index cfeec8e2ce..b6adc6ebc4 100644
--- a/pyscf/mp/mp2.py
+++ b/pyscf/mp/mp2.py
@@ -19,7 +19,6 @@
 '''
 
 
-import copy
 import numpy
 from pyscf import gto
 from pyscf import lib
@@ -636,7 +635,7 @@ def density_fit(self, auxbasis=None, with_df=None):
         if with_df is not None:
             mymp.with_df = with_df
         if mymp.with_df.auxbasis != auxbasis:
-            mymp.with_df = copy.copy(mymp.with_df)
+            mymp.with_df = mymp.with_df.copy()
             mymp.with_df.auxbasis = auxbasis
         return mymp
 
diff --git a/pyscf/mp/test/test_dfmp2_native.py b/pyscf/mp/test/test_dfmp2_native.py
index 6b26dc9f53..3399e2ada3 100644
--- a/pyscf/mp/test/test_dfmp2_native.py
+++ b/pyscf/mp/test/test_dfmp2_native.py
@@ -15,7 +15,6 @@
 
 
 import unittest
-import copy
 import numpy
 from pyscf import gto
 from pyscf import scf
@@ -60,7 +59,7 @@ class KnownValues(unittest.TestCase):
 
     def setUp(self):
         self.assertTrue(mf.converged)
-        self.mf = copy.copy(mf)
+        self.mf = mf.copy()
         self.mf.mol = mf.mol.copy()
         self.mf.mo_coeff = mf.mo_coeff.copy()
         self.mf.mo_occ = mf.mo_occ.copy()
diff --git a/pyscf/mp/test/test_dfump2_native.py b/pyscf/mp/test/test_dfump2_native.py
index 50841d4e52..dd9d79ece3 100644
--- a/pyscf/mp/test/test_dfump2_native.py
+++ b/pyscf/mp/test/test_dfump2_native.py
@@ -15,7 +15,6 @@
 
 
 import unittest
-import copy
 import numpy
 from pyscf import gto
 from pyscf import scf
@@ -57,7 +56,7 @@ class KnownValues(unittest.TestCase):
 
     def setUp(self):
         self.assertTrue(mf.converged)
-        self.mf = copy.copy(mf)
+        self.mf = mf.copy()
         self.mf.mol = mf.mol.copy()
         self.mf.mo_coeff = mf.mo_coeff.copy()
         self.mf.mo_occ = mf.mo_occ.copy()
diff --git a/pyscf/pbc/cc/test/test_eom_krccsd.py b/pyscf/pbc/cc/test/test_eom_krccsd.py
index d30b64152c..ec35fc8822 100644
--- a/pyscf/pbc/cc/test/test_eom_krccsd.py
+++ b/pyscf/pbc/cc/test/test_eom_krccsd.py
@@ -1,4 +1,3 @@
-import copy
 import pyscf.pbc.tools.make_test_cell as make_test_cell
 import numpy
 import numpy as np
@@ -317,7 +316,7 @@ def test_n3_ee(self):
     @unittest.skip('Highly sensitive to numerical noise')
     def test_t3p2_imds_complex_slow(self):
         '''Test `_slow` t3p2 implementation.'''
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         # Round to make this insensitive to small changes between PySCF versions
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
@@ -341,7 +340,7 @@ def test_t3p2_imds_complex_slow(self):
     @unittest.skip('Highly sensitive to numerical noise')
     def test_t3p2_imds_complex(self):
         '''Test t3p2 implementation.'''
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         # Round to make this insensitive to small changes between PySCF versions
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
@@ -366,7 +365,7 @@ def test_t3p2_imds_complex(self):
     def test_t3p2_imds_complex_against_so(self):
         '''Test t3[2] implementation against spin-orbital implementation.'''
         from pyscf.pbc.scf.addons import convert_to_ghf
-        kmf = copy.copy(rand_kmf2)
+        kmf = rand_kmf2.copy()
         # Round to make this insensitive to small changes between PySCF versions
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
@@ -412,7 +411,7 @@ def test_t3p2_imds_complex_against_so(self):
     def test_t3p2_imds_complex_against_so_frozen(self):
         '''Test t3[2] implementation against spin-orbital implementation with frozen orbitals.'''
         from pyscf.pbc.scf.addons import convert_to_ghf
-        kmf = copy.copy(rand_kmf2)
+        kmf = rand_kmf2.copy()
         # Round to make this insensitive to small changes between PySCF versions
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
diff --git a/pyscf/pbc/cc/test/test_krccsd.py b/pyscf/pbc/cc/test/test_krccsd.py
index 8016d57807..32d7ab3d45 100644
--- a/pyscf/pbc/cc/test/test_krccsd.py
+++ b/pyscf/pbc/cc/test/test_krccsd.py
@@ -17,7 +17,6 @@
 #          Timothy Berkelbach <tim.berkelbach@gmail.com>
 #
 
-import copy
 import unittest
 import numpy as np
 
@@ -713,7 +712,7 @@ def test_ccsd_t_high_cost(self):
     def test_rccsd_t_vs_gccsd_t(self):
         '''Test rccsd(t) vs gccsd(t) with k-points.'''
         from pyscf.pbc.scf.addons import convert_to_ghf
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
         kmf.get_veff = lambda *x: mat_veff
@@ -747,7 +746,7 @@ def test_rccsd_t_vs_gccsd_t(self):
 
     def test_rand_ccsd(self):
         '''Single (eom-)ccsd iteration with random t1/t2.'''
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
         kmf.get_veff = lambda *x: mat_veff
@@ -779,7 +778,7 @@ def test_rand_ccsd(self):
     def test_rand_ccsd_frozen0(self):
         '''Single (eom-)ccsd iteration with random t1/t2 and lowest lying orbital
         at multiple k-points frozen.'''
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
         kmf.get_veff = lambda *x: mat_veff
@@ -820,7 +819,7 @@ def test_rand_ccsd_frozen0(self):
     def test_rand_ccsd_frozen1(self):
         '''Single (eom-)ccsd iteration with random t1/t2 and single frozen occupied
         orbital.'''
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
         kmf.get_veff = lambda *x: mat_veff
@@ -863,7 +862,7 @@ def test_rand_ccsd_frozen1(self):
     def test_rand_ccsd_frozen2(self):
         '''Single (eom-)ccsd iteration with random t1/t2 and full occupied frozen
         at a single k-point.'''
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
         kmf.get_veff = lambda *x: mat_veff
@@ -906,7 +905,7 @@ def test_rand_ccsd_frozen2(self):
     def test_rand_ccsd_frozen3(self):
         '''Single (eom-)ccsd iteration with random t1/t2 and single frozen virtual
         orbital.'''
-        kmf = copy.copy(rand_kmf)
+        kmf = rand_kmf.copy()
         mat_veff = kmf.get_veff().round(4)
         mat_hcore = kmf.get_hcore().round(4)
         kmf.get_veff = lambda *x: mat_veff
diff --git a/pyscf/pbc/df/aft.py b/pyscf/pbc/df/aft.py
index 90d2dfbc74..fd6da4af9b 100644
--- a/pyscf/pbc/df/aft.py
+++ b/pyscf/pbc/df/aft.py
@@ -19,7 +19,6 @@
 '''Density expansion on plane waves'''
 
 
-import copy
 import contextlib
 import numpy as np
 from pyscf import lib
@@ -231,7 +230,7 @@ def weighted_coulG(mydf, kpt=np.zeros(3), exx=False, mesh=None, omega=None):
 def _fake_nuc(cell, with_pseudo=True):
     '''A fake cell with steep gaussians to mimic nuclear density
     '''
-    fakenuc = copy.copy(cell)
+    fakenuc = cell.copy(deep=False)
     fakenuc._atm = cell._atm.copy()
     fakenuc._atm[:,gto.PTR_COORD] = np.arange(gto.PTR_ENV_START,
                                               gto.PTR_ENV_START+cell.natm*3,3)
@@ -536,7 +535,7 @@ def range_coulomb(self, omega):
         if key in self._rsh_df:
             rsh_df = self._rsh_df[key]
         else:
-            rsh_df = self._rsh_df[key] = copy.copy(self).reset()
+            rsh_df = self._rsh_df[key] = self.copy().reset()
             if hasattr(self, '_dataname'):
                 rsh_df._dataname = f'{self._dataname}/lr/{key}'
             logger.info(self, 'Create RSH-DF object %s for omega=%s', rsh_df, omega)
@@ -686,7 +685,7 @@ def get_jk(self, dm, hermi=1, kpts=None, kpts_band=None,
     get_mo_pairs_G = get_mo_pairs = aft_ao2mo.get_mo_pairs_G
 
     def update_mf(self, mf):
-        mf = copy.copy(mf)
+        mf = mf.copy()
         mf.with_df = self
         return mf
 
diff --git a/pyscf/pbc/df/df.py b/pyscf/pbc/df/df.py
index 1c0ad54779..fad11b04f2 100644
--- a/pyscf/pbc/df/df.py
+++ b/pyscf/pbc/df/df.py
@@ -30,7 +30,6 @@
 '''
 
 import os
-import copy
 import ctypes
 import warnings
 import tempfile
@@ -447,7 +446,7 @@ def get_jk(self, dm, hermi=1, kpts=None, kpts_band=None,
     ao2mo_7d = df_ao2mo.ao2mo_7d
 
     def update_mp(self):
-        mf = copy.copy(self)
+        mf = self.copy()
         mf.with_df = self
         return mf
 
diff --git a/pyscf/pbc/df/df_jk.py b/pyscf/pbc/df/df_jk.py
index 0f60943d2e..535db5ea18 100644
--- a/pyscf/pbc/df/df_jk.py
+++ b/pyscf/pbc/df/df_jk.py
@@ -23,7 +23,6 @@
 '''
 
 
-import copy
 from functools import reduce
 import numpy
 from pyscf import lib
@@ -64,7 +63,7 @@ def density_fit(mf, auxbasis=None, mesh=None, with_df=None):
         if mesh is not None:
             with_df.mesh = mesh
 
-    mf = copy.copy(mf)
+    mf = mf.copy()
     mf.with_df = with_df
     mf._eri = None
     return mf
diff --git a/pyscf/pbc/df/fft.py b/pyscf/pbc/df/fft.py
index 0f99ef73eb..9c24059457 100644
--- a/pyscf/pbc/df/fft.py
+++ b/pyscf/pbc/df/fft.py
@@ -18,7 +18,6 @@
 
 '''Density expansion on plane waves'''
 
-import copy
 import numpy
 from pyscf import lib
 from pyscf import gto
@@ -322,7 +321,7 @@ def get_jk(self, dm, hermi=1, kpts=None, kpts_band=None,
     get_mo_pairs_G = get_mo_pairs = fft_ao2mo.get_mo_pairs_G
 
     def update_mf(self, mf):
-        mf = copy.copy(mf)
+        mf = mf.copy()
         mf.with_df = self
         return mf
 
diff --git a/pyscf/pbc/df/ft_ao.py b/pyscf/pbc/df/ft_ao.py
index c00c4c07e3..2fb27f2180 100644
--- a/pyscf/pbc/df/ft_ao.py
+++ b/pyscf/pbc/df/ft_ao.py
@@ -21,7 +21,6 @@
 '''
 
 import ctypes
-import copy
 import numpy as np
 from pyscf import lib
 from pyscf import gto
@@ -426,7 +425,7 @@ def smooth_basis_cell(self):
 
     def compact_basis_cell(self):
         '''Construct a cell with only the smooth part of the AO basis'''
-        cell_c = copy.copy(self)
+        cell_c = self.copy(deep=False)
         mask = self.bas_type != SMOOTH_BASIS
         cell_c._bas = self._bas[mask]
         cell_c.bas_map = cell_c.bas_map[mask]
diff --git a/pyscf/pbc/df/gdf_builder.py b/pyscf/pbc/df/gdf_builder.py
index 25879f64ba..5fdb64a89c 100644
--- a/pyscf/pbc/df/gdf_builder.py
+++ b/pyscf/pbc/df/gdf_builder.py
@@ -25,7 +25,6 @@
 '''
 
 import os
-import copy
 import ctypes
 import tempfile
 import numpy as np
@@ -585,7 +584,7 @@ def _int_nuc_vloc(self, fakenuc, intor='int3c2e', aosym='s2', comp=None,
         charge = -cell.atom_charges()
         if cell.dimension > 0:
             mod_cell = self.modchg_cell
-            fakenuc = copy.copy(fakenuc)
+            fakenuc = fakenuc.copy(deep=False)
             fakenuc._atm, fakenuc._bas, fakenuc._env = \
                     gto.conc_env(mod_cell._atm, mod_cell._bas, mod_cell._env,
                                  fakenuc._atm, fakenuc._bas, fakenuc._env)
@@ -767,7 +766,7 @@ def auxbar(fused_cell):
 def make_modchg_basis(auxcell, smooth_eta):
     # * chgcell defines smooth gaussian functions for each angular momentum for
     #   auxcell. The smooth functions may be used to carry the charge
-    chgcell = copy.copy(auxcell)  # smooth model density for coulomb integral to carry charge
+    chgcell = auxcell.copy(deep=False)  # smooth model density for coulomb integral to carry charge
     half_sph_norm = .5/np.sqrt(np.pi)
     chg_bas = []
     chg_env = [smooth_eta]
@@ -806,7 +805,7 @@ def fuse(Lpq, axis=0):
         return auxcell, fuse
 
     chgcell = make_modchg_basis(auxcell, eta)
-    fused_cell = copy.copy(auxcell)
+    fused_cell = auxcell.copy(deep=False)
     fused_cell._atm, fused_cell._bas, fused_cell._env = \
             gto.conc_env(auxcell._atm, auxcell._bas, auxcell._env,
                          chgcell._atm, chgcell._bas, chgcell._env)
@@ -917,7 +916,7 @@ def _guess_eta(cell, kpts=None, mesh=None):
 
 def _compensate_nuccell(cell, eta):
     '''A cell of the compensated Gaussian charges for nucleus'''
-    modchg_cell = copy.copy(cell)
+    modchg_cell = cell.copy(deep=False)
     half_sph_norm = .5/np.sqrt(np.pi)
     norm = half_sph_norm/gto.gaussian_int(2, eta)
     chg_env = [eta, norm]
diff --git a/pyscf/pbc/df/incore.py b/pyscf/pbc/df/incore.py
index ab3abd5b2f..c250590e29 100644
--- a/pyscf/pbc/df/incore.py
+++ b/pyscf/pbc/df/incore.py
@@ -17,7 +17,6 @@
 #
 
 import ctypes
-import copy
 import numpy as np
 from pyscf import lib
 from pyscf import gto
diff --git a/pyscf/pbc/df/mdf_jk.py b/pyscf/pbc/df/mdf_jk.py
index 94b089ea44..de305a736c 100644
--- a/pyscf/pbc/df/mdf_jk.py
+++ b/pyscf/pbc/df/mdf_jk.py
@@ -22,7 +22,6 @@
 J. Chem. Phys. 147, 164119 (2017)
 '''
 
-import copy
 import numpy
 from pyscf.lib import logger
 from pyscf.pbc.df import df_jk
@@ -63,7 +62,7 @@ def density_fit(mf, auxbasis=None, mesh=None, with_df=None):
         if mesh is not None:
             with_df.mesh = mesh
 
-    mf = copy.copy(mf)
+    mf = mf.copy()
     mf.with_df = with_df
     mf._eri = None
     return mf
diff --git a/pyscf/pbc/df/rsdf_builder.py b/pyscf/pbc/df/rsdf_builder.py
index a1639f14ef..1696dbfdda 100644
--- a/pyscf/pbc/df/rsdf_builder.py
+++ b/pyscf/pbc/df/rsdf_builder.py
@@ -28,7 +28,6 @@
 '''
 
 import os
-import copy
 import ctypes
 import warnings
 import tempfile
diff --git a/pyscf/pbc/df/rsdf_helper.py b/pyscf/pbc/df/rsdf_helper.py
index 5c1fdfab22..bc93221ecc 100644
--- a/pyscf/pbc/df/rsdf_helper.py
+++ b/pyscf/pbc/df/rsdf_helper.py
@@ -17,7 +17,6 @@
 #
 
 import ctypes
-import copy
 import h5py
 import numpy as np
 from scipy.special import gamma, gammaincc, comb
@@ -918,7 +917,7 @@ def intor_j2c(cell, omega, precision=None, kpts=None, hermi=1, shls_slice=None,
     fintor = getattr(mol_gto.moleintor.libcgto, intor)
     cintopt = lib.c_null_ptr()
 
-    pcell = copy.copy(cell)
+    pcell = cell.copy(deep=False)
     pcell.precision = min(cell.precision, cell.precision)
     pcell._atm, pcell._bas, pcell._env = \
             atm, bas, env = mol_gto.conc_env(cell._atm, cell._bas, cell._env,
@@ -1179,7 +1178,7 @@ def wrap_int3c_nospltbas(cell, auxcell, omega, shlpr_mask, prescreening_data,
 
 # GTO data
     intor = cell._add_suffix(intor)
-    pcell = copy.copy(cell)
+    pcell = cell.copy(deep=False)
     pcell._atm, pcell._bas, pcell._env = \
             atm, bas, env = mol_gto.conc_env(cell._atm, cell._bas, cell._env,
                                              cell._atm, cell._bas, cell._env)
diff --git a/pyscf/pbc/df/rsdf_jk.py b/pyscf/pbc/df/rsdf_jk.py
index 057b172911..08b330aeda 100644
--- a/pyscf/pbc/df/rsdf_jk.py
+++ b/pyscf/pbc/df/rsdf_jk.py
@@ -17,7 +17,6 @@
 #
 
 import numpy
-import copy
 
 def density_fit(mf, auxbasis=None, mesh=None, with_df=None):
     '''Generate density-fitting SCF object
@@ -47,7 +46,7 @@ def density_fit(mf, auxbasis=None, mesh=None, with_df=None):
         if mesh is not None:
             with_df.mesh = mesh
 
-    mf = copy.copy(mf)
+    mf = mf.copy()
     mf.with_df = with_df
     mf._eri = None
     return mf
diff --git a/pyscf/pbc/dft/multigrid.py b/pyscf/pbc/dft/multigrid.py
index 2b79c3d686..248db41f75 100644
--- a/pyscf/pbc/dft/multigrid.py
+++ b/pyscf/pbc/dft/multigrid.py
@@ -19,7 +19,6 @@
 '''Multigrid to compute DFT integrals'''
 
 import ctypes
-import copy
 import numpy
 import scipy.linalg
 
@@ -661,7 +660,7 @@ def _eval_rho_bra(cell, dms, shls_slice, hermi, xctype, kpts, grids,
     b = numpy.linalg.inv(a.T)
     ish0, ish1, jsh0, jsh1 = shls_slice
     nshells_j = jsh1 - jsh0
-    pcell = copy.copy(cell)
+    pcell = cell.copy(deep=False)
     rest_dms = []
     rest_bas = []
     i1 = 0
@@ -740,7 +739,7 @@ def _eval_rho_ket(cell, dms, shls_slice, hermi, xctype, kpts, grids,
     b = numpy.linalg.inv(a.T)
     ish0, ish1, jsh0, jsh1 = shls_slice
     nshells_i = ish1 - ish0
-    pcell = copy.copy(cell)
+    pcell = cell.copy(deep=False)
     rest_dms = []
     rest_bas = []
     j1 = 0
@@ -1544,7 +1543,7 @@ def multi_grids_tasks_for_rcut(cell, fft_mesh=None, verbose=None):
     ao_loc = cell.ao_loc_nr()
 
     def make_cell_dense_exp(shls_dense, r0, r1):
-        cell_dense = copy.copy(cell)
+        cell_dense = cell.copy(deep=False)
         cell_dense._bas = cell._bas.copy()
         cell_dense._env = cell._env.copy()
 
@@ -1575,7 +1574,7 @@ def make_cell_dense_exp(shls_dense, r0, r1):
         return cell_dense, ao_idx, ke_cutoff, rcut_atom
 
     def make_cell_sparse_exp(shls_sparse, r0, r1):
-        cell_sparse = copy.copy(cell)
+        cell_sparse = cell.copy(deep=False)
         cell_sparse._bas = cell._bas.copy()
         cell_sparse._env = cell._env.copy()
 
@@ -1663,7 +1662,7 @@ def multi_grids_tasks_for_ke_cut(cell, fft_mesh=None, verbose=None):
 
     # cell that needs dense integration grids
     def make_cell_dense_exp(shls_dense, ke0, ke1):
-        cell_dense = copy.copy(cell)
+        cell_dense = cell.copy(deep=False)
         cell_dense._bas = cell._bas.copy()
         cell_dense._env = cell._env.copy()
 
@@ -1695,7 +1694,7 @@ def make_cell_dense_exp(shls_dense, ke0, ke1):
 
     # cell that needs sparse integration grids
     def make_cell_sparse_exp(shls_sparse, ke0, ke1):
-        cell_sparse = copy.copy(cell)
+        cell_sparse = cell.copy(deep=False)
         cell_sparse._bas = cell._bas.copy()
         cell_sparse._env = cell._env.copy()
 
diff --git a/pyscf/pbc/gto/cell.py b/pyscf/pbc/gto/cell.py
index 56eec7878a..4a3deea101 100644
--- a/pyscf/pbc/gto/cell.py
+++ b/pyscf/pbc/gto/cell.py
@@ -199,7 +199,6 @@ def intor_cross(intor, cell1, cell2, comp=None, hermi=0, kpts=None, kpt=None,
 
         \langle \mu | intor | \nu \rangle, \mu \in cell1, \nu \in cell2
     '''
-    import copy
     intor, comp = moleintor._get_intor_and_comp(cell1._add_suffix(intor), comp)
 
     if kpts is None:
@@ -211,7 +210,7 @@ def intor_cross(intor, cell1, cell2, comp=None, hermi=0, kpts=None, kpt=None,
         kpts_lst = np.reshape(kpts, (-1,3))
     nkpts = len(kpts_lst)
 
-    pcell = copy.copy(cell1)
+    pcell = cell1.copy(deep=False)
     pcell.precision = min(cell1.precision, cell2.precision)
     pcell._atm, pcell._bas, pcell._env = \
             atm, bas, env = conc_env(cell1._atm, cell1._bas, cell1._env,
@@ -795,7 +794,6 @@ def _split_basis(cell, delimiter=EXP_DELIMITER):
     Split the contracted basis to small segmant.  The new basis has more
     shells.  Each shell has less primitive basis and thus is more local.
     '''
-    import copy
     _bas = []
     _env = cell._env.copy()
     contr_coeff = []
@@ -830,7 +828,7 @@ def _split_basis(cell, delimiter=EXP_DELIMITER):
                 count += 1
         contr_coeff.append(np.vstack([np.eye(degen*nc)] * count))
 
-    pcell = copy.copy(cell)
+    pcell = cell.copy(deep=False)
     pcell._bas = np.asarray(np.vstack(_bas), dtype=np.int32)
     pcell._env = _env
     return pcell, scipy.linalg.block_diag(*contr_coeff)
diff --git a/pyscf/pbc/gto/ecp.py b/pyscf/pbc/gto/ecp.py
index 73c51f8b4a..a3eff1f2f1 100644
--- a/pyscf/pbc/gto/ecp.py
+++ b/pyscf/pbc/gto/ecp.py
@@ -21,7 +21,6 @@
 '''
 
 from functools import reduce
-import copy
 import numpy
 from pyscf import lib
 from pyscf.pbc import gto
@@ -39,7 +38,7 @@ def ecp_int(cell, kpts=None):
     cell, contr_coeff = gto.cell._split_basis(cell)
     lib.logger.debug1(cell, 'nao %d -> nao %d', *(contr_coeff.shape))
 
-    ecpcell = copy.copy(cell)
+    ecpcell = cell.copy(deep=False)
     # append a fake s function to mimic the auxiliary index in pbc.incore.
     exp_ptr = cell._ecpbas[-1,PTR_EXP]
     ecpcell._bas = numpy.array([[0, 0, 1, 1, 0, exp_ptr, 0, 0]], dtype=numpy.int32)
diff --git a/pyscf/pbc/gto/pseudo/pp_int.py b/pyscf/pbc/gto/pseudo/pp_int.py
index f230365ca0..e6c62c6781 100644
--- a/pyscf/pbc/gto/pseudo/pp_int.py
+++ b/pyscf/pbc/gto/pseudo/pp_int.py
@@ -24,7 +24,6 @@
 '''
 
 import ctypes
-import copy
 import numpy
 import scipy.special
 from pyscf import lib
@@ -191,7 +190,7 @@ def fake_cell_vloc(cell, cn=0):
                 fake_bas.append([ia, 0, 1, 1, 0, ptr, ptr+1, 0])
                 ptr += 2
 
-    fakecell = copy.copy(cell)
+    fakecell = cell.copy(deep=False)
     fakecell._atm = numpy.asarray(fake_atm, dtype=numpy.int32).reshape(-1, gto.ATM_SLOTS)
     fakecell._bas = numpy.asarray(fake_bas, dtype=numpy.int32).reshape(-1, gto.BAS_SLOTS)
     fakecell._env = numpy.asarray(numpy.hstack(fake_env), dtype=numpy.double)
@@ -248,7 +247,7 @@ def fake_cell_vnl(cell):
                     hl_blocks.append(hl)
                     ptr += 2
 
-    fakecell = copy.copy(cell)
+    fakecell = cell.copy(deep=False)
     fakecell._atm = numpy.asarray(fake_atm, dtype=numpy.int32)
     fakecell._bas = numpy.asarray(fake_bas, dtype=numpy.int32)
     fakecell._env = numpy.asarray(numpy.hstack(fake_env), dtype=numpy.double)
diff --git a/pyscf/pbc/mpicc/kccsd_rhf.py b/pyscf/pbc/mpicc/kccsd_rhf.py
index 7741ae0e9f..bb1234fb14 100644
--- a/pyscf/pbc/mpicc/kccsd_rhf.py
+++ b/pyscf/pbc/mpicc/kccsd_rhf.py
@@ -17,7 +17,6 @@
 #          Timothy Berkelbach <tim.berkelbach@gmail.com>
 #
 
-import copy
 from functools import reduce
 
 import numpy
@@ -987,7 +986,7 @@ def energy_tril(cc, t1, t2, eris):
 
 def _update_procs_mf(mf):
     '''Update mean-field objects to be the same on all processors'''
-    mf1 = copy.copy(mf)
+    mf1 = mf.copy()
 
     mo_coeff  = comm.bcast(mf.mo_coeff, root=0)
     mo_energy = comm.bcast(mf.mo_energy, root=0)
diff --git a/pyscf/pbc/scf/addons.py b/pyscf/pbc/scf/addons.py
index 5ad6f30a08..bc4e5c4f8d 100644
--- a/pyscf/pbc/scf/addons.py
+++ b/pyscf/pbc/scf/addons.py
@@ -17,7 +17,6 @@
 #         Timothy Berkelbach <tim.berkelbach@gmail.com>
 #
 
-import copy
 from functools import reduce
 import numpy
 import scipy.linalg
@@ -287,7 +286,7 @@ def convert_to_uhf(mf, out=None):
 
     if out is None:
         if isinstance(mf, (scf.uhf.UHF, scf.kuhf.KUHF)):
-            return copy.copy(mf)
+            return mf.copy()
         else:
             if isinstance(mf, scf.kghf.KGHF):
                 raise NotImplementedError(
@@ -342,13 +341,13 @@ def convert_to_rhf(mf, out=None):
 
     elif nelec[0] != nelec[1] and isinstance(mf, scf.rohf.ROHF):
         if getattr(mf, '_scf', None):
-            return mol_addons._update_mf_without_soscf(mf, copy.copy(mf._scf), False)
+            return mol_addons._update_mf_without_soscf(mf, mf._scf.copy(), False)
         else:
-            return copy.copy(mf)
+            return mf.copy()
 
     else:
         if isinstance(mf, (scf.hf.RHF, scf.khf.KRHF)):
-            return copy.copy(mf)
+            return mf.copy()
         else:
             if isinstance(mf, scf.kghf.KGHF):
                 raise NotImplementedError(
@@ -400,7 +399,7 @@ def convert_to_ghf(mf, out=None):
 
     if isinstance(mf, scf.ghf.GHF):
         if out is None:
-            return copy.copy(mf)
+            return mf.copy()
         else:
             out.__dict__.update(mf.__dict__)
             return out
diff --git a/pyscf/pbc/scf/rsjk.py b/pyscf/pbc/scf/rsjk.py
index e0510d048d..d0cb9593bf 100644
--- a/pyscf/pbc/scf/rsjk.py
+++ b/pyscf/pbc/scf/rsjk.py
@@ -24,7 +24,6 @@
     Q. Sun, arXiv:2302.11307
 '''
 
-import copy
 import ctypes
 import numpy as np
 import scipy.linalg
diff --git a/pyscf/pbc/tools/pbc.py b/pyscf/pbc/tools/pbc.py
index 447f752590..681a54a5b9 100644
--- a/pyscf/pbc/tools/pbc.py
+++ b/pyscf/pbc/tools/pbc.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import warnings
-import copy
 import numpy as np
 import scipy.linalg
 from pyscf import lib
@@ -442,7 +441,7 @@ def precompute_exx(cell, kpts):
 
 def madelung(cell, kpts):
     Nk = get_monkhorst_pack_size(cell, kpts)
-    ecell = copy.copy(cell)
+    ecell = cell.copy(deep=False)
     ecell._atm = np.array([[1, cell._env.size, 0, 0, 0, 0]])
     ecell._env = np.append(cell._env, [0., 0., 0.])
     ecell.unit = 'B'
@@ -585,7 +584,7 @@ def super_cell(cell, ncopy, wrap_around=False):
         zs[(ncopy[2]+1)//2:] -= ncopy[2]
     Ts = lib.cartesian_prod((xs, ys, zs))
     Ls = np.dot(Ts, a)
-    supcell = copy.copy(cell)
+    supcell = cell.copy(deep=False)
     supcell.a = np.einsum('i,ij->ij', ncopy, a)
     mesh = np.asarray(ncopy) * np.asarray(cell.mesh)
     supcell.mesh = (mesh // 2) * 2 + 1
@@ -611,7 +610,7 @@ def cell_plus_imgs(cell, nimgs):
                              np.arange(-nimgs[1], nimgs[1]+1),
                              np.arange(-nimgs[2], nimgs[2]+1)))
     Ls = np.dot(Ts, a)
-    supcell = copy.copy(cell)
+    supcell = cell.copy(deep=False)
     supcell.a = np.einsum('i,ij->ij', nimgs, a)
     supcell.mesh = np.array([(nimgs[0]*2+1)*cell.mesh[0],
                              (nimgs[1]*2+1)*cell.mesh[1],
diff --git a/pyscf/pbc/x2c/sfx2c1e.py b/pyscf/pbc/x2c/sfx2c1e.py
index 95903161a3..cd1b1a26c2 100644
--- a/pyscf/pbc/x2c/sfx2c1e.py
+++ b/pyscf/pbc/x2c/sfx2c1e.py
@@ -23,7 +23,6 @@
 
 
 from functools import reduce
-import copy
 import numpy
 import scipy.linalg
 from pyscf import lib
diff --git a/pyscf/pbc/x2c/x2c1e.py b/pyscf/pbc/x2c/x2c1e.py
index 4e04489e13..a255b93c53 100644
--- a/pyscf/pbc/x2c/x2c1e.py
+++ b/pyscf/pbc/x2c/x2c1e.py
@@ -26,7 +26,6 @@
 
 
 from functools import reduce
-import copy
 import numpy
 import scipy.linalg
 from pyscf import lib
diff --git a/pyscf/scf/addons.py b/pyscf/scf/addons.py
index bd20c0dc7f..c3498d3360 100644
--- a/pyscf/scf/addons.py
+++ b/pyscf/scf/addons.py
@@ -18,7 +18,6 @@
 #          Junzi Liu <latrix1247@gmail.com>
 #          Susi Lehtola <susi.lehtola@gmail.com>
 
-import copy
 from functools import reduce
 import numpy
 import scipy.linalg
@@ -785,9 +784,9 @@ def convert_to_uhf(mf, out=None, remove_df=False):
         # returned SCF object, whether with_df exists in SOSCF has no effects on the
         # mean-field energy and other properties.
         if getattr(mf, '_scf', None):
-            return _update_mf_without_soscf(mf, copy.copy(mf._scf), remove_df)
+            return _update_mf_without_soscf(mf, mf._scf.copy(), remove_df)
         else:
-            return copy.copy(mf)
+            return mf.copy()
 
     else:
         known_cls = {
@@ -890,9 +889,9 @@ def convert_to_rhf(mf, out=None, remove_df=False):
     elif (isinstance(mf, scf.hf.RHF) or
           (nelec[0] != nelec[1] and isinstance(mf, scf.rohf.ROHF))):
         if getattr(mf, '_scf', None):
-            return _update_mf_without_soscf(mf, copy.copy(mf._scf), remove_df)
+            return _update_mf_without_soscf(mf, mf._scf.copy(), remove_df)
         else:
-            return copy.copy(mf)
+            return mf.copy()
 
     else:
         if nelec[0] == nelec[1]:
@@ -948,9 +947,9 @@ def convert_to_ghf(mf, out=None, remove_df=False):
 
     elif isinstance(mf, scf.ghf.GHF):
         if getattr(mf, '_scf', None):
-            return _update_mf_without_soscf(mf, copy.copy(mf._scf), remove_df)
+            return _update_mf_without_soscf(mf, mf._scf.copy(), remove_df)
         else:
-            return copy.copy(mf)
+            return mf.copy()
 
     else:
         known_cls = {
@@ -1119,7 +1118,6 @@ def fast_newton(mf, mo_coeff=None, mo_occ=None, dm0=None,
     Newton solver attributes [max_cycle_inner, max_stepsize, ah_start_tol,
     ah_conv_tol, ah_grad_trust_region, ...] can be passed through **newton_kwargs.
     '''
-    import copy
     from pyscf.lib import logger
     from pyscf import df
     from pyscf.soscf import newton_ah
@@ -1141,7 +1139,7 @@ def fast_newton(mf, mo_coeff=None, mo_occ=None, dm0=None,
         approx_grids.level = max(0, mf.grids.level-3)
         mf1.grids = approx_grids
 
-        approx_numint = copy.copy(mf._numint)
+        approx_numint = mf._numint.copy()
         mf1._numint = approx_numint
     for key in newton_kwargs:
         setattr(mf1, key, newton_kwargs[key])
@@ -1156,7 +1154,7 @@ def fast_newton(mf, mo_coeff=None, mo_occ=None, dm0=None,
         logger.note(mf, 'Generating initial guess with DIIS-SCF for newton solver')
         logger.note(mf, '========================================================')
         if dual_basis:
-            mf0 = copy.copy(mf)
+            mf0 = mf.copy()
             mf0.mol = pmol
             mf0 = mf0.density_fit(auxbasis)
         else:
diff --git a/pyscf/scf/atom_hf.py b/pyscf/scf/atom_hf.py
index 7657863b12..aac6e44316 100644
--- a/pyscf/scf/atom_hf.py
+++ b/pyscf/scf/atom_hf.py
@@ -16,7 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
-import copy
 import numpy
 from pyscf import gto
 from pyscf.lib import logger
@@ -29,7 +28,7 @@ def get_atm_nrhf(mol, atomic_configuration=elements.NRSRHF_CONFIGURATION):
     elements = set([a[0] for a in mol._atom])
     logger.info(mol, 'Spherically averaged atomic HF for %s', elements)
 
-    atm_template = copy.copy(mol)
+    atm_template = mol.copy(deep=False)
     atm_template.charge = 0
     atm_template.symmetry = False  # TODO: enable SO3 symmetry here
     atm_template.atom = atm_template._atom = []
diff --git a/pyscf/scf/atom_ks.py b/pyscf/scf/atom_ks.py
index 77ad838b84..d259a0968d 100644
--- a/pyscf/scf/atom_ks.py
+++ b/pyscf/scf/atom_ks.py
@@ -16,7 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #         Susi Lehtola <susi.lehtola@gmail.com>
 
-import copy
 import numpy
 from pyscf import gto
 from pyscf.lib import logger
@@ -29,7 +28,7 @@ def get_atm_nrks(mol, atomic_configuration=elements.NRSRHFS_CONFIGURATION, xc='s
     elements = set([a[0] for a in mol._atom])
     logger.info(mol, 'Spherically averaged atomic KS for %s', elements)
 
-    atm_template = copy.copy(mol)
+    atm_template = mol.copy(deep=False)
     atm_template.charge = 0
     atm_template.symmetry = False  # TODO: enable SO3 symmetry here
     atm_template.atom = atm_template._atom = []
diff --git a/pyscf/scf/test/test_addons.py b/pyscf/scf/test/test_addons.py
index 3c5eaed2a1..53633f5420 100644
--- a/pyscf/scf/test/test_addons.py
+++ b/pyscf/scf/test/test_addons.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import unittest
-import copy
 import numpy
 import scipy.linalg
 
@@ -307,7 +306,7 @@ def test_convert_to_scf(self):
         self.assertTrue(isinstance(addons.convert_to_uhf(sym_mf_u), scf.uhf_symm.UHF))
         self.assertTrue(isinstance(addons.convert_to_ghf(sym_mf_u), scf.ghf_symm.GHF))
 
-        mf1 = copy.copy(mf)
+        mf1 = mf.copy()
         self.assertTrue(isinstance(mf1.convert_from_(mf), scf.rhf.RHF))
         self.assertTrue(isinstance(mf1.convert_from_(mf_u), scf.rhf.RHF))
         self.assertFalse(isinstance(mf1.convert_from_(mf_u), scf.rohf.ROHF))
@@ -316,7 +315,7 @@ def test_convert_to_scf(self):
         self.assertFalse(isinstance(mf1.convert_from_(sym_mf_u), scf.rohf.ROHF))
         self.assertFalse(isinstance(mf1.convert_from_(sym_mf), scf.hf_symm.RHF))
         self.assertFalse(isinstance(mf1.convert_from_(sym_mf_u), scf.hf_symm.RHF))
-        mf1 = copy.copy(mf_u)
+        mf1 = mf_u.copy()
         self.assertTrue(isinstance(mf1.convert_from_(mf), scf.uhf.UHF))
         self.assertTrue(isinstance(mf1.convert_from_(mf_u), scf.uhf.UHF))
         self.assertTrue(isinstance(mf1.convert_from_(sym_mf), scf.uhf.UHF))
diff --git a/pyscf/scf/test/test_dhf.py b/pyscf/scf/test/test_dhf.py
index 5651eef297..0fbf3dc22f 100644
--- a/pyscf/scf/test/test_dhf.py
+++ b/pyscf/scf/test/test_dhf.py
@@ -16,7 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
-import copy
 import numpy
 import unittest
 from pyscf import gto
@@ -125,7 +124,7 @@ def test_get_veff(self):
         v = mf.get_veff(mol, dm)
         self.assertAlmostEqual(lib.fp(v), (-21.613084684028077-28.50754366262467j), 8)
 
-        mf1 = copy.copy(mf)
+        mf1 = mf.copy()
         mf1.direct_scf = False
         v1 = mf1.get_veff(mol, dm)
         self.assertAlmostEqual(abs(v-v1).max(), 0, 9)
diff --git a/pyscf/scf/test/test_ghf.py b/pyscf/scf/test/test_ghf.py
index cfc74e7370..c029870ea7 100644
--- a/pyscf/scf/test/test_ghf.py
+++ b/pyscf/scf/test/test_ghf.py
@@ -16,7 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
-import copy
 import unittest
 import tempfile
 import numpy
@@ -243,7 +242,7 @@ def test_canonicalize(self):
         self.assertAlmostEqual(abs(e - mfsym.mo_energy).max(), 0, 6)
 
     def test_get_occ(self):
-        mf1 = copy.copy(mfsym)
+        mf1 = mfsym.copy()
         mf1.irrep_nelec = {}
         mf1.irrep_nelec['B1'] = 1
         occ = mf1.get_occ(mf.mo_energy, mf.mo_coeff+0j)
@@ -346,7 +345,7 @@ def test_get_irrep_nelec(self):
         self.assertEqual(irrep_nelec['B1'], 2)
         self.assertEqual(irrep_nelec['B2'], 2)
 
-        mf1 = copy.copy(mfsym)
+        mf1 = mfsym.copy()
         mf1.irrep_nelec = irrep_nelec
         mf1.irrep_nelec['A1'] = 2
         mf1.irrep_nelec['A2'] = 2
diff --git a/pyscf/scf/test/test_jk.py b/pyscf/scf/test/test_jk.py
index 31cfa3eae1..c379a8425d 100644
--- a/pyscf/scf/test/test_jk.py
+++ b/pyscf/scf/test/test_jk.py
@@ -16,7 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
-import copy
 import numpy
 import unittest
 from pyscf import lib
@@ -102,7 +101,7 @@ def test_shls_slice1(self):
         self.assertAlmostEqual(abs(k1part - k1ref).max(), 0, 12)
 
     def test_mols(self):
-        pmol = copy.copy(mol)
+        pmol = mol.copy(deep=False)
         mols = (mol, pmol, pmol, mol)
         dm = mf.make_rdm1()
         vj0 = jk.get_jk(mols, dm, 'ijkl,lk->ij')
diff --git a/pyscf/scf/test/test_rhf.py b/pyscf/scf/test/test_rhf.py
index edeca8a0b2..048192187d 100644
--- a/pyscf/scf/test/test_rhf.py
+++ b/pyscf/scf/test/test_rhf.py
@@ -16,7 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
-import copy
 import numpy
 import unittest
 import tempfile
@@ -328,11 +327,11 @@ def test_analyze(self):
         popandchg, dip = mf1.analyze(with_meta_lowdin=False)
         self.assertAlmostEqual(numpy.linalg.norm(popandchg[0]), 3.2031790129016922, 6)
 
-        mf1 = copy.copy(n2mf)
+        mf1 = n2mf.copy()
         (pop, chg), dip = n2mf.analyze()
         self.assertAlmostEqual(numpy.linalg.norm(pop), 4.5467414321488357, 6)
         self.assertAlmostEqual(numpy.linalg.norm(dip), 0, 9)
-        mf1 = copy.copy(n2mf)
+        mf1 = n2mf.copy()
         mf1.mo_coeff = numpy.array(n2mf.mo_coeff)
         popandchg, dip = mf1.analyze(with_meta_lowdin=False)
         self.assertAlmostEqual(numpy.linalg.norm(popandchg[0]), 3.8893148995392353, 6)
@@ -747,7 +746,7 @@ def test_uniq_var(self):
         self.assertAlmostEqual(abs(f1 + f1.T).max(), 0, 12)
 
     def test_check_convergence(self):
-        mf1 = copy.copy(n2mf)
+        mf1 = n2mf.copy()
         mf1.diis = False
         count = [0]
         def check_convergence(envs):
@@ -781,7 +780,7 @@ def test_get_irrep_nelec(self):
         self.assertEqual(irrep_nelec['A1g'], 6)
         self.assertEqual(irrep_nelec['E1ux'], 2)
         self.assertEqual(irrep_nelec['E1uy'], 2)
-        n2_rhf = copy.copy(n2mf)
+        n2_rhf = n2mf.copy()
         n2_rhf.irrep_nelec = irrep_nelec
         n2_rhf.irrep_nelec['A2g'] = 0
         n2_rhf.irrep_nelec['E2gx'] = 2
diff --git a/pyscf/scf/test/test_uhf.py b/pyscf/scf/test/test_uhf.py
index 13d14c51be..fc7607d3d5 100644
--- a/pyscf/scf/test/test_uhf.py
+++ b/pyscf/scf/test/test_uhf.py
@@ -16,7 +16,6 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
-import copy
 import numpy
 import scipy.linalg
 import unittest
@@ -318,7 +317,7 @@ def test_analyze(self):
         (pop, chg), dip = mf.analyze(with_meta_lowdin=False)
         self.assertAlmostEqual(numpy.linalg.norm(pop[0]+pop[1]), 3.2031790129016922, 6)
 
-        mf1 = copy.copy(n2mf)
+        mf1 = n2mf.copy()
         (pop, chg), dip = n2mf.analyze()
         self.assertAlmostEqual(numpy.linalg.norm(pop[0]+pop[1]), 4.5467414321488357, 6)
         self.assertAlmostEqual(numpy.linalg.norm(dip), 0, 9)
diff --git a/pyscf/sgx/sgx.py b/pyscf/sgx/sgx.py
index 7f8153da3e..6eada197bd 100644
--- a/pyscf/sgx/sgx.py
+++ b/pyscf/sgx/sgx.py
@@ -20,7 +20,6 @@
 Pseudo-spectral methods (COSX, PS, SN-K)
 '''
 
-import copy
 import numpy
 from pyscf import lib
 from pyscf import gto
@@ -83,7 +82,7 @@ def sgx_fit(mf, auxbasis=None, with_df=None, pjs=False):
             mf.with_df = with_df
         elif mf.with_df.auxbasis != auxbasis:
             #logger.warn(mf, 'DF might have been initialized twice.')
-            mf = copy.copy(mf)
+            mf = mf.copy()
             mf.with_df = with_df
         return mf
 
@@ -348,7 +347,7 @@ def get_jk(self, dm, hermi=1, vhfopt=None, with_j=True, with_k=True,
             if key in self._rsh_df:
                 rsh_df = self._rsh_df[key]
             else:
-                rsh_df = copy.copy(self)
+                rsh_df = self.copy()
                 rsh_df._rsh_df = None  # to avoid circular reference
                 # Not all attributes need to be reset. Resetting _vjopt
                 # because it is used by get_j method of regular DF object.
diff --git a/pyscf/tdscf/test/test_tddks.py b/pyscf/tdscf/test/test_tddks.py
index 366768a66a..7b8b093330 100644
--- a/pyscf/tdscf/test/test_tddks.py
+++ b/pyscf/tdscf/test/test_tddks.py
@@ -19,7 +19,6 @@
 import unittest
 import tempfile
 import numpy
-import copy
 from pyscf import lib, gto, scf, dft
 from pyscf import tdscf
 try:
diff --git a/pyscf/tdscf/test/test_tdgks.py b/pyscf/tdscf/test/test_tdgks.py
index 0847813d7b..c0177149b7 100644
--- a/pyscf/tdscf/test/test_tdgks.py
+++ b/pyscf/tdscf/test/test_tdgks.py
@@ -19,7 +19,6 @@
 import unittest
 import tempfile
 import numpy
-import copy
 from pyscf import lib, gto, scf, dft
 from pyscf import tdscf
 try:
diff --git a/pyscf/tdscf/test/test_tdrks.py b/pyscf/tdscf/test/test_tdrks.py
index fef61007cc..febc276e9f 100644
--- a/pyscf/tdscf/test/test_tdrks.py
+++ b/pyscf/tdscf/test/test_tdrks.py
@@ -18,7 +18,6 @@
 
 import unittest
 import numpy
-import copy
 from pyscf import lib, gto, scf, dft
 from pyscf.tdscf import rhf, rks
 from pyscf import tdscf
@@ -295,7 +294,7 @@ def test_nto(self):
         self.assertAlmostEqual(w[0], 0.99997335352278072, 7)
         self.assertAlmostEqual(lib.fp(w), 0.99998775067586554, 7)
 
-        pmol = copy.copy(mol)
+        pmol = mol.copy(deep=False)
         pmol.symmetry = True
         pmol.build(0, 0)
         mf = scf.RHF(pmol).run()
diff --git a/pyscf/tdscf/test/test_tduks.py b/pyscf/tdscf/test/test_tduks.py
index b4c4a42581..9679fc5ea5 100644
--- a/pyscf/tdscf/test/test_tduks.py
+++ b/pyscf/tdscf/test/test_tduks.py
@@ -18,7 +18,6 @@
 
 import unittest
 import numpy
-import copy
 from pyscf import lib, gto, scf, dft
 from pyscf import tdscf
 
@@ -351,7 +350,7 @@ def test_nto(self):
         self.assertAlmostEqual(lib.fp(w[0]), 0.00027305600430816, 9)
         self.assertAlmostEqual(lib.fp(w[1]), 0.99964370569529093, 9)
 
-        pmol = copy.copy(mol)
+        pmol = mol.copy(deep=False)
         pmol.symmetry = True
         pmol.build(0, 0)
         mf = scf.UHF(pmol).run()
diff --git a/pyscf/tools/test/test_cubegen.py b/pyscf/tools/test/test_cubegen.py
index 3272571621..6a383a7e1f 100644
--- a/pyscf/tools/test/test_cubegen.py
+++ b/pyscf/tools/test/test_cubegen.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import unittest
-import copy
 import tempfile
 from pyscf import lib, gto, scf
 from pyscf.tools import cubegen
diff --git a/pyscf/tools/test/test_molden.py b/pyscf/tools/test/test_molden.py
index f7664b1b50..ce74c2f0bb 100644
--- a/pyscf/tools/test/test_molden.py
+++ b/pyscf/tools/test/test_molden.py
@@ -14,7 +14,6 @@
 # limitations under the License.
 
 import unittest
-import copy
 import tempfile
 from pyscf import lib, gto, scf
 from pyscf.tools import molden
diff --git a/pyscf/x2c/test/test_tdscf.py b/pyscf/x2c/test/test_tdscf.py
index b98d9a2c59..ad14d4c6f4 100644
--- a/pyscf/x2c/test/test_tdscf.py
+++ b/pyscf/x2c/test/test_tdscf.py
@@ -19,7 +19,6 @@
 import unittest
 import tempfile
 import numpy
-import copy
 from pyscf import lib, gto, scf
 from pyscf.x2c import x2c, dft, tdscf
 try:
diff --git a/pyscf/x2c/x2c.py b/pyscf/x2c/x2c.py
index b0b8f7b8f5..d7f6e38a63 100644
--- a/pyscf/x2c/x2c.py
+++ b/pyscf/x2c/x2c.py
@@ -18,7 +18,6 @@
 '''
 
 from functools import reduce
-import copy
 import numpy
 import scipy.linalg
 from pyscf import lib
@@ -58,7 +57,7 @@ def get_xmol(self, mol=None):
             mol = self.mol
 
         if self.basis is not None:
-            xmol = copy.copy(mol)
+            xmol = mol.copy(deep=False)
             xmol.build(False, False, basis=self.basis)
             return xmol, None
         elif self.xuncontract: