Dear all.
I am looking for a way to calculate net neutral molecules with a non-integer nuclear charge for a perturbative approach. So I am looking for self-consistent electron densities calculated as if a subset of the nuclei had a different charge, e.g. N2 with the basis set of N2 but with nuclear charges (6.95, 7.05). This is similar in spirit to a question [1] asked here before, which states it can be done, but I don’t see a way how to do it.
This is what I tried:
A) QMMM interface to add a point charge at the nuclei (example input below)
The problem is that all nuclear interactions are included in the total energy, also the ones between the newly added point charge and the other atoms. Since the point charges are at exactly the position of the nuclei, the divergence gives a total energy of nan (and does not converge).
B) QMMM interface with point charge at the nuclei but ghost atom (example input below)
Still same issue even though the ghost atom should not interact electrostatically. Wouldn’t this consititute a bug in the QMMM code?
Does anybody know how to calculate this case in psi4 properly?
Thank you for your time,
Guido
[1] Use fractional charge in calculation
Example inputs:
A)
molecule mol {
N 0. 0. 1.
N 0. 0. 0.
no_reorient
no_com
symmetry c1
}
Chrgfield = QMMM()
Chrgfield.extern.addCharge(0.05, 0., 0., 1.)
psi4.set_global_option_python('EXTERN', Chrgfield.extern)
set basis 6-31G
set reference rhf
e, wfn = energy('scf', return_wfn=True)
B)
molecule mol {
@N 0. 0. 1.
N 0. 0. 0.
no_reorient
no_com
symmetry c1
-7 1
}
Chrgfield = QMMM()
Chrgfield.extern.addCharge(7.05, 0., 0., 1.)
psi4.set_global_option_python('EXTERN', Chrgfield.extern)
set basis 6-31G
set reference rhf
e, wfn = energy('scf', return_wfn=True)