I need to calculate interaction energy for a dimer using (mid)bond functions
which is known to give results closer to the convergence limit.
First I took a simple example of Ar2 dimer with aug-cc-pvdz basis
and tried to reproduce the results by comparison with other codes.
It worked well when there were no bond functions used, but did NOT work
with the bond functions used.
Here are the inputs
$molecule
0 1
@Ar 0 0 0.000000000
@Be 0 0 1.875
Ar 0 0 3.75
$end
$basis
Ar 1
aug-cc-pvdz
****
Be 2
S 3 1.00
0.9 1.000000000
0.3 1.000000000
0.1 1.000000000
P 3 1.00
0.9 1.000000000
0.3 1.000000000
0.1 1.000000000
D 2 1.00
0.6 1.000000000
0.2 1.000000000
F 1 1.00
0.3000000 1.00000000
G 1 1.00
0.25000000 1.00000000
****
Ar 3
aug-cc-pvdz
****
$end
$rem
BASIS gen
METHOD ccsd(t)
N_FROZEN_CORE fc
PRINT_ORBITALS
$end
$molecule
0 1
Ar 0 0 0.000000000
@Be 0 0 1.875
Ar 0 0 3.75
$end
$basis
Ar 1
aug-cc-pvdz
****
Be 2
S 3 1.00
0.9 1.000000000
0.3 1.000000000
0.1 1.000000000
P 3 1.00
0.9 1.000000000
0.3 1.000000000
0.1 1.000000000
D 2 1.00
0.6 1.000000000
0.2 1.000000000
F 1 1.00
0.3000000 1.00000000
G 1 1.00
0.25000000 1.00000000
****
Ar 3
aug-cc-pvdz
****
$end
$rem
BASIS gen
METHOD ccsd(t)
N_FROZEN_CORE fc
PRINT_ORBITALS
$end
What is wrong in my input? I tried to used Gh instead of @Be, tried both
mixed and gen basis, but it did NOT help.
Another question raised by me in another topic, is there any way to
increase the number of digits in the output, e.g. in this case
CCSD(T) total energy = -1053.94949414
8 decimal places in this case means an error of 0.2 cm-1 for Ar2 interaction
energy at the equilibrium configuration which I cannot ignore.
$molecule
0 1
Ar 0 0 0.000000000
GH 0 0 1.875
Ar 0 0 3.75
$end
$rem
method hf
basis mixed
print_general_basis true
$end
$basis
Ar 1
aug-cc-pvdz
****
Gh 2
S 3 1.00
0.9 1.000000000
0.3 1.000000000
0.1 1.000000000
P 3 1.00
0.9 1.000000000
0.3 1.000000000
0.1 1.000000000
D 2 1.00
0.6 1.000000000
0.2 1.000000000
F 1 1.00
0.3000000 1.00000000
G 1 1.00
0.25000000 1.00000000
****
Ar 3
aug-cc-pvdz
****
$end
PRINT_GENERAL_BASIS = TRUE regurgitates the basis info in the output and it looks like what I intended, therefore seems to work.
As for the digits in CCSD(T) total energy, I’m not sure this can be increased but I also question whether 1e-9 Hartree differences (or 13 sig figs in total energy) are meaningful.
Why is it that you think the original is not working? What I was saying is that with my input file, the basis info that was returned seems like what I was expecting, including the mid-bond functions. If I do the calculation with or without that ghost atom, I do get a different Hartree-Fock energy:
-1053.6024059317 (with bond functions)
-1053.6015769549 (without)
As to the accuracy, my bad, it’s 0.001 cm-1 indeed. Anyway, more digits might be important e.g. for finite-difference formulas when calculating anharmonic force constants of PES.
Yes, I confirm that midbond functions are working in Q-Chem. The problem of disagreement with other codes was that I didn’t define the midbond basis set properly.
It is supposed to be
S 1 1.00
0.9000000 1.0
S 1 1.00
0.3000000 1.0
S 1 1.00
0.1000000 1.0
P 1 1.00
0.9000000 1.0
P 1 1.00
0.3000000 1.0
P 1 1.00
0.1000000 1.0
D 1 1.00
0.6000000 1.0
D 1 1.00
0.2000000 1.0
F 1 1.00
0.30000 1.0
G 1 1.00
0.25000 1.0
I’m not sure though in which cases I have to use the “purecard” option for compatibility. Thanks for your help.
It’s hard to say about PURECART because I don’t know what other codes use as defaults (and Q-Chem defaults are different for different basis sets). The good news is that it likely doesn’t matter very much, only in cases where you are trying to reproduce a different code to all digits is this important. For correlated wave function methods (as opposed to DFT), I would set PURECART=1111, meaning to use spherical harmonics for all angular momentum classes. The extra s-type function that you get from using Cartesian d functions doesn’t buy you very much variational flexibility but it does create one additional virtual orbital, and CCSD (for example) scales as o^2 v^4. For DFT, you might as well use PURECART=2222 because the code always evaluates integrals using Cartesian functions (6d not 5d), so even if it doesn’t buy you much in the way of variational flexibility it doesn’t cost you anything (in DFT) to leave that extra function in there.