I was trying to calculate the hyperpolarizability $\beta$ of a cluster of two water molecules in Gaussian16 using the route card:
#t cam-b3lyp/daug-cc-pvtz Polar nosymm scf=tight int=grid=superfine
#t cam-b3lyp/daug-cc-pvtz Polar=doubleNumer nosymm scf=tight int=grid=superfine
The files which can be used for testing are attached in these links [1], [2].
In case (1) I got a value of 83.44 au for the $\beta_{xxx}$. For (2) I got a value of 23.66 au for $\beta_{xxx}$. I also altered the distance between two water molecules by $\pm$ 0.5 angstroms and I got an average $\beta_{xxx}$ value of 23.44 au .
I am not sure about the capabilities of other electronic structure codes (I know ORCA does not have an option to calculate hyperpolarizability; Dalton Turbomole does, but I don't have access to it), but is this error reproducible in other codes or this is just a problem of G16?
From Gaussian16 support, I got the following response:
We have reviewed this further and it is a numerical error which is caused by near linear dependencies in the basis set which the analytic second derivative for beta is not handling well. You can see this if you use aug-cc-pvtz in place of daug-cc-pvtz.
It is also possible to get a better analytic value using Polar=(Cubic,Fourpoint) which is more numerically stable but also significantly more time consuming computationally. With that method the results are:
#p cam-b3lyp/daug-cc-pvtz Polar=(Cubic,fourpoint) nosymm scf=conver=10 int=grid=superfine
$\beta_{xxx}$=23.473904
#p cam-b3lyp/daug-cc-pvtz Polar=DoubleNumer nosymm scf=Conver=10 int=grid=superfine
$\beta_{xxx}$=23.5511654
where you see the comparison with numerical is about 5 significant figures.