diff --git a/README.md b/README.md
index 69da2288..172669f3 100644
--- a/README.md
+++ b/README.md
@@ -51,6 +51,7 @@ Features
 - Nonlocal functional correction (vv10) for SCF and gradient;
 - ECP is supported and calculated on CPU;
 - PCM solvent models, analytical gradients, and semi-analytical Hessian matrix;
+- SMD solvent models and solvation free energy
 
 Limitations
 --------
diff --git a/examples/16-smd.py b/examples/16-smd.py
new file mode 100644
index 00000000..28b864f4
--- /dev/null
+++ b/examples/16-smd.py
@@ -0,0 +1,38 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import numpy
+import pyscf
+from pyscf import gto, df
+from gpu4pyscf import scf, dft
+from gpu4pyscf.solvent import smd
+
+atom = '''
+O       0.0000000000    -0.0000000000     0.1174000000
+H      -0.7570000000    -0.0000000000    -0.4696000000
+H       0.7570000000     0.0000000000    -0.4696000000
+'''
+mol = pyscf.M(atom=atom, basis='def2-tzvpp', verbose=4)
+
+mf = dft.rks.RKS(mol, xc='HYB_GGA_XC_B3LYP')#.density_fit()
+mf = mf.SMD()
+mf.verbose = 4
+mf.grids.atom_grid = (99,590)
+mf.small_rho_cutoff = 1e-10
+mf.with_solvent.lebedev_order = 29 # 302 Lebedev grids
+mf.with_solvent.method = 'SMD'
+mf.with_solvent.solvent = 'water'
+e_tot = mf.kernel()
+print(e_tot)
diff --git a/gpu4pyscf/solvent/__init__.py b/gpu4pyscf/solvent/__init__.py
index 157e7129..5227c986 100644
--- a/gpu4pyscf/solvent/__init__.py
+++ b/gpu4pyscf/solvent/__init__.py
@@ -13,7 +13,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from gpu4pyscf.solvent import pcm
+from gpu4pyscf.solvent import pcm, smd
 
 def PCM(method_or_mol, solvent_obj=None, dm=None):
     '''Initialize PCM model.
@@ -35,3 +35,24 @@ def PCM(method_or_mol, solvent_obj=None, dm=None):
         return pcm.pcm_for_scf(method_or_mol, solvent_obj, dm)
     else:
         raise NotImplementedError('PCM model only support SCF')
+
+def SMD(method_or_mol, solvent_obj=None, dm=None):
+    '''Initialize SMD model.
+
+    Examples:
+
+    >>> mf = PCM(scf.RHF(mol))
+    >>> mf.kernel()
+    >>> sol = PCM(mol)
+    >>> mc = PCM(CASCI(mf, 6, 6), sol)
+    >>> mc.kernel()
+    '''
+    from pyscf import gto
+    from pyscf import scf
+
+    if isinstance(method_or_mol, gto.mole.Mole):
+        return smd.SMD(method_or_mol)
+    elif isinstance(method_or_mol, scf.hf.SCF):
+        return pcm.pcm_for_scf(method_or_mol, solvent_obj, dm)
+    else:
+        raise NotImplementedError('SMD model only support SCF')
\ No newline at end of file
diff --git a/gpu4pyscf/solvent/_attach_solvent.py b/gpu4pyscf/solvent/_attach_solvent.py
index 23e08af2..1df93bab 100644
--- a/gpu4pyscf/solvent/_attach_solvent.py
+++ b/gpu4pyscf/solvent/_attach_solvent.py
@@ -87,10 +87,21 @@ def energy_elec(self, dm=None, h1e=None, vhf=None):
                 dm = self.make_rdm1()
             if getattr(vhf, 'e_solvent', None) is None:
                 vhf = self.get_veff(self.mol, dm)
+
             e_tot, e_coul = oldMF.energy_elec(self, dm, h1e, vhf)
             e_tot += vhf.e_solvent
             self.scf_summary['e_solvent'] = vhf.e_solvent.real
-            logger.debug(self, 'Solvent Energy = %.15g', vhf.e_solvent)
+
+            if self.with_solvent.method.upper() == 'SMD':
+                if self.with_solvent.e_cds is None:
+                    e_cds = self.with_solvent.get_cds()
+                    self.with_solvent.e_cds = e_cds
+                else:
+                    e_cds = self.with_solvent.e_cds
+                e_tot += e_cds
+                self.scf_summary['e_cds'] = e_cds
+                logger.info(self, f'CDS correction = {e_cds:.15f}')
+            logger.info(self, 'Solvent Energy = %.15g', vhf.e_solvent)
             return e_tot, e_coul
 
         def nuc_grad_method(self):
diff --git a/gpu4pyscf/solvent/grad/pcm.py b/gpu4pyscf/solvent/grad/pcm.py
index 8298f434..be469508 100644
--- a/gpu4pyscf/solvent/grad/pcm.py
+++ b/gpu4pyscf/solvent/grad/pcm.py
@@ -71,9 +71,9 @@ def get_dF_dA(surface):
         riJ = cupy.linalg.norm(ri_rJ, axis=-1)
         diJ = (riJ - R_in_J) / R_sw_J
         diJ[:,ia] = 1.0
-        diJ[diJ < 1e-12] = 0.0
+        diJ[diJ < 1e-8] = 0.0
         ri_rJ[:,ia,:] = 0.0
-        ri_rJ[diJ < 1e-12] = 0.0
+        ri_rJ[diJ < 1e-8] = 0.0
 
         fiJ = switch_h(diJ)
         dfiJ = grad_switch_h(diJ) / (fiJ * riJ * R_sw_J)
diff --git a/gpu4pyscf/solvent/pcm.py b/gpu4pyscf/solvent/pcm.py
index 0604b046..428e5a08 100644
--- a/gpu4pyscf/solvent/pcm.py
+++ b/gpu4pyscf/solvent/pcm.py
@@ -40,7 +40,7 @@ def pcm_for_scf(mf, solvent_obj=None, dm=None):
 
 # Inject PCM to SCF, TODO: add it to other methods later
 from gpu4pyscf import scf
-scf.hf.RHF.PCM    = scf.hf.RHF.PCM    = pcm_for_scf
+scf.hf.RHF.PCM = pcm_for_scf
 
 # TABLE II,  J. Chem. Phys. 122, 194110 (2005)
 XI = {
@@ -75,8 +75,8 @@ def pcm_for_scf(mf, solvent_obj=None, dm=None):
     5810: 4.90792902522,
 }
 
-Bondi = radii.VDW
-Bondi[1] = 1.1/radii.BOHR      # modified version
+modified_Bondi = radii.VDW.copy()
+modified_Bondi[1] = 1.1/radii.BOHR      # modified version
 #radii_table = bondi * 1.2
 PI = numpy.pi
 
@@ -91,7 +91,7 @@ def switch_h(x):
     y[x>1] = 1.0
     return y
 
-def gen_surface(mol, ng=302, vdw_scale=1.2):
+def gen_surface(mol, ng=302, rad=modified_Bondi, vdw_scale=1.2, r_probe=0.0):
     '''J. Phys. Chem. A 1999, 103, 11060-11079'''
     unit_sphere = numpy.empty((ng,4))
     libdft.MakeAngularGrid(unit_sphere.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(ng))
@@ -100,7 +100,7 @@ def gen_surface(mol, ng=302, vdw_scale=1.2):
     atom_coords = cupy.asarray(mol.atom_coords(unit='B'))
     charges = mol.atom_charges()
     N_J = ng * cupy.ones(mol.natm)
-    R_J = cupy.asarray([vdw_scale*Bondi[chg] for chg in charges])
+    R_J = cupy.asarray([rad[chg] for chg in charges])
     R_sw_J = R_J * (14.0 / N_J)**0.5
     alpha_J = 1.0/2.0 + R_J/R_sw_J - ((R_J/R_sw_J)**2 - 1.0/28)**0.5
     R_in_J = R_J - alpha_J * R_sw_J
@@ -117,20 +117,20 @@ def gen_surface(mol, ng=302, vdw_scale=1.2):
     for ia in range(mol.natm):
         symb = mol.atom_symbol(ia)
         chg = gto.charge(symb)
-        r_vdw = vdw_scale*Bondi[chg]
+        r_vdw = rad[chg]
 
         atom_grid = r_vdw * unit_sphere[:,:3] + atom_coords[ia,:]
         #riJ = scipy.spatial.distance.cdist(atom_grid[:,:3], atom_coords)
         riJ = cupy.sum((atom_grid[:,None,:] - atom_coords[None,:,:])**2, axis=2)**0.5
         diJ = (riJ - R_in_J) / R_sw_J
         diJ[:,ia] = 1.0
-        diJ[diJ<1e-12] = 0.0
+        diJ[diJ<1e-8] = 0.0
 
         fiJ = switch_h(diJ)
 
         w = unit_sphere[:,3] * 4.0 * PI
         swf = cupy.prod(fiJ, axis=1)
-        idx = w*swf > 1e-16
+        idx = w*swf > 1e-12
 
         p0, p1 = p1, p1+sum(idx).get()
         gslice_by_atom.append([p0,p1])
@@ -208,10 +208,15 @@ def get_D_S(surface, with_S=True, with_D=False):
     return D, S
 
 class PCM(ddcosmo.DDCOSMO):
+    _keys = {
+        'method', 'vdw_scale', 'surface'
+    }
     def __init__(self, mol):
         ddcosmo.DDCOSMO.__init__(self, mol)
         self.method = 'C-PCM'
         self.vdw_scale = 1.2 # default value in qchem
+        self.r_probe = 0.0
+        self.radii_table = None
         self.surface = {}
         self._intermediates = {}
 
@@ -231,13 +236,14 @@ def dump_flags(self, verbose=None):
         return self
 
     def build(self, ng=None):
-        vdw_scale = self.vdw_scale
-        self.radii_table = vdw_scale * Bondi
+        if self.radii_table is None:
+            vdw_scale = self.vdw_scale
+            self.radii_table = vdw_scale * modified_Bondi + self.r_probe
         mol = self.mol
         if ng is None:
             ng = gen_grid.LEBEDEV_ORDER[self.lebedev_order]
 
-        self.surface = gen_surface(mol, ng=ng, vdw_scale=vdw_scale)
+        self.surface = gen_surface(mol, rad=self.radii_table, ng=ng)
         self._intermediates = {}
         F, A = get_F_A(self.surface)
         D, S = get_D_S(self.surface, with_S=True, with_D=True)
@@ -251,7 +257,7 @@ def build(self, ng=None):
             f_epsilon = (epsilon - 1.0)/(epsilon + 1.0/2.0)
             K = S
             R = -f_epsilon * cupy.eye(K.shape[0])
-        elif self.method.upper() in ['IEF-PCM', 'IEFPCM']:
+        elif self.method.upper() in ['IEF-PCM', 'IEFPCM', 'SMD']:
             f_epsilon = (epsilon - 1.0)/(epsilon + 1.0)
             DA = D*A
             DAS = cupy.dot(DA, S)
diff --git a/gpu4pyscf/solvent/smd.py b/gpu4pyscf/solvent/smd.py
new file mode 100644
index 00000000..5fb5a8e0
--- /dev/null
+++ b/gpu4pyscf/solvent/smd.py
@@ -0,0 +1,577 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+'''
+SMD solvent model
+'''
+
+import numpy as np
+import cupy
+from pyscf import lib
+from pyscf.data import radii
+from pyscf.dft import gen_grid
+from gpu4pyscf.solvent import pcm, _attach_solvent
+from gpu4pyscf.lib import logger
+
+@lib.with_doc(_attach_solvent._for_scf.__doc__)
+def smd_for_scf(mf, solvent_obj=None, dm=None):
+    if solvent_obj is None:
+        solvent_obj = SMD(mf.mol)
+    return _attach_solvent._for_scf(mf, solvent_obj, dm)
+
+# Inject PCM to SCF, TODO: add it to other methods later
+from gpu4pyscf import scf
+scf.hf.RHF.SMD = smd_for_scf
+
+hartree2kcal = 627.5
+# see https://pubs.acs.org/doi/epdf/10.1021/jp810292n
+
+sigma_water = {
+    ('H'): 48.69,
+    ('C'): 129.74,
+    ('H','C'): -60.77,
+    ('C','C'): -72.95,
+    ('O','C'): 68.69,
+    ('N','C'): -48.22,
+    ('N','C3'): 84.10,
+    ('O','N'): 121.98,
+    ('F'): 38.18,
+    ('Cl'): 9.82,
+    ('Br'): -8.72,
+    ('S'): -9.10,
+    ('O','P'): 68.85}
+
+sigma_n = {
+    ('C'): 58.10,
+    ('H','C'): -36.37,
+    ('C','C'): -62.05,
+    ('O'): -17.56,
+    ('H','O'): -19.39,
+    ('O','C'): -15.70,
+    ('N'): 32.62,
+    ('C','N'): -99.76,
+    ('Cl'): -24.31,
+    ('Br'): -35.42,
+    ('S'): -33.17,
+    ('Si'): -18.04
+}
+
+sigma_alpha = {
+    ('C'): 48.10,
+    ('O'): 193.06,
+    ('O','C'): 95.99,
+    ('C','N'): 152.20,
+    ('N','C'): -41.00
+}
+
+sigma_beta = {
+    ('C'): 32.87,
+    ('O'): -43.79,
+    ('O','O'):-128.16,
+    ('O','N'):79.13
+}
+
+# Molecular surface tension (cal/mol*AA^-2)
+sigma_gamma = 0.35
+sigma_phi2 = -4.19
+sigma_psi2 = -6.68
+sigma_beta2 = 0.0
+gamma0 = 1.0
+
+# rzz and delta r_zz in AA
+r_zz = {
+    ('H','C'): [1.55, 0.3],
+    ('H','O'): [1.55, 0.3],
+    ('C','H'): [1.55, 0.3],
+    ('C','C'): [1.84, 0.3],
+    ('C','N'): [1.84, 0.3],
+    ('C','O'): [1.84, 0.3],
+    ('C','F'): [1.84, 0.3],
+    ('C','P'): [2.2, 0.3],
+    ('C','S'): [2.2, 0.3],
+    ('C','Cl'): [2.1, 0.3],
+    ('C','Br'): [2.3, 0.3],
+    ('C','I'): [2.6, 0.3],
+    ('N','C'): [1.84, 0.3],
+    ('N','C3'): [1.225, 0.065],
+    ('O','C'): [1.33, 0.1],
+    ('O','N'): [1.5, 0.3],
+    ('O','O'): [1.8, 0.3],
+    ('O','P'): [2.1, 0.3]
+}
+
+# database from https://comp.chem.umn.edu/solvation/mnsddb.pdf
+# solvent name: [n, n25, alpha, beta, gamma, epsilon, phi, psi)
+solvent_db = {
+    '1,1,1-trichloroethane':[1.4379, 1.4313, 0.0, 0.09, 36.24, 7.0826, 0.0, 0.60],
+    '1,1,2-trichloroethane':[1.4717, 1.4689, 0.13, 0.13, 48.97, 7.1937, 0.0, 0.60],
+    '1,2,4-trimethylbenzene':[1.5048, 1.5024, 0.0, 0.19, 42.03, 2.3653, 0.667, 0.0],
+    '1,2-dibromoethane':[1.5387, 1.5364, 0.10, 0.17, 56.93, 4.9313, 0.0, 0.5],
+    '1,2-dichloroethane':[1.4448, 1.4425, 0.10, 0.11, 45.86, 10.125, 0.0, 0.5],
+    '1,2-ethanediol':[1.4318, 1.4306, 0.58, 0.78, 69.07, 40.245, 0.0, 0.5],
+    '1,4-dioxane':[1.4224, 1.4204, 0.00, 0.64, 47.14, 2.2099, 0.0, 0.0],
+    '1-bromo-2-methylpropane':[1.4348, 1.4349, 0.00, 0.12, 34.69, 7.7792, 0.0, 0.2],
+    '1-bromooctane':[1.4524, 1.4500, 0.0, 0.12, 41.28, 5.0244, 0.0, 0.111],
+    '1-bromopentane':[1.4447, 1.4420, 0.00, 0.12, 38.7, 6.269, 0.0, 0.167],
+    '1-bromopropane':[1.4343, 1.4315, 0.0, 0.12, 36.36, 8.0496, 0.0, 0.250],
+    '1-butanol':[1.3993, 1.3971, 0.37, 0.48, 35.88, 17.332, 0.0, 0.0],
+    '1-chlorohexane':[1.4199, -1, 0.0, 0.10, 37.03, 5.9491, 0.0, 0.143],
+    '1-chloropentane':[1.4127, 1.4104, 0.0, 0.1, 35.12, 6.5022, 0.0, 0.167],
+    '1-chloropropane':[1.3879, 1.3851, 0.0, 0.1, 30.66, 8.3548, 0.0, 0.25],
+    '1-decanol':[1.4372, 1.4353, 0.37, 0.48, 41.04, 7.5305, 0.0, 0.0],
+    '1-fluorooctane':[1.3935, 1.3927, 0.0, 0.10, 33.92, 3.89, 0.0, 0.111],
+    '1-heptanol':[1.4249, 1.4224, 0.37, 0.48, 38.5, 11.321, 0.0, 0.0],
+    '1-hexanol':[1.4178, 1.4162, 0.37, 0.48, 37.15, 12.51, 0.0, 0.0],
+    '1-hexene':[1.3837, 1.385, 0.00, 0.07, 25.76, 2.0717, 0.0, 0.0],
+    '1-hexyne':[1.3989, 1.3957, 0.12, 0.10, 28.79, 2.615, 0.0, 0.0],
+    '1-iodobutane':[1.5001, 1.4958, 0.00, 0.15, 40.65, 6.173, 0.0, 0.0],
+    '1-iodohexadecane':[1.4806, -1, 0.00, 0.15, 46.48, 3.5338, 0.0, 0.0],
+    '1-iodopentane':[1.4959, -1, 0.00, 0.15, 41.56, 5.6973, 0.0, 0.0],
+    '1-iodopropane':[1.5058, 1.5027, 0.00, 0.15, 41.45, 6.9626, 0.0, 0.0],
+    '1-nitropropane':[1.4018, 1.3996, 0.00, 0.31, 43.32, 23.73, 0.0, 0.0],
+    '1-nonanol':[1.4333, 1.4319, 0.37, 0.48, 40.14, 8.5991, 0.0, 0.0],
+    '1-octanol':[1.4295, 1.4279, 0.37, 0.48, 39.01, 9.8629, 0.0, 0.0],
+    '1-pentanol':[1.4101, 1.4080, 0.37, 0.48, 36.5, 15.13, 0.0, 0.0],
+    '1-pentene':[1.3715, 1.3684, 0.00, 0.07, 22.24, 1.9905, 0.0, 0.0],
+    '1-propanol':[1.3850, 1.3837, 0.37, 0.48, 33.57, 20.524, 0.0, 0.0],
+    '2,2,2-trifluoroethanol':[1.2907, -1, 0.57, 0.25, 42.02, 26.726, 0.0, 0.5],
+    '2,2,4-trimethylpentane':[1.3915, 1.3889, 0.00, 0.00, 26.38, 1.9358, 0.0, 0.0],
+    '2,4-dimethylpentane':[1.3815, 1.3788, 0.0, 0.00, 25.42, 1.8939, 0.0, 0.0],
+    '2,4-dimethylpyridine':[1.5010, 1.4985, 0.0, 0.63, 46.86, 9.4176, 0.625, 0.0],
+    '2,6-dimethylpyridine':[1.4953, 1.4952, 0.0, 0.63, 44.64, 7.1735, 0.625, 0.0],
+    '2-bromopropane':[1.4251, 1.4219, 0.00, 0.14, 33.46, 9.3610, 0.0, 0.25],
+    '2-butanol':[1.3978, 1.3949, 0.33, 0.56, 32.44, 15.944, 0.0, 0.0],
+    '2-chlorobutane':[1.3971, 1.3941, 0.00, 0.12, 31.1, 8.3930, 0.0, 0.2],
+    '2-heptanone':[1.4088, 1.4073, 0.0, 0.51, 37.6, 11.658, 0.0, 0.0],
+    '2-hexanone':[1.4007, 1.3987, 0.0, 0.51, 36.63, 14.136, 0.0, 0.0],
+    '2-methoxyethanol':[1.4024, 1.4003, 0.30, 0.84, 44.39, 17.2, 0.0, 0.0],
+    '2-methyl-1-propanol':[1.3955, 1.3938, 0.37, 0.48, 32.38, 16.777, 0.0, 0.0],
+    '2-methyl-2-propanol':[1.3878, 1.3852, 0.31, 0.60, 28.73, 12.47, 0.0, 0.0],
+    '2-methylpentane':[1.3715, 1.3687, 0.0, 0.00, 24.3, 1.89, 0.0, 0.0],
+    '2-methylpyridine':[1.4957, 1.4984, 0.0, 0.58, 47.5, 9.9533, 0.714, 0.0],
+    '2-nitropropane':[1.3944, 1.3923, 0.00, 0.33, 42.16, 25.654, 0.00, 0.0],
+    '2-octanone':[1.4151, 1.4133, 0.00, 0.51, 37.29, 9.4678, 0.0, 0.0],
+    '2-pentanone':[1.3895, 1.3885, 0.00, 0.51, 33.46, 15.200, 0.0, 0.0],
+    '2-propanol':[1.3776, 1.3752, 0.33, 0.56, 30.13, 19.264, 0.0, 0.0],
+    '2-propen-1-ol':[1.4135, 0.00, 0.38, 0.48, 36.39, 19.011, 0.0, 0.0],
+    'E-2-pentene':[1.3793, 1.3761, 0.00, 0.07, 23.62, 2.051, 0.0, 0.0],
+    '3-methylpyridine':[1.5040, 1.5043, 0.0, 0.54, 49.61, 11.645, 0.714, 0.0],
+    '3-pentanone':[1.3924, 1.3905, 0.00, 0.51, 35.61, 16.78, 0.0, 0.0],
+    '4-heptanone':[1.4069, 1.4045, 0.00, 0.51, 35.98, 12.257, 0.0, 0.0],
+    '4-methyl-2-pentanone':[1.3962, 1.394, 0.0, 0.51, 33.83, 12.887, 0.0, 0.0],
+    '4-methylpyridine':[1.5037, 1.503, 0.0, 0.54, 50.17, 11.957, 0.714, 0.0],
+    '5-nonanone':[1.4195, 0.00, 0.0, 0.51, 37.83, 10.6, 0.0, 0.0],
+    'acetic acid':[1.3720, 1.3698, 0.61, 0.44, 39.01, 6.2528, 0.0, 0.0],
+    'acetone':[1.3588, 1.3559, 0.04, 0.49, 33.77, 20.493, 0.0, 0.0],
+    'acetonitrile':[1.3442, 1.3416, 0.07, 0.32, 41.25, 35.688, 0.0, 0.0],
+    'acetophenone':[1.5372, 1.5321, 0.00, 0.48, 56.19, 17.44, 0.667, 0.0],
+    'aniline':[1.5863, 1.5834, 0.26, 0.41, 60.62, 6.8882, 0.857, 0.0],
+    'anisole':[1.5174, 1.5143, 0.0, 0.29, 50.52, 4.2247, 0.75, 0.0],
+    'benzaldehyde':[1.5463, 1.5433, 0.0, 0.39, 54.69, 18.220, 0.857, 0.0],
+    'benzene':[1.5011, 1.4972, 0.0, 0.14, 40.62, 2.2706, 1.0, 0.0],
+    'benzonitrile':[1.5289, 1.5257, 0.0, 0.33, 55.83, 25.592, 0.75, 0.0],
+    'benzylalcohol':[1.5396, 1.5384, 0.33, 0.56, 52.96, 12.457, 0.75, 0.0],
+    'bromobenzene':[1.5597, 1.5576, 0.0, 0.09, 50.72, 5.3954, 0.857, 0.143],
+    'bromoethane':[1.4239, 1.4187, 0.0, 0.12, 34.0, 9.01, 0.0, 0.333],
+    'bromoform':[1.6005, 1.5956, 0.15, 0.06, 64.58, 4.2488, 0.0, 0.75],
+    'butanal':[1.3843, 1.3766, 0.0, 0.45, 35.06, 13.45, 0.0, 0.0],
+    'butanoic acid':[1.3980, 1.3958, 0.60, 0.45, 37.49, 2.9931, 0.0, 0.0],
+    'butanone':[1.3788, 1.3764, 0.00, 0.51, 34.5, 18.246, 0.0, 0.0],
+    'butanonitrile':[1.3842, 1.382, 0.0, 0.36, 38.75, 24.291, 0.0, 0.0],
+    'butylethanoate':[1.3941, 1.3923, 0.0, 0.45, 35.81, 4.9941, 0.0, 0.0],
+    'butylamine':[1.4031, 1.3987, 0.16, 0.61, 33.74, 4.6178, 0.0, 0.0],
+    'n-butylbenzene':[1.4898, 1.4874, 0.0, 0.15, 41.33, 2.36, 0.6, 0.0],
+    'sec-butylbenzene':[1.4895, 1.4878, 0.0, 0.16, 40.35, 2.3446, 0.60, 0.0],
+    'tert-butylbenzene':[1.4927, 1.4902, 0.0, 0.16, 39.78, 2.3447, 0.6, 0.0],
+    'carbon disulfide':[1.6319, 1.6241, 0.0, 0.07, 45.45, 2.6105, 0.0, 0.0],
+    'carbon tetrachloride':[1.4601, 1.4574, 0.00, 0.00, 38.04, 2.2280, 0.0, 0.8],
+    'chlorobenzene':[1.5241, 1.5221, 0.0, 0.07, 47.48, 5.6968, 0.857, 0.143],
+    'chloroform':[1.4459, 1.4431, 0.15, 0.02, 38.39, 4.7113, 0.0, 0.75],
+    'a-chlorotoluene':[1.5391, 0.0, 0.0, 0.33, 53.04, 6.7175, 0.75, 0.125],
+    'o-chlorotoluene':[1.5268, 1.5233, 0.00, 0.07, 47.43, 4.6331, 0.75, 0.125],
+    'm-cresol':[1.5438, 1.5394, 0.57, 0.34, 51.37, 12.44, 0.75, 0.0],
+    'o-cresol':[1.5361, 1.5399, 0.52, 0.30, 53.11, 6.76, 0.75, 0.0],
+    'cyclohexane':[1.4266, 1.4235, 0.00, 0.00, 35.48, 2.0165, 0.0, 0.0],
+    'cyclohexanone':[1.4507, 1.4507, 0.00, 0.56, 49.76, 15.619, 0.00, 0.0],
+    'cyclopentane':[1.4065, 1.4036, 0.00, 0.00, 31.49, 1.9608, 0.0, 0.0],
+    'cyclopentanol':[1.4530, -1, 0.32, 0.56, 46.8, 16.989, 0.0, 0.0],
+    'cyclopentanone':[1.4366, 1.4347, 0.00, 0.52, 47.21, 13.58, 0.0, 0.0],
+    'decalin (cis/trans mixture)':[1.4753, 1.472, 0.00, 0.00, 43.82, 2.196, 0.0, 0.0],
+    'cis-decalin':[1.4810, 1.4788, 0.00, 0.00, 45.45, 2.2139, 0.0, 0.0],
+    'n-decane':[1.4102, 1.4094, 0.00, 0.00, 33.64, 1.9846, 0.0, 0.0],
+    'dibromomethane':[1.5420, 1.5389, 0.10, 0.10, 56.21, 7.2273, 0.0, 0.667],
+    'butylether':[1.3992, 1.3968, 0.00, 0.45, 35.98, 3.0473, 0.0, 0.0],
+    'o-dichlorobenzene':[1.5515, 1.5491, 0.00, 0.04, 52.72, 9.9949, 0.75, 0.25],
+    'E-1,2-dichloroethene':[1.4454, 1.4435, 0.09, 0.05, 37.13, 2.14, 0.0, 0.5],
+    'Z-1,2-dichloroethene':[1.4490, 1.4461, 0.11, 0.05, 39.8, 9.2, 0.0, 0.5],
+    'dichloromethane':[1.4242, 1.4212, 0.10, 0.05, 39.15, 8.93, 0.0, 0.667],
+    'diethylether':[1.3526, 1.3496, 0.00, 0.41, 23.96, 4.2400, 0.0, 0.0],
+    'diethylsulfide':[1.4430, 1.4401, 0.00, 0.32, 35.36, 5.723, 0.0, 0.0],
+    'diethylamine':[1.3864, 1.3825, 0.08, 0.69, 28.57, 3.5766, 0.0, 0.0],
+    'diiodomethane':[1.7425, 1.738, 0.05, 0.23, 95.25, 5.32, 0.0, 0.0],
+    'diisopropyl ether':[1.3679, 1.3653, 0.00, 0.41, 24.86, 3.38, 0.0, 0.0],
+    'cis-1,2-dimethylcyclohexane':[1.4360, 1.4336, 0.00, 0.00, 36.28, 2.06, 0.0, 0.0],
+    'dimethyldisulfide':[1.5289, 1.522, 0.00, 0.28, 48.06, 9.6, 0.0, 0.0],
+    'N,N-dimethylacetamide':[1.4380, 1.4358, 0.00, 0.78, 47.62, 37.781, 0.0, 0.0],
+    'N,N-dimethylformamide':[1.4305, 1.4280, 0.00, 0.74, 49.56, 37.219, 0.0, 0.0],
+    'dimethylsulfoxide':[1.4783, 1.4783, 0.00, 0.88, 61.78, 46.826, 0.0, 0.0],
+    'diphenylether':[1.5787, -1, 0.00, 0.20, 38.5, 3.73, 0.923, 0.0],
+    'dipropylamine':[1.4050, 1.4018, 0.08, 0.69, 32.11, 2.9112, 0.0, 0.0],
+    'n-dodecane':[1.4216, 1.4151, 0.00, 0.00, 35.85, 2.0060, 0.0, 0.0],
+    'ethanethiol':[1.4310, 1.4278, 0.00, 0.24, 33.22, 6.667, 0.0, 0.0],
+    'ethanol':[1.3611, 1.3593, 0.37, 0.48, 31.62, 24.852, 0.0, 0.0],
+    'ethylethanoate':[1.3723, 1.3704, 0.00, 0.45, 33.67, 5.9867, 0.0, 0.0],
+    'ethylmethanoate':[1.3599, 1.3575, 0.00, 0.38, 33.36, 8.3310, 0.0, 0.0],
+    'ethylphenylether':[1.5076, 1.5254, 0.00, 0.32, 46.65, 4.1797, 0.667, 0.0],
+    'ethylbenzene':[1.4959, 1.4932, 0.00, 0.15, 41.38, 2.4339, 0.75, 0.0],
+    'fluorobenzene':[1.4684, 1.4629, 0.00, 0.10, 38.37, 5.42, 0.857, 0.143],
+    'formamide':[1.4472, 1.4468, 0.62, 0.60, 82.08, 108.94, 0.0, 0.0],
+    'formicacid':[1.3714, 1.3693, 0.75, 0.38, 53.44, 51.1, 0.0, 0.0],
+    'n-heptane':[1.3878, 1.3855, 0.00, 0.00, 28.28, 1.9113, 0.0, 0.0],
+    'n-hexadecane':[1.4345, 1.4325, 0.00, 0.00, 38.93, 2.0402, 0.0, 0.0],
+    'n-hexane':[1.3749, 1.3722, 0.00, 0.00, 25.75, 1.8819, 0.0, 0.0],
+    'hexanoicacid':[1.4163, 1.4146, 0.60, 0.45, 39.65, 2.6, 0.0, 0.0],
+    'iodobenzene':[1.6200, 1.6172, 0.00, 0.12, 55.72, 4.5470, 0.857, 0.0],
+    'iodoethane':[1.5133, 1.5100, 0.00, 0.15, 40.96, 7.6177, 0.0, 0.0],
+    'iodomethane':[1.5380, 1.5270, 0.00, 0.13, 43.67, 6.8650, 0.0, 0.0],
+    'isopropylbenzene':[1.4915, 1.4889, 0.00, 0.16, 39.85, 2.3712, 0.667, 0.0],
+    'p-isopropyltoluene':[1.4909, 1.4885, 0.00, 0.19, 38.34, 2.2322, 0.600, 0.0],
+    'mesitylene':[1.4994, 1.4968, 0.00, 0.19, 39.65, 2.2650, 0.667, 0.0],
+    'methanol':[1.3288, 1.3265, 0.43, 0.47, 31.77, 32.613, 0.0, 0.0],
+    'methylbenzoate':[1.5164, 1.5146, 0.00, 0.46, 53.5, 6.7367, 0.600, 0.0],
+    'methylbutanoate':[1.3878, 1.3847, 0.00, 0.45, 35.44, 5.5607, 0.0, 0.0],
+    'methylethanoate':[1.3614, 1.3589, 0.00, 0.45, 35.59, 6.8615, 0.0, 0.0],
+    'methylmethanoate':[1.3433, 1.3415, 0.00, 0.38, 35.06, 8.8377, 0.0, 0.0],
+    'methylpropanoate':[1.3775, 1.3742, 0.00, 0.45, 35.18, 6.0777, 0.0, 0.0],
+    'N-methylaniline':[1.5684, 1.5681, 0.17, 0.43, 53.11, 5.9600, 0.75, 0.0],
+    'methylcyclohexane':[1.4231, 1.4206, 0.00, 0.00, 33.52, 2.024, 0.0, 0.0],
+    'N-methylformamide(E/Zmixture)':[1.4319, 1.4310, 0.40, 0.55, 55.44, 181.56, 0.0, 0.0],
+    'nitrobenzene':[1.5562, 1.5030, 0.00, 0.28, 57.54, 34.809, 0.667, 0.0],
+    'nitroethane':[1.3917, 1.3897, 0.02, 0.33, 46.25, 28.29, 0.0, 0.0],
+    'nitromethane':[1.3817, 1.3796, 0.06, 0.31, 52.58, 36.562, 0.0, 0.0],
+    'o-nitrotoluene':[1.5450, 1.5474, 0.0, 0.27, 59.12, 25.669, 0.6, 0.0],
+    'n-nonane':[1.4054, 1.4031, 0.0, 0.0, 32.21, 1.9605, 0.0, 0.0],
+    'n-octane':[1.3974, 1.3953, 0.0, 0.0, 30.43, 1.9406, 0.0, 0.0],
+    'n-pentadecane':[1.4315, 1.4298, 0.0, 0.0, 38.34, 2.0333, 0.0, 0.0],
+    'pentanal':[1.3944, 1.3917, 0.0, 0.4, 36.62, 10.0, 0.0, 0.0],
+    'n-pentane':[1.3575, 1.3547, 0.0, 0.0, 22.3, 1.8371, 0.0, 0.0],
+    'pentanoic acid':[1.4085, 1.4060, 0.60, 0.45, 38.4, 2.6924, 0.0, 0.0],
+    'pentyl ethanoate':[1.4023, -1.0, 0.0, 0.45, 36.23, 4.7297, 0.0, 0.0],
+    'pentylamine':[1.448, 1.4093, 0.16, 0.61, 35.54, 4.2010, 0.0, 0.0],
+    'perfluorobenzene':[1.3777, 1.3761, 0.00, 0.00, 31.74, 2.029, 0.5, 0.5],
+    'propanal':[1.3636, 1.3593, 0.00, 0.45, 32.48, 18.5, 0.0, 0.0],
+    'propanoic acid':[1.3869, 1.3848, 0.60, 0.45, 37.71, 3.44, 0.0, 0.0],
+    'propanonitrile':[1.3655, 1.3633, 0.02, 0.36, 38.5, 29.324, 0.0, 0.0],
+    'propyl ethanoate':[1.3842, 1.3822, 0.0, 0.45, 34.26, 5.5205, 0.0, 0.0],
+    'propylamine':[1.3870, 1.3851, 0.16, 0.61, 31.31, 4.9912, 0.0, 0.0],
+    'pyridine':[1.5095, 1.5073, 0.0, 0.52, 52.62, 12.978, 0.833, 0.0],
+    'tetrachloroethene':[1.5053, 1.5055, 0.0, 0.0, 45.19, 2.268, 0.0, 0.667],
+    'tetrahydrofuran':[1.4050, 1.4044, 0.0, 0.48, 39.44, 7.4257, 0.0, 0.0],
+    'tetrahydrothiophene-S,S-dioxide':[1.4833, -1.0, 0.0, 0.88, 87.49, 43.962, 0.0, 0.0],
+    'tetralin':[1.5413, 1.5392, 0.0, 0.19, 47.74, 2.771, 0.6, 0.0],
+    'thiophene':[1.5289, 1.5268, 0.0, 0.15, 44.16, 2.7270, 0.8, 0.0],
+    'thiophenol':[1.5893, 1.580, 0.09, 0.16, 55.24, 4.2728, 0.857, 0.0],
+    'toluene':[1.4961, 1.4936, 0.0, 0.14, 40.2, 2.3741, 0.857, 0.0],
+    'trans-decalin':[1.4695, 1.4671, 0.0, 0.0, 42.19, 2.1781, 0.0, 0.0],
+    'tributylphosphate':[1.4224, 1.4215, 0.0, 1.21, 27.55, 8.1781, 0.0, 0.0],
+    'trichloroethene':[1.4773, 1.4556, 0.08, 0.03, 41.45, 3.422, 0.0, 0.6],
+    'triethylamine':[1.4010, 1.3980, 0.0, 0.79, 29.1, 2.3832, 0.0, 0.0],
+    'n-undecane':[1.4398, 1.4151, 0.0, 0.0, 34.85, 1.991, 0.0, 0.0],
+    'water':[1.3328, 1.3323, 0.82, 0.35, -1.0, 78.355, -1.0, -1.0],
+    'xylene (mixture)':[1.4995, 1.4969, 0.0, 0.16, 41.38, 2.3879, 0.75, 0.0],
+    'm-xylene':[1.4972, 1.4946, 0.0, 0.16, 40.98, 2.3478, 0.75, 0.0],
+    'o-xylene':[1.5055, 1.5029, 0.0, 0.16, 42.83, 2.5454, 0.75, 0.0],
+    'p-xylene':[1.4958, 1.4932, 0.0, 0.16, 40.32, 2.2705, 0.75, 0.0],
+    '':[0]*8
+}
+
+def smd_radii(alpha):
+    '''
+    eq. (16)
+    use smd radii if defined
+    use Bondi radii if defined
+    use 2.0 otherwise
+    '''
+    radii_table = radii.VDW.copy() * radii.BOHR
+    radii_table[1] = 1.20
+    radii_table[6] = 1.85
+    radii_table[7] = 1.89
+    if alpha >= 0.43:
+        r = 1.52
+    else:
+        r = 1.52 + 1.8 * (0.43 - alpha)
+    radii_table[8] = r
+    radii_table[9] = 1.73
+    radii_table[14] = 2.47
+    radii_table[15] = 2.12
+    radii_table[16] = 2.49
+    radii_table[17] = 2.38
+    #radii_table[35] = 3.06 # original SMD
+    # following value from SMD18
+    # https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/chem.201803652
+    radii_table[35] = 2.60
+    radii_table[53] = 2.74
+    return radii_table/radii.BOHR
+
+def swtich_function(R, r, dr):
+    return np.exp(dr/(R-dr-r)) if R<r+dr else 0
+
+def atomic_surface_tension(symbols, coords, n, alpha, beta, water=True):
+    '''
+    - list of atomic symbols
+    - atomic coordinates in Anstrong
+    - solvent descriptors: n, alpha, beta
+    '''
+
+    def get_bond_tension(bond):
+        if water:
+            return sigma_water.get(bond, 0.0)
+        t = 0.0
+        t += sigma_n.get(bond, 0.0) * n
+        t += sigma_alpha.get(bond, 0.0) * alpha
+        t += sigma_beta.get(bond, 0.0) * beta
+        return t
+
+    def get_atom_tension(sym_i):
+        if water:
+            return sigma_water.get(sym_i, 0.0)
+        t = 0.0
+        t += sigma_n.get(sym_i, 0.0) * n
+        t += sigma_alpha.get(sym_i, 0.0) * alpha
+        t += sigma_beta.get(sym_i, 0.0) * beta
+        return t
+
+    rij = cupy.sum((coords[:,None,:] - coords[None,:,:])**2, axis=2)**0.5
+    tensions = []
+    for i, sym_i in enumerate(symbols):
+        if sym_i not in ['H', 'C', 'N', 'O', 'F', 'Si', 'S', 'Cl', 'Br']:
+            tensions.append(0)
+            continue
+
+        #sig_n = sigma_n.get(sym_i, 0.0)
+        #sig_a = sigma_alpha.get(sym_i, 0.0)
+        #sig_b = sigma_beta.get(sym_i, 0.0)
+        #tension = sig_n * n + sig_a * alpha + sig_b * beta
+        tension = get_atom_tension(sym_i)
+        if sym_i in ['F', 'Si', 'S', 'Cl', 'Br']:
+            tensions.append(tension)
+            continue
+
+        if sym_i == 'H':
+            t_HC = 0.0
+            t_HO = 0.0
+            for j, sym_j in enumerate(symbols):
+                if sym_j == 'C':
+                    r, dr = r_zz['H','C']
+                    t_HC += swtich_function(rij[i,j], r, dr)
+                if sym_j == 'O':
+                    r, dr = r_zz['H','O']
+                    t_HO += swtich_function(rij[i,j], r, dr)
+            sig_HC = get_bond_tension(('H','C'))
+            sig_HO = get_bond_tension(('H','O'))
+            tension += sig_HC * t_HC + sig_HO * t_HO
+            tensions.append(tension)
+
+        if sym_i == 'C':
+            t_CC = 0.0
+            t_CN = 0.0
+            for j, sym_j in enumerate(symbols):
+                if sym_j == 'C' and i != j:
+                    r, dr = r_zz['C', 'C']
+                    t_CC += swtich_function(rij[i,j], r, dr)
+                if sym_j == 'N':
+                    r, dr = r_zz['C', 'N']
+                    t_CN += swtich_function(rij[i,j], r, dr)
+            sig_CC = get_bond_tension(('C','C'))
+            sig_CN = get_bond_tension(('C','N'))
+            tension += sig_CC * t_CC + sig_CN * t_CN**2
+            tensions.append(tension)
+
+        if sym_i == 'N':
+            t_NC = 0.0
+            t_NC3 = 0.0
+            for j, sym_j in enumerate(symbols):
+                if sym_j == 'C':
+                    r, dr = r_zz['N','C']
+                    tk = 0.0
+                    for k, sym_k in enumerate(symbols):
+                        if k != i and k != j:
+                            rjk, drjk = r_zz['C', sym_k]
+                            tk += swtich_function(rij[j,k], rjk, drjk)
+                    t_NC += swtich_function(rij[i,j], r, dr) * tk**2
+
+                    r, dr = r_zz['N','C3']
+                    t_NC3 += swtich_function(rij[i,j], r, dr)
+            sig_NC = get_bond_tension(('N','C'))
+            sig_NC3= get_bond_tension(('N','C3'))
+            tension += sig_NC * t_NC**1.3 + sig_NC3 * t_NC3
+            tensions.append(tension)
+
+        if sym_i == 'O':
+            t_OC = 0.0
+            t_ON = 0.0
+            t_OO = 0.0
+            t_OP = 0.0
+            for j, sym_j in enumerate(symbols):
+                if sym_j == 'C':
+                    r, dr = r_zz['O','C']
+                    t_OC += swtich_function(rij[i,j], r, dr)
+                if sym_j == 'N':
+                    r, dr = r_zz['O','N']
+                    t_ON += swtich_function(rij[i,j], r, dr)
+                if sym_j == 'O' and j != i:
+                    r, dr = r_zz['O','O']
+                    t_OO += swtich_function(rij[i,j], r, dr)
+                if sym_j == 'P':
+                    r, dr = r_zz['O','P']
+                    t_OP += swtich_function(rij[i,j], r, dr)
+            sig_OC = get_bond_tension(('O','C'))
+            sig_ON = get_bond_tension(('O','N'))
+            sig_OO = get_bond_tension(('O','O'))
+            sig_OP = get_bond_tension(('O','P'))
+            tension += sig_OC * t_OC + sig_ON * t_ON + sig_OO * t_OO + sig_OP * t_OP
+            tensions.append(tension)
+    return cupy.asarray(tensions)
+
+def molecular_surface_tension(beta, gamma, phi, psi):
+    sig_gamma = sigma_gamma * gamma / gamma0
+    sig_phi = sigma_phi2 * phi**2
+    sig_psi = sigma_psi2 * psi**2
+    sig_beta= sigma_beta2 * beta**2
+    return sig_gamma + sig_phi + sig_psi + sig_beta
+
+def naive_sasa(mol, rad):
+    coords = mol.atom_coords(unit='A')
+    charges = mol.atom_charges()
+    radius = [rad[ch] for ch in charges]
+    sasa = []
+    for i in range(mol.natm):
+        area = 4 * np.pi * radius[i]
+        for j in range(mol.natm):
+            if i != j:
+                dr = coords[i] - coords[j]
+                r = (dr[0]**2 + dr[1]**2 + dr[2]**2)**0.5
+                r1 = radius[i]
+                r2 = radius[j]
+                if r < r1 + r2:
+                    overlap = (r1 + r2 - r) / (r1 + r2)
+                    area -= overlap * area
+        sasa.append(area)
+    return cupy.asarray(sasa)
+
+def get_cds(smdobj):
+    mol = smdobj.mol
+    n, _, alpha, beta, gamma, _, phi, psi = smdobj.solvent_descriptors
+    symbols = [mol.atom_symbol(ia) for ia in range(mol.natm)]
+    coords = mol.atom_coords(unit='A')
+    if smdobj._solvent.lower() != 'water':
+        atm_tension = atomic_surface_tension(symbols, coords, n, alpha, beta, water=False)
+        mol_tension = molecular_surface_tension(beta, gamma, phi, psi)
+    else:
+        logger.info(mol, 'no solvent descriptor is needed for water')
+        atm_tension = atomic_surface_tension(symbols, coords, n, alpha, beta, water=True)
+        mol_tension = 0.0
+
+    # generate surface for calculating SASA
+    rad = radii.VDW + 0.4/radii.BOHR
+    surface = pcm.gen_surface(mol, ng=302, rad=rad)
+    area = surface['area']
+    gridslice = surface['gslice_by_atom']
+    SASA = cupy.asarray([cupy.sum(area[p0:p1], axis=0) for p0,p1, in gridslice])
+    SASA *= radii.BOHR**2
+    mol_cds = mol_tension * cupy.sum(SASA) / 1000 # in kcal/mol
+    atm_cds = cupy.sum(SASA * atm_tension) / 1000 # in kcal/mol
+
+    return (mol_cds + atm_cds)/hartree2kcal # hartree
+
+class SMD(pcm.PCM):
+    _keys = {
+        'intopt', 'method', 'e_cds', 'solvent_descriptors', 'r_probe'
+    }
+    def __init__(self, mol, solvent=''):
+        super().__init__(mol)
+        self.vdw_scale = 1.0
+        self.r_probe = 0.4/radii.BOHR
+        self.method = 'SMD'  # use IEFPCM for electrostatic
+        if solvent not in solvent_db:
+            raise RuntimeError(f'{solvent} is not available in SMD')
+        self._solvent = solvent
+        self.solvent_descriptors = solvent_db[solvent]
+        self.radii_table = smd_radii(self.solvent_descriptors[2])
+        self.e_cds = None
+
+    @property
+    def solvent(self):
+        return self._solvent
+
+    @solvent.setter
+    def solvent(self, solvent):
+        self._solvent = solvent
+        self.solvent_descriptors = solvent_db[solvent]
+        self.radii_table = smd_radii(self.solvent_descriptors[2])
+        self.eps = self.solvent_descriptors[5]
+
+    @property
+    def sol_desc(self):
+        return self.solvent_descriptors
+
+    @sol_desc.setter
+    def sol_desc(self, values):
+        '''
+        format of sol desc
+        [n, n25, alpha, beta, gamma, epsilon, phi, psi]
+        '''
+        assert len(values) == 8
+        self.solvent_descriptors = values
+        self.radii_table = smd_radii(self.solvent_descriptors[2])
+        self.eps = values[5]
+
+    def dump_flags(self, verbose=None):
+        n, _, alpha, beta, gamma, _, phi, psi = self.solvent_descriptors
+        logger.info(self, '******** %s ********', self.__class__)
+        logger.info(self, 'lebedev_order = %s (%d grids per sphere)',
+                    self.lebedev_order, gen_grid.LEBEDEV_ORDER[self.lebedev_order])
+        logger.info(self, 'lmax = %s'         , self.lmax)
+        logger.info(self, 'eta = %s'          , self.eta)
+        logger.info(self, 'eps = %s'          , self.eps)
+        logger.info(self, 'frozen = %s'       , self.frozen)
+        logger.info(self, '---------- SMD solvent descriptors -------')
+        logger.info(self, f'n     = {n}')
+        logger.info(self, f'alpha = {alpha}')
+        logger.info(self, f'beta  = {beta}')
+        logger.info(self, f'gamma = {gamma}')
+        logger.info(self, f'phi   = {phi}')
+        logger.info(self, f'psi   = {psi}')
+        logger.info(self, '--------------------- end ----------------')
+        logger.info(self, 'equilibrium_solvation = %s', self.equilibrium_solvation)
+        logger.info(self, 'radii_table %s', self.radii_table*radii.BOHR)
+        if self.atom_radii:
+            logger.info(self, 'User specified atomic radii %s', str(self.atom_radii))
+        self.grids.dump_flags(verbose)
+        return self
+
+    def get_cds(self):
+        return get_cds(self)
+
+    def nuc_grad_method(self, grad_method):
+        raise RuntimeError('SMD gradient is not implemented')
+
+    def Hessian(self, hess_method):
+        raise RuntimeError('SMD Hessian is not implemented')
+
+    def reset(self, mol=None):
+        super().reset(mol)
+        self.e_cds = None
+        return self
\ No newline at end of file
diff --git a/gpu4pyscf/solvent/tests/test_smd.py b/gpu4pyscf/solvent/tests/test_smd.py
new file mode 100644
index 00000000..fca5f64f
--- /dev/null
+++ b/gpu4pyscf/solvent/tests/test_smd.py
@@ -0,0 +1,69 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import unittest
+import numpy
+from pyscf import gto, df
+from gpu4pyscf import scf
+from gpu4pyscf.solvent import smd
+
+def setUpModule():
+    global mol, epsilon, lebedev_order
+    mol = gto.Mole()
+    mol.atom = '''
+O       0.0000000000    -0.0000000000     0.1174000000
+H      -0.7570000000    -0.0000000000    -0.4696000000
+H       0.7570000000     0.0000000000    -0.4696000000
+    '''
+    mol.basis = 'def2-tzvpp'
+    #mol.output = '/dev/null'
+    mol.build()
+    lebedev_order = 29
+
+def tearDownModule():
+    global mol
+    mol.stdout.close()
+    del mol
+
+class KnownValues(unittest.TestCase):
+    def test_cds_solvent(self):
+        smdobj = smd.SMD(mol)
+        smdobj.solvent = 'toluene'
+        e_cds = smdobj.get_cds()
+        assert numpy.abs(e_cds - -0.0013476060879874362) < 1e-8
+
+    def test_cds_water(self):
+        smdobj = smd.SMD(mol)
+        smdobj.solvent = 'water'
+        e_cds = smdobj.get_cds()
+        assert numpy.abs(e_cds - 0.0022847142144050057) < 1e-8
+
+    def test_smd_solvent(self):
+        mf = scf.RHF(mol)
+        mf = mf.SMD()
+        mf.with_solvent.solvent = 'ethanol'
+        e_tot = mf.kernel()
+        assert numpy.abs(e_tot - -76.075066568) < 2e-4
+
+    def test_smd_water(self):
+        mf = scf.RHF(mol)
+        mf = mf.SMD()
+        mf.with_solvent.solvent = 'water'
+        e_tot = mf.kernel()
+        assert numpy.abs(e_tot - -76.0756052903) < 2e-4
+
+if __name__ == "__main__":
+    print("Full Tests for SMDs")
+    unittest.main()
\ No newline at end of file