def METHOD_NAME(ia, inc):
    coord = mol.atom_coord(ia).copy()
    ptr = mol._atm[ia,gto.PTR_COORD]
    de = []
    for i in range(3):
        mol._env[ptr+i] = coord[i] + inc
        e1a = mfs(mol.copy())[1]
        mol._env[ptr+i] = coord[i] - inc
        e1b = mfs(mol.copy())[1]
        mol._env[ptr+i] = coord[i]
        de.append((e1a-e1b)/(2*inc))
    return de