Nuclear gradients for PolEmbed

[Unathi-Skosana]

I'm trying to optimize a molecular geometry within the polarizable embedding model

import pyscf
from pyscf import solvent

mol = pyscf.M(atom='''
C        0.000000    0.000000             -0.542500
O        0.000000    0.000000              0.677500
H        0.000000    0.9353074360871938   -1.082500
H        0.000000   -0.9353074360871938   -1.082500
''')

pe_options = {
    "max_iter": 50,
    "potfile":  "4NP_in_water/4NP_in_water.pot"
}
pe = solvent.pol_embed.PolEmbed(mol, pe_options)
pe.equilibrium_solvation = True
pe.verbose = 4

mf = mol.UKS()
mf.xc = 'b3lyp'
mf_with_pe = solvent.PE(mf, pe)
mf_with_pe.kernel()
new_mol_eq = mf_with_pe.Gradients().optimizer().kernel()

Looks like Gradients aren't yet implemented for PolEmbed. Will this, by any chance, be implemented in the future?