# Calculating hyperpolarizability in NWChem

I'm trying to make a pipeline geometry optimization -> calculation of polarizability/hyperpolarizability. And NWChem seems pretty good for it, the only thing I can't understand its how to calculate hyperpolarizability.

Currently I'm using this input, but it gives me only polarizability (alpha as 3x3 matrix) and I need hyperpolarizability (beta as 3x3x3 tensor).

Full input:

echo

start molecule

title "MeCN"
charge 0

geometry units angstroms print xyz autosym
C       -6.60969        1.24918       -0.00000
C       -6.19457       -0.15256       -0.00000
H       -7.63700        1.35035       -0.36238
H       -5.96068        1.84678       -0.64699
H       -6.56251        1.66855        1.00937
N       -5.86512       -1.26505        0.00000
end

basis
* library 6-31G*
end

dft
xc b3lyp
mult 1
end

property
response  2 7.73178E-2   # response order and frequency in Hartree energy units
end


• Please share the full input so I can reproduce. Thanks. – Jeff Hammond Nov 30 '20 at 4:06
• I’m surprised by this issue, by the way, so I’m trying to figure out if it’s a bug or if beta is zero by symmetry and that’s why it’s omitted. – Jeff Hammond Nov 30 '20 at 4:08
• I've added full input for MeCN, the choice of molecule, basis and functional is arbitrary. – roma ichenko Nov 30 '20 at 10:04
• It’s not the cause of the issue but 6-31G* is likely not going to provide meaningful results for hyperpolarizabilities. – Jeff Hammond Nov 30 '20 at 15:30
• in NWChem, there is no difference between static and dynamic, because we use the response formalism. The CC response module does not include orbital response, which is critical for geometry derivatives (e.g. gradients) but of limited value for electronic response properties, and is only valid for omega=0 anyways. – Jeff Hammond Nov 30 '20 at 19:43

# Major Issue

The major issue here is that NWChem doesn't support second-order response in the DFT code: "Response calculations are currently supported only for order 1 (linear response), single frequency, electric field and mixed electric-magnetic field perturbations." [1]

The NWChem coupled-cluster response property capability supports hyperpolarizabilities (I wrote it) and I guess that you got the impression the DFT code supported this too from the paper:

J. R. Hammond and K. Kowalski, J. Chem. Phys. 130, 194108 (2008). "Parallel computation of coupled-cluster hyperpolarizabilities." [2]

However, that paper used Dalton for DFT hyperpolarizabilities (see the end of Section II). The Dalton response property code was amazing and the primary motivation for using NWChem is to compute coupled-cluster response properties in parallel and/or with ROHF/UHF references (because NWChem TCE uses the spin-orbital formalism).

If it helps at all, below are the results for methylcyanide from that paper:

# Minor Issue

There is a relatively simple issue in the input file above. You need to specify the property task rather than the energy task for the input parser to look at the property .. end block.

echo

start molecule

title "MeCN"
charge 0

geometry units angstroms print xyz autosym
C       -6.60969        1.24918       -0.00000
C       -6.19457       -0.15256       -0.00000
H       -7.63700        1.35035       -0.36238
H       -5.96068        1.84678       -0.64699
H       -6.56251        1.66855        1.00937
N       -5.86512       -1.26505        0.00000
end

basis
* library 6-31G*
end

dft
xc b3lyp
mult 1
end

property
response  2 7.73178E-2   # response order and frequency in Hartree energy units
end



# CCSD Results

It isn't going to be practical for larger molecules, but I was able to obtain the CCSD result using NWChem just now:

echo

start CH3CN_ccsd_aug-cc-pvdz_hyperpolar

permanent_dir .
scratch_dir .

memory stack 1600 mb heap 200 mb global 1400 mb noverify

geometry units bohr
symmetry cs
N      2.5028407872      0.0000000000      0.0000000000
C      0.3138164240      0.0000000000      0.0000000000
C     -2.4431799386      0.0000000000      0.0000000000
H     -3.1404599851      0.9650062580     -1.6714398685
H     -3.1404599851      0.9650062580      1.6714398685
H     -3.1404599851     -1.9300125160      0.0000000000
end

basis spherical
* library aug-cc-pvdz
end

#set lindep:tol 1d-9
set int:acc_std 1d-14

scf
singlet
rhf
direct
end

tce
scf
freeze atomic
ccsd
maxiter 100
thresh 1.0e-7
io ga
2eorb
2emet 13
tilesize 18
attilesize 30
end

set tce:lineresp T
set tce:leftresp T
set tce:respaxis T T T
set tce:afreq 0.0
set tce:bfreq 0.0



This is the output:

CCSD Quadratic Response
Static Hyperpolarizability
-----------------------------------------------
beta(X,X,X) =      -34.504019056 / au
beta(X,Y,Y) =        3.997809593 / au
beta(X,Z,Z) =        3.997810512 / au
beta(Y,X,Y) =        3.997809593 / au
beta(Y,Y,X) =        3.997809593 / au
beta(Y,Y,Y) =        4.544908021 / au
beta(Y,Z,Z) =       -4.544912885 / au
beta(Z,X,Z) =        3.997810512 / au
beta(Z,Y,Z) =       -4.544912885 / au
beta(Z,Z,X) =        3.997810512 / au
beta(Z,Z,Y) =       -4.544912885 / au
-----------------------------------------------


This result was obtained using my Mac laptop with one core in approximately 10 minutes using a binary compiled with GCC 10 and the NWChem internal BLAS, which is the slowest but most reliable build. Larger basis sets will increase the computational requirements quickly.

• Thank you a lot for your detailed answer! Unfortunately CCSD is too expensive to calculate hyperpolarizabilities for the molecules I need. What software would you recommend for this task? – roma ichenko Dec 1 '20 at 3:53
• Use Dalton for DFT response properties. I haven’t used it in a decade but it was solid then. But how big are your molecules? I could do para-nitroaniline with CCSD/d-aug-cc-pVTZ in 2008. Computers are a lot bigger now. – Jeff Hammond Dec 1 '20 at 4:02
• I'm going to calculate properties for organic chromophores with formula like C40 F10 and some nitrogens, oxygens and hydrogens. So the only option is something like M06-2X/aug-cc-pVTZ and maybe MP2. – roma ichenko Dec 1 '20 at 4:06
• MP2 and CC2 are not that useful here. MP2 doesn’t have a response function so you can’t get dynamic responses while CC2 is inaccurate. I’ve not need M06 benchmarked for hyperpolarizabilities. Is it better than the range-corrected hybrids like LC-PBE? – Jeff Hammond Dec 1 '20 at 4:13
• In benchmarks that I’ve seen m06 and long range corrected functionals perform roughly equally using mp2 or experimental data as a refernce. – roma ichenko Dec 1 '20 at 6:31