diff --git a/pyscf/df/addons.py b/pyscf/df/addons.py
index 38ab630c37..620d900306 100644
--- a/pyscf/df/addons.py
+++ b/pyscf/df/addons.py
@@ -16,6 +16,7 @@
 # Author: Qiming Sun <osirpt.sun@gmail.com>
 #
 
+
 import sys
 import numpy
 from pyscf.lib import logger
@@ -119,8 +120,8 @@ def aug_etb_for_dfbasis(mol, dfbasis=DFBASIS, beta=ETB_BETA,
             emax_by_l = [emax[liljsum==ll].max() for ll in range(l_max1*2-1)]
             emin_by_l = [emin[liljsum==ll].min() for ll in range(l_max1*2-1)]
             # Tune emin and emax
-            emin_by_l = numpy.array(emin_by_l) * 2  # *2 for alpha+alpha on same center
-            emax_by_l = numpy.array(emax_by_l) * 2  #/ (numpy.arange(l_max1*2-1)*.5+1)
+            emin_by_l = numpy.array(emin_by_l) * 2
+            emax_by_l = numpy.array(emax_by_l) * 2
 
             ns = numpy.log((emax_by_l+emin_by_l)/emin_by_l) / numpy.log(beta)
             etb = []
@@ -187,6 +188,7 @@ def make_auxbasis(mol, mp2fit=False):
                 logger.debug(mol, '  ETB auxbasis for %s  %s', k, auxbasis[k])
     return auxbasis
 
+# TODO: add auxbasis keyword etb and auto
 def make_auxmol(mol, auxbasis=None):
     '''Generate a fake Mole object which uses the density fitting auxbasis as
     the basis sets.  If auxbasis is not specified, the optimized auxiliary fitting
@@ -197,7 +199,6 @@ def make_auxmol(mol, auxbasis=None):
     even-tempered Gaussian basis set will be generated.
 
     See also the paper JCTC, 13, 554 about generating auxiliary fitting basis.
-    JCP, 127, 074102
 
     Kwargs:
         auxbasis : str, list, tuple
diff --git a/pyscf/df/auto_basis.py b/pyscf/df/auto_basis.py
new file mode 100644
index 0000000000..7b2f1b5532
--- /dev/null
+++ b/pyscf/df/auto_basis.py
@@ -0,0 +1,25 @@
+'''
+JCTC, 13, 554
+'''
+
+f = [20, 4.0, 4.0, 3.5, 2.5, 2.0, 2.0]
+beta_big = [1.8, 2.0, 2.2, 2.2, 2.2, 2.3, 3.0, 3.0]
+beta_small = 1.8
+
+def auto_basis(mol, l_inc=1):
+    l_val = highest_l_in_occupied_shell
+    l_max = max(mol._bas[:,0])
+    l_max_aux = max(l_val*2, l_max+l_inc)
+
+    a_min_aux = a_min[:,None] + a_min
+    a_max_aux = a_min[:,None] + a_min
+    r_exp = _gaussian_int(a, l)
+    a_max_eff = 2 * k(l)**2 / (np.pi * r_exp)
+    a_max_eff_aux = a_max_eff[:,None] + a_max_eff
+    a_max_aux = ([max(f[l] * a_max_eff_aux, a_max_aux) for l in range(2*l_val+1)] +
+                 a_max_aux[2*l_val+1:])
+
+    beta = list(beta_big)
+    beta[:2*l_val] = beta_small
+
+    ns = np.log(a_max_aux/a_min_aux) / beta