wB97X-V and wB97M-V setup

[stanos4]

Dear all,

I’m trying to set-up the aforementioned functionals to compare the results to CP2K. I try to follow this publication for settings of the functionals in Psi4:
DOI: 10.1021/acs.jpca.8b03345
( Density Functional Theory for Microwave Spectroscopy of Noncovalent Complexes: A Benchmark Study)

to reproduce Ar-Ar interaction energies of the original wB97X-V paper (and possibly other energies):
DOI: 10.1063/1.4952647
(ω B97M-V: A combinatorially optimized, range-separated hybrid, meta-GGA density functional with VV10 nonlocal correlation)

here is the Python script that I use:

import numpy as np
import psi4
psi4.set_num_threads(8)
templ = """
0 1
Ar 0.0 0.0 0.0
Ar XXXX 0.0 0.0
nocom
noreorient
"""
Ar1 = """
0 1
Ar 0.0 0.0 0.0
nocom
noreorient
"""
method = "WB97M-V/def2-TZVPD"
method = "WB97X-V/def2-QZVPD"
psi4.set_options({"basis_guess": 0, "full_hess_every": hessian})
psi4.set_options({'scf_type':'df',"cc_type": "df"})
psi4.set_options({"dft_spherical_points": 1454, "dft_radial_points": 500})
psi4.set_options({"damping_percentage": 40})
psi4.set_options({"dft_bs_radius_alpha": 1.3})
psi4.set_options({"g_convergence": "interfrag_tight"})
energies = []
for x in range(0, 31):
    cx = 3.25 + x*0.1
    Ar2 = templ.replace("XXXX", str(cx))
    mol = psi4.geometry(Ar2)
    e2,wfn = psi4.energy(method, molecule=mol, return_wfn=True)
    energies.append(e2 * psi4.constants.hartree2kcalmol)
energies = np.array(energies)
mol = psi4.geometry(Ar1)
e1,wfn = psi4.energy(method, molecule=mol, return_wfn=True)
e1 = e1 * psi4.constants.hartree2kcalmol
energies - 2*e1

the results are:
[0.26051278, -0.04387099, -0.22114193, -0.33644522, -0.39150897,
-0.40188769, -0.40330375, -0.00466373, -0.36437939, -0.34776137,
-0.30978748, -0.28900184, -0.25777778, -0.2334124 , -0.2097273 ,
-0.18432839, -0.17304105, -0.15599188, -0.14514117, -0.12967334,
-0.1173416 , -0.10737023, -0.099812 , -0.02796177, -0.10710261,
-0.07656742, -0.07045404, -0.06649192, -0.06076612, -0.0561549 ,
-0.05177564]

compared to the original wB97X-V paper, this minimum is roughly twice as low as the original one. Do you please have got suggestions what am I doing wrong?

using Psi4 v 1.3rc2 installed with anaconda.
Thank you.
Kind regards,
stanislav.

[loriab]

some of those options are irrelevant
{“basis_guess”: 0, “full_hess_every”: hessian,“cc_type”: “df”, “g_convergence”: “interfrag_tight”}

are you sure the reference you’re comparing to used density-fitting (sorry, haven’t followed the links)?

[stanos4]

as far sa I can see, those are the settings used in the “…A Benchmark Study” - in supporting info they attached python files as well. Switching these off gives ~the same results.

I’m interested in reproducing results from the original wB97X-V and wB97M-v papers. That means, I’m in principle interested in settings “as tight as possible”, but my experience with psi4 is limited, so I will be very happy for any suggestions. Thank you very much.

[hokru]

The settings are indeed very very tight for psi4. Much tighter than the Q-chem settings, I assume. Could you let us know what graph you want to reproduce? I dont see Ar2 in the paper.
Could there be a CP correction applied?

[jmisiewicz]

The ωB97M-V paper tells you what grid settings they used in Section III.

I’m not seeing the argon dimer in the paper either, only the benzene argon dimer.

[stanos4]

Problem solved - my basis set was evidently not sufficient. Shame on me and thank you all for your feedback. I’m also trying to recompute results with the VV10 grid used in the wb97x-v paper, but it takes some time. I will post the settings later on.

Attached is graph I’m trying to reproduce - from original wB97X-V paper.

[stanos4]

shall it be of any use to anybody, here are settings that seemingly reproduce results in the wB97X-V paper up to <0.05kcal/mol, for the Ar-Ar binding energy - the plot in the previous post. I’m not sure if I still miss some settings, but the DFT and VV10 grids should be as in the original paper, as well as the basis set.

import numpy as np
import psi4
psi4.set_num_threads(16)
psi4.set_memory("30 GB")
templ = """
0 1
Ar 0.0 0.0 0.0
Ar XXXX 0.0 0.0
nocom
noreorient
"""
Ar1 = """
0 1
Ar 0.0 0.0 0.0
nocom
noreorient
"""
#method = "WB97M-V/def2-TZVPD"
#method = "WB97X-V/def2-QZVPPD"
#method = "WB97X-V/def2-QZVPD"
method = "WB97X-V/aug-cc-pVQZ"
psi4.set_options({"basis_guess": 0, "full_hess_every": hessian})
psi4.set_options({"dft_spherical_points": 974, "dft_radial_points": 500})
psi4.set_options({"damping_percentage": 40})
psi4.set_options({"dft_bs_radius_alpha": 1.3})
psi4.set_options({"g_convergence": "interfrag_tight"})
psi4.set_options({"DFT_VV10_SPHERICAL_POINTS": 590, "DFT_VV10_RADIAL_POINTS": 99})

energies = []
for x in range(0, 31):
    cx = 3.25 + x*0.1
    Ar2 = templ.replace("XXXX", str(cx))
    mol = psi4.geometry(Ar2)
    e2,wfn = psi4.energy(method, molecule=mol, return_wfn=True)
    energies.append(e2 * psi4.constants.hartree2kcalmol)

energies = np.array(energies)
mol = psi4.geometry(Ar1)
e1,wfn = psi4.energy(method, molecule=mol, return_wfn=True)
e1 = e1 * psi4.constants.hartree2kcalmol
for id, xx in enumerate(energies - 2*e1):
    print("{:.2f} A  : {:6.3f} kcal/mol".format(3.25+id*0.1, xx))