def METHOD_NAME(self):
    from pyscf import qmmm
    mol = gto.Mole()
    mol.atom = ''' O                  0.00000000    0.00000000   -0.11081188
                   H                 -0.00000000   -0.84695236    0.59109389
                   H                 -0.00000000    0.89830571    0.52404783 '''
    mol.verbose = 0
    mol.basis = '6-31g'
    mol.build()
    coords = [(0.5,0.6,0.1)]
    #coords = [(0.0,0.0,0.0)]
    charges = [-0.1]
    mf = qmmm.add_mm_charges(mol.RHF, coords, charges)
    ps = mf.MP2().as_scanner()
    g = ps.nuc_grad_method().as_scanner()(mol)[1]
    e1 = ps(''' O                  0.00100000    0.00000000   -0.11081188
             H                 -0.00000000   -0.84695236    0.59109389
             H                 -0.00000000    0.89830571    0.52404783 ''')
    e2 = ps(''' O                 -0.00100000    0.00000000   -0.11081188
             H                 -0.00000000   -0.84695236    0.59109389
             H                 -0.00000000    0.89830571    0.52404783 ''')
    ref = (e1 - e2)/0.002 * lib.param.BOHR
    self.assertAlmostEqual(g[0,0], ref, 4)