11
$\begingroup$

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
$\endgroup$
6
  • $\begingroup$ Please share the full input so I can reproduce. Thanks. $\endgroup$ Nov 30 '20 at 4:06
  • $\begingroup$ 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. $\endgroup$ Nov 30 '20 at 4:08
  • $\begingroup$ I've added full input for MeCN, the choice of molecule, basis and functional is arbitrary. $\endgroup$ Nov 30 '20 at 10:04
  • $\begingroup$ It’s not the cause of the issue but 6-31G* is likely not going to provide meaningful results for hyperpolarizabilities. $\endgroup$ Nov 30 '20 at 15:30
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 30 '20 at 19:43
8
$\begingroup$

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: CH3CN hyperpolarizabilities with DFT and CC methods CH3CN hyperpolarizabilities with CC and many basis sets

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.

$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$ Dec 1 '20 at 3:53
  • $\begingroup$ 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. $\endgroup$ Dec 1 '20 at 4:02
  • $\begingroup$ 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. $\endgroup$ Dec 1 '20 at 4:06
  • $\begingroup$ 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? $\endgroup$ Dec 1 '20 at 4:13
  • $\begingroup$ In benchmarks that I’ve seen m06 and long range corrected functionals perform roughly equally using mp2 or experimental data as a refernce. $\endgroup$ Dec 1 '20 at 6:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.