Calculating hyperpolarizability in NWChem

Question

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

task dft energy

Answer 1

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

task dft property

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

task tce energy

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.