diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 3cc91ba1..dc7367d6 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -21,7 +21,7 @@ jobs:
       - name: test
         run: ./.github/workflows/run_test.sh
       - name: Upload coverage to codecov
-        uses: codecov/codecov-action@v2
+        uses: codecov/codecov-action@v3
         with:
           token: ${{secrets.CODECOV_TOKEN}}
           files: ./pyscfad/coverage.xml
diff --git a/.pylintrc b/.pylintrc
index a9996da1..4193f161 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -437,7 +437,7 @@ valid-metaclass-classmethod-first-arg=mcs
 
 # Exceptions that will emit a warning when being caught. Defaults to
 # "Exception"
-overgeneral-exceptions=StandardError,
-                       Exception,
-                       BaseException
+overgeneral-exceptions=builtins.StandardError,
+                       builtins.Exception,
+                       builtins.BaseException
 
diff --git a/README.md b/README.md
index 2b0bf11b..1ccb8414 100644
--- a/README.md
+++ b/README.md
@@ -30,13 +30,19 @@ pip install 'pyscf-properties @ git+https://github.com/fishjojo/properties.git@a
 pip install pyscfad
 ```
 
-* To install the development version of pyscfad, use the following command instead:
+---
+* To install the development version, use the following command instead:
 ```
 pip install git+https://github.com/fishjojo/pyscfad.git
 ```
----
 
-* To Install from source and install dependencies manually
+* The dependencies can be installed via a predefined conda environment
+```
+conda env create -f environment.yml
+conda activate pyscfad_env
+```
+
+* Alternatively, the dependencies can be installed from source
 ```
 pip install numpy scipy h5py
 pip install jax jaxlib jaxopt
@@ -51,7 +57,7 @@ export PYTHONPATH=$HOME/pyscf:$PYTHONPATH
 ```
 
 ---
-* Running pyscfad inside a docker container:
+* One can also run PySCFAD inside a docker container:
 ```
 docker pull fishjojo/pyscfad:latest
 docker run -rm -t -i fishjojo/pyscfad:latest /bin/bash
@@ -60,7 +66,9 @@ docker run -rm -t -i fishjojo/pyscfad:latest /bin/bash
 Running examples
 ----------------
 
-* Add the following lines to the PySCF configure file ($HOME/.pyscf\_conf.py)
+* In order to perform AD calculations, 
+the following lines need to be added to 
+the PySCF configure file($HOME/.pyscf\_conf.py)
 ```
 pyscfad = True
 pyscf_numpy_backend = 'jax'
diff --git a/codecov.yml b/codecov.yml
new file mode 100644
index 00000000..322b2974
--- /dev/null
+++ b/codecov.yml
@@ -0,0 +1,11 @@
+comment: false
+coverage:
+  status:
+    project:
+      default:
+        threshold: 1%
+        target: 80%
+    patch:
+      default:
+        threshold: 1%
+        informational: true
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 00000000..c8743ca6
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,14 @@
+name: pyscfad_env
+channels:
+  - defaults
+dependencies:
+  - python=3.9
+  - pip=23.0
+  - pip:
+    - numpy>=1.17
+    - scipy
+    - jax>=0.1.65,<0.3.14
+    - jaxlib>=0.1.65,<0.3.14
+    - jaxopt>=0.2
+    - -e git+https://github.com/fishjojo/pyscf.git@ad#egg=pyscf
+    - -e git+https://github.com/fishjojo/properties.git@ad#egg=pyscf-properties
diff --git a/examples/scf/03-rohf.py b/examples/scf/03-rohf.py
new file mode 100644
index 00000000..917998a8
--- /dev/null
+++ b/examples/scf/03-rohf.py
@@ -0,0 +1,34 @@
+import pyscf
+from pyscfad import gto, scf
+
+"""
+Analytic nuclear gradient for ROHF computed by auto-differentiation
+
+Reference results from PySCF:
+converged SCF energy = -75.578154312784
+--------------- ROHF gradients ---------------
+         x                y                z
+0 O     0.0000000000    -0.0000000000    -0.0023904882
+1 H     0.0000000000    -0.0432752607     0.0011952441
+2 H    -0.0000000000     0.0432752607     0.0011952441
+----------------------------------------------
+"""
+
+mol = gto.Mole()
+mol.atom = '''
+O        0.000000    0.000000    0.117790
+H        0.000000    0.755453   -0.471161
+H        0.000000   -0.755453   -0.471161'''
+mol.basis = '631g'
+mol.charge = 1
+mol.spin = 1  # = 2S = spin_up - spin_down
+mol.verbose = 5
+mol.build()
+
+mf = scf.ROHF(mol)
+mf.kernel()
+
+jac = mf.energy_grad()
+print(f'Nuclaer gradient:\n{jac.coords}')
+print(f'Gradient wrt basis exponents:\n{jac.exp}')
+print(f'Gradient wrt basis contraction coefficients:\n{jac.ctr_coeff}')
diff --git a/pyscfad/gto/mole.py b/pyscfad/gto/mole.py
index 1d7c5b29..1055094b 100644
--- a/pyscfad/gto/mole.py
+++ b/pyscfad/gto/mole.py
@@ -64,11 +64,11 @@ def __init__(self, **kwargs):
         self.exp = None
         self.ctr_coeff = None
         self.r0 = None
-        gto.Mole.__init__(self, **kwargs)
+        super().__init__(**kwargs)
 
     def atom_coords(self, unit='Bohr'):
         if self.coords is None:
-            return gto.Mole.atom_coords(self, unit)
+            return super().atom_coords(unit)
         else:
             if unit[:3].upper() == 'ANG':
                 return self.coords * param.BOHR
@@ -81,7 +81,7 @@ def build(self, *args, **kwargs):
         trace_ctr_coeff = kwargs.pop('trace_ctr_coeff', True)
         trace_r0 = kwargs.pop('trace_r0', False)
 
-        gto.Mole.build(self, *args, **kwargs)
+        super().build(*args, **kwargs)
 
         if trace_coords:
             self.coords = np.asarray(self.atom_coords())
@@ -99,8 +99,8 @@ def intor(self, intor, comp=None, hermi=0, aosym='s1', out=None,
               shls_slice=None, grids=None):
         if (self.coords is None and self.exp is None
                 and self.ctr_coeff is None and self.r0 is None):
-            return gto.Mole.intor(self, intor, comp=comp, hermi=hermi,
-                                  aosym=aosym, out=out, shls_slice=shls_slice)
+            return super().intor(intor, comp=comp, hermi=hermi,
+                                 aosym=aosym, out=out, shls_slice=shls_slice)
         else:
             return moleintor.getints(self, intor, shls_slice,
                                      comp, hermi, aosym, out=None)
diff --git a/pyscfad/scf/__init__.py b/pyscfad/scf/__init__.py
index f034ab6d..6e7d5b11 100644
--- a/pyscfad/scf/__init__.py
+++ b/pyscfad/scf/__init__.py
@@ -1,8 +1,12 @@
 from pyscfad.scf import hf
 from pyscfad.scf import uhf
+from pyscfad.scf import rohf
 
 def RHF(mol, **kwargs):
     return hf.RHF(mol, **kwargs)
 
 def UHF(mol, **kwargs):
     return uhf.UHF(mol, **kwargs)
+
+def ROHF(mol, **kwargs):
+    return rohf.ROHF(mol, **kwargs)
diff --git a/pyscfad/scf/hf.py b/pyscfad/scf/hf.py
index e73b996d..070dda11 100644
--- a/pyscfad/scf/hf.py
+++ b/pyscfad/scf/hf.py
@@ -1,4 +1,4 @@
-from functools import partial
+from functools import partial, reduce
 import numpy
 from jax import jacrev, jacfwd
 from jaxopt import linear_solve
@@ -138,6 +138,8 @@ def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,
         fock = mf.get_fock(h1e, s1e, vhf, dm)
         mo_energy, mo_coeff = mf.eig(fock, s1e)
         mo_occ = mf.get_occ(stop_grad(mo_energy), stop_grad(mo_coeff))
+        # hack for ROHF
+        mo_energy = getattr(mo_energy, 'mo_energy', mo_energy)
         return scf_conv, e_tot, mo_energy, mo_coeff, mo_occ
 
     if isinstance(mf.diis, lib.diis.DIIS):
@@ -204,6 +206,9 @@ def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,
         if dump_chk:
             mf.dump_chk(locals())
 
+    # hack for ROHF
+    mo_energy = getattr(mo_energy, 'mo_energy', mo_energy)
+
     log.timer('scf_cycle', *cput0)
     del log
     # A post-processing hook before return
@@ -235,6 +240,11 @@ def _dot_eri_dm_nosymm(eri, dm, with_j, with_k):
         vk = vk.reshape(dm.shape)
     return vj, vk
 
+def level_shift(s, d, f, factor):
+    dm_vir = s - reduce(np.dot, (s, d, s))
+    return f + dm_vir * factor
+
+
 @util.pytree_node(Traced_Attributes, num_args=1)
 class SCF(pyscf_hf.SCF):
     '''
diff --git a/pyscfad/scf/rohf.py b/pyscfad/scf/rohf.py
new file mode 100644
index 00000000..f84c2896
--- /dev/null
+++ b/pyscfad/scf/rohf.py
@@ -0,0 +1,228 @@
+from functools import reduce
+import numpy
+from pyscf import numpy as np
+from pyscf.lib import stop_grad
+from pyscf.lib import logger
+from pyscf.scf import chkfile
+from pyscf.scf import rohf as pyscf_rohf
+from pyscfad import util
+from pyscfad.scf import hf
+
+@util.pytree_node(['fock', 'focka', 'fockb'], num_args=1)
+class _FockMatrix():
+    def __init__(self, fock, focka=None, fockb=None):
+        self.fock = fock
+        self.focka = focka
+        self.fockb = fockb
+
+    def __repr__(self):
+        return self.fock.__repr__()
+
+@util.pytree_node(['mo_energy',], num_args=1)
+class _OrbitalEnergy():
+    def __init__(self, mo_energy, mo_ea=None, mo_eb=None):
+        self.mo_energy = mo_energy
+        self.mo_ea = mo_ea
+        self.mo_eb = mo_eb
+
+    def __repr__(self):
+        return self.mo_energy.__repr__()
+
+def get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None,
+             diis_start_cycle=None, level_shift_factor=None, damp_factor=None):
+    if h1e is None: h1e = mf.get_hcore()
+    if s1e is None: s1e = mf.get_ovlp()
+    if vhf is None: vhf = mf.get_veff(mf.mol, dm)
+    if dm is None: dm = mf.make_rdm1()
+    if getattr(dm, 'ndim', None) == 2:
+        dm = np.array((dm*.5, dm*.5))
+# To Get orbital energy in get_occ, we saved alpha and beta fock, because
+# Roothaan effective Fock cannot provide correct orbital energy with `eig`
+# TODO, check other treatment  J. Chem. Phys. 133, 141102
+    focka = h1e + vhf[0]
+    fockb = h1e + vhf[1]
+    f = get_roothaan_fock((focka,fockb), dm, s1e)
+    if cycle < 0 and diis is None:  # Not inside the SCF iteration
+        return _FockMatrix(f, focka, fockb)
+
+    if diis_start_cycle is None:
+        diis_start_cycle = mf.diis_start_cycle
+    if level_shift_factor is None:
+        level_shift_factor = mf.level_shift
+    if damp_factor is None:
+        damp_factor = mf.damp
+
+    dm_tot = dm[0] + dm[1]
+    if 0 <= cycle < diis_start_cycle-1 and abs(damp_factor) > 1e-4:
+        raise NotImplementedError('ROHF Fock-damping')
+    if diis and cycle >= diis_start_cycle:
+        f = diis.update(s1e, dm_tot, f, mf, h1e, vhf)
+    if abs(level_shift_factor) > 1e-4:
+        f = hf.level_shift(s1e, dm_tot*.5, f, level_shift_factor)
+    return _FockMatrix(f, focka, fockb)
+
+def get_roothaan_fock(focka_fockb, dma_dmb, s):
+    nao = s.shape[0]
+    focka, fockb = focka_fockb
+    dma, dmb = dma_dmb
+    fc = (focka + fockb) * .5
+# Projector for core, open-shell, and virtual
+    pc = np.dot(dmb, s)
+    po = np.dot(dma-dmb, s)
+    pv = np.eye(nao) - np.dot(dma, s)
+    fock  = reduce(np.dot, (pc.conj().T, fc, pc)) * .5
+    fock += reduce(np.dot, (po.conj().T, fc, po)) * .5
+    fock += reduce(np.dot, (pv.conj().T, fc, pv)) * .5
+    fock += reduce(np.dot, (po.conj().T, fockb, pc))
+    fock += reduce(np.dot, (po.conj().T, focka, pv))
+    fock += reduce(np.dot, (pv.conj().T, fc, pc))
+    fock = fock + fock.conj().T
+    return fock
+
+def get_grad(mo_coeff, mo_occ, fock):
+    occidxa = mo_occ > 0
+    occidxb = mo_occ == 2
+    viridxa = ~occidxa
+    viridxb = ~occidxb
+    uniq_var_a = viridxa.reshape(-1,1) & occidxa
+    uniq_var_b = viridxb.reshape(-1,1) & occidxb
+
+    if getattr(fock, 'focka', None) is not None:
+        focka = fock.focka
+        fockb = fock.fockb
+    elif isinstance(fock, (tuple, list)) or getattr(fock, 'ndim', None) == 3:
+        focka, fockb = fock
+    else:
+        focka = fockb = fock
+    focka = reduce(numpy.dot, (mo_coeff.conj().T, focka, mo_coeff))
+    fockb = reduce(numpy.dot, (mo_coeff.conj().T, fockb, mo_coeff))
+
+    g = numpy.zeros_like(focka)
+    g[uniq_var_a]  = focka[uniq_var_a]
+    g[uniq_var_b] += fockb[uniq_var_b]
+    return g[uniq_var_a | uniq_var_b].ravel()
+
+def make_rdm1(mo_coeff, mo_occ, **kwargs):
+    if getattr(mo_occ, 'ndim', None) == 1:
+        mo_occa = mo_occ > 0
+        mo_occb = mo_occ == 2
+    else:
+        mo_occa, mo_occb = mo_occ
+    dm_a = np.dot(mo_coeff*mo_occa, mo_coeff.conj().T)
+    dm_b = np.dot(mo_coeff*mo_occb, mo_coeff.conj().T)
+    return np.array((dm_a, dm_b))
+
+def get_occ(mf, mo_energy=None, mo_coeff=None):
+    from pyscf.scf.rohf import _fill_rohf_occ
+    if mo_energy is None: mo_energy = mf.mo_energy
+    if getattr(mo_energy, 'mo_ea', None) is not None:
+        mo_ea = numpy.asarray(mo_energy.mo_ea)
+        mo_eb = numpy.asarray(mo_energy.mo_eb)
+        mo_eab = numpy.asarray(mo_energy.mo_energy)
+    else:
+        mo_ea = mo_eb = mo_eab = numpy.asarray(mo_energy)
+    nmo = mo_ea.size
+    mo_occ = numpy.zeros(nmo)
+    if getattr(mf, 'nelec', None) is None:
+        nelec = mf.mol.nelec
+    else:
+        nelec = mf.nelec
+    if nelec[0] > nelec[1]:
+        nocc, ncore = nelec
+    else:
+        ncore, nocc = nelec
+    nopen = nocc - ncore
+    mo_occ = _fill_rohf_occ(mo_eab, mo_ea, mo_eb, ncore, nopen)
+
+    if mf.verbose >= logger.INFO and nocc < nmo and ncore > 0:
+        ehomo = max(mo_eab[mo_occ> 0])
+        elumo = min(mo_eab[mo_occ==0])
+        if ehomo+1e-3 > elumo:
+            logger.warn(mf, 'HOMO %.15g >= LUMO %.15g', ehomo, elumo)
+        else:
+            logger.info(mf, '  HOMO = %.15g  LUMO = %.15g', ehomo, elumo)
+        if nopen > 0 and mf.verbose >= logger.DEBUG:
+            core_idx = mo_occ == 2
+            open_idx = mo_occ == 1
+            vir_idx = mo_occ == 0
+            logger.debug(mf, '                  Roothaan           | alpha              | beta')
+            logger.debug(mf, '  Highest 2-occ = %18.15g | %18.15g | %18.15g',
+                         max(mo_eab[core_idx]),
+                         max(mo_ea[core_idx]), max(mo_eb[core_idx]))
+            logger.debug(mf, '  Lowest 0-occ =  %18.15g | %18.15g | %18.15g',
+                         min(mo_eab[vir_idx]),
+                         min(mo_ea[vir_idx]), min(mo_eb[vir_idx]))
+            for i in numpy.where(open_idx)[0]:
+                logger.debug(mf, '  1-occ =         %18.15g | %18.15g | %18.15g',
+                             mo_eab[i], mo_ea[i], mo_eb[i])
+
+        if mf.verbose >= logger.DEBUG:
+            numpy.set_printoptions(threshold=nmo)
+            logger.debug(mf, '  Roothaan mo_energy =\n%s', mo_eab)
+            logger.debug1(mf, '  alpha mo_energy =\n%s', mo_ea)
+            logger.debug1(mf, '  beta  mo_energy =\n%s', mo_eb)
+            numpy.set_printoptions(threshold=1000)
+    return mo_occ
+
+@util.pytree_node(hf.Traced_Attributes, num_args=1)
+class ROHF(hf.SCF, pyscf_rohf.ROHF):
+    def __init__(self, mol, **kwargs):
+        pyscf_rohf.ROHF.__init__(self, mol)
+        self.__dict__.update(kwargs)
+
+    def eig(self, fock, s, x0=None):
+        focka = getattr(fock, 'focka', None)
+        fockb = getattr(fock, 'fockb', None)
+        fockab = getattr(fock, 'fock', fock)
+        e, c = self._eigh(fockab, s, x0)
+        if focka is not None:
+            c_copy = numpy.asarray(stop_grad(c))
+            focka_copy = numpy.asarray(stop_grad(focka))
+            fockb_copy = numpy.asarray(stop_grad(fockb))
+            mo_ea = numpy.einsum('pi,pi->i', c_copy.conj(), focka_copy.dot(c_copy)).real
+            mo_eb = numpy.einsum('pi,pi->i', c_copy.conj(), fockb_copy.dot(c_copy)).real
+            e = _OrbitalEnergy(e, mo_ea, mo_eb)
+        return e, c
+
+    def get_grad(self, mo_coeff, mo_occ, fock=None):
+        if fock is None:
+            dm1 = self.make_rdm1(mo_coeff, mo_occ)
+            fock = self.get_hcore(self.mol) + self.get_veff(self.mol, dm1)
+        return get_grad(stop_grad(mo_coeff), stop_grad(mo_occ), stop_grad(fock))
+
+    def make_rdm1(self, mo_coeff=None, mo_occ=None, **kwargs):
+        if mo_coeff is None: mo_coeff = self.mo_coeff
+        if mo_occ is None: mo_occ = self.mo_occ
+        if self.mol.spin < 0:
+            # Flip occupancies of alpha and beta orbitals
+            mo_occ = (mo_occ == 2), (mo_occ > 0)
+        return make_rdm1(mo_coeff, mo_occ, **kwargs)
+
+    def get_veff(self, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
+        if mol is None: mol = self.mol
+        if dm is None: dm = self.make_rdm1()
+        if getattr(dm, 'ndim', None) == 2:
+            dm = np.array((dm*.5, dm*.5))
+
+        if self._eri is not None or not self.direct_scf:
+            vj, vk = self.get_jk(mol, dm, hermi)
+            vhf = vj[0] + vj[1] - vk
+        else:
+            ddm = dm - np.asarray(dm_last)
+            vj, vk = self.get_jk(mol, ddm, hermi)
+            vhf = vj[0] + vj[1] - vk
+            vhf += np.asarray(vhf_last)
+        return vhf
+
+    def dump_chk(self, envs):
+        if self.chkfile:
+            mo_energy = getattr(envs['mo_energy'], 'mo_energy',
+                                envs['mo_energy'])
+            chkfile.dump_scf(self.mol, self.chkfile,
+                             envs['e_tot'], mo_energy,
+                             envs['mo_coeff'], envs['mo_occ'],
+                             overwrite_mol=False)
+        return self
+
+    get_occ = get_occ
+    get_fock = get_fock
diff --git a/pyscfad/scf/test/test_rohf.py b/pyscfad/scf/test/test_rohf.py
new file mode 100644
index 00000000..5129f6b7
--- /dev/null
+++ b/pyscfad/scf/test/test_rohf.py
@@ -0,0 +1,11 @@
+from pyscfad import scf
+
+def test_nuc_grad(get_h2o_plus):
+    mol = get_h2o_plus
+    mf = scf.ROHF(mol)
+    g1 = mf.energy_grad().coords
+    mf.kernel()
+    g2 = mf.energy_grad().coords
+    g0 = mf.nuc_grad_method().kernel()
+    assert abs(g1-g0).max() < 1e-6
+    assert abs(g2-g0).max() < 1e-6
diff --git a/pyscfad/version.py b/pyscfad/version.py
index 3dc1f76b..b3f47562 100644
--- a/pyscfad/version.py
+++ b/pyscfad/version.py
@@ -1 +1 @@
-__version__ = "0.1.0"
+__version__ = "0.1.2"