Contents
Psi4: Ar$_2$ interaction energies
SAPT
Supermolecular: CCSD(T), MP2, HF
Method 1: CP correction within Psi4
This is a method I borrowed from the Psi4 Wiki.
<code | Ar2_CCSDt_scan_CPpsi4.psi4>
#! Example potential energy surface scan and CP-correction for Ar2
#Based on a Psi4 example from the wiki
#http://psicode.org/psi4manual/master/tutorial.html#analysis-of-intermolecular-interactions
memory 10 gb
molecule dimer {
Ar
--
Ar 1 R
}
Rvals=[3.0, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]
set basis aug-cc-pVDZ
set freeze_core True
# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
ecp = {}
for R in Rvals:
dimer.R = R
ecp[R] = energy('ccsd(t)', bsse_type = 'cp')
psi4.print_out("\n")
psi4.print_out("CP-corrected CCSD(T)/aug-cc-pVDZ interaction energies\n\n")
psi4.print_out(" R [Ang] E_int [kJ/mol] \n")
psi4.print_out("-----------------------------------------------------\n")
for R in Rvals:
e = ecp[R] * psi_hartree2kJmol
psi4.print_out(" %3.1f %10.6f\n" % (R, e))This can be run using:
$ psi4 -n 1 Ar2_CCSDt_scan_CPpsi4.psi4 Ar2_CCSDt_scan_CPpsi4_aDZ.out &
In this case, I have used only one processor and 10 GB of memory. The aug-cc-pVDZ basis is a small one, so this is OK. The output is rather large, but the important part is the bottom bit which looks like this:
<code | Ar2_CCSDt_scan_CPpsi4_aDZ.out>
CP-corrected CCSD(T)/aug-cc-pVDZ interaction energies
R [Ang] E_int [kJ/mol]
-----------------------------------------------------
3.0 12.559936
3.5 0.706474
3.8 -0.292682
4.0 -0.515145
4.2 -0.487592
4.5 -0.398161
4.8 -0.308207
5.0 -0.233998
5.5 -0.134115
6.0 -0.078839
*** Psi4 exiting successfully. Buy a developer a beer!Here are some results in larger basis sets. Here's the aug-cc-pVTZ basis. This time I used 4 cores for the calculation, but this was not necessary; a single core would have been OK as this system is small:
<code | Ar2_CCSDt_scan_CPpsi4_aTZ.out>
CP-corrected CCSD(T)/aug-cc-pVTZ interaction energies
R [Ang] E_int [kJ/mol]
-----------------------------------------------------
3.0 10.025281
3.5 -0.092815
3.8 -0.790347
4.0 -0.831467
4.2 -0.687915
4.5 -0.524758
4.8 -0.389159
5.0 -0.287042
5.5 -0.159122
6.0 -0.092021And here are the results for the aug-cc-pVQZ basis:
<code | Ar2_CCSDt_scan_CPpsi4_aQZ.out>
CP-corrected CCSD(T)/aug-cc-pVQZ interaction energies
R [Ang] E_int [kJ/mol]
-----------------------------------------------------
3.0 8.695992
3.5 -0.453211
3.8 -0.975822
4.0 -0.933543
4.2 -0.749688
4.5 -0.564690
4.8 -0.415487
5.0 -0.304355
5.5 -0.166692
6.0 -0.095592Important
Has the interaction energy converged with respect to basis set? It doesn't appear to have, has it? How do you tell? Check Helgaker's book for more on basis-set convergence.
Important
Compare these results with those from NWChem. Pay attention to the units!
Method 2: Explicit CP correction
Important
This method has now been corrected using a suggestion from Stephen Rayner. Charlie Woollard may also have mentioned the correction: We need to use a clean() command between different calculations or Psi4 messes up stuff.
I wanted to try out these calculations in another way: one in which the dimer and monomer calculations are performed explicitly. This has the advantage that we see what's going on, and that we can extract the HF and MP2 energies too. I could not get the CCSD(T) to work this way until I inserted the clean() command between calculations. More on this later. First the command file:
Ar2_CCSDt_scan_CPexplicit.psi4
<code | Ar2_CCSDt_scan_CPexplicit.psi4>
#! Example potential energy surface scan and CP-correction for Ar2
#Based on a Psi4 example from the wiki
#http://psicode.org/psi4manual/master/tutorial.html#analysis-of-intermolecular-interactions
memory 40 gb
molecule dimer {
Ar
Ar 1 R
}
#Define a monomer with a ghost atom. We need this for the CP correction.
#Only one monomer need be defined as, for this system, the other is
#identical by symmetry.
molecule mono {
Ar
@Ar 1 R
}
Rvals=[3.0, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]
# For testing only:
#Rvals=[3.0, 3.5]
#molecules=[dimer,mono]
set basis aug-cc-pVDZ
set freeze_core True
# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_hf = {} # for dimer HF energies
e_d_mp2 = {} # dimer MP2
e_d_cc = {} # dimer CCSD(T)
e_m_hf = {} # And same for the monomer
e_m_mp2 = {}
e_m_cc = {}
# The CCSD(T) lines have been commented out below as I could not get that to work.
for R in Rvals:
dimer.R = R
energy('ccsd(t)',molecule=dimer)
e_d_hf[R] = get_variable('SCF TOTAL ENERGY')
e_d_mp2[R] = get_variable('MP2 TOTAL ENERGY')
e_d_cc[R] = get_variable('CCSD(T) TOTAL ENERGY')
#Remove tmp files etc or else the mono jobs will fail.
clean()
mono.R = R
energy('ccsd(t)',molecule=mono)
e_m_hf[R] = get_variable('SCF TOTAL ENERGY')
e_m_mp2[R] = get_variable('MP2 TOTAL ENERGY')
e_m_cc[R] = get_variable('CCSD(T) TOTAL ENERGY')
#And another clean or else the next job will fail.
clean()
psi4.print_out("\n")
psi4.print_out("CP-corrected HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies\n\n")
psi4.print_out(" R [Ang] HF MP2 CCSD(T) [kJ/mol] \n")
psi4.print_out("--------------------------------------------------------------\n")
for R in Rvals:
e_hf = (e_d_hf[R] -2.0*e_m_hf[R]) * psi_hartree2kJmol
e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2[R]) * psi_hartree2kJmol
e_cc = (e_d_cc[R] -2.0*e_m_cc[R]) * psi_hartree2kJmol
psi4.print_out(" %3.1f %10.6f %10.6f %10.6f \n" % (R, e_hf,e_mp2,e_cc))
Results for HF & MP2 in various basis sets
<code | aug-vv-pVDZ basis>
CP-corrected HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies
R [Ang] HF MP2 CCSD(T) [kJ/mol]
--------------------------------------------------------------
3.0 17.222357 11.528070 12.573556
3.5 2.849476 0.345935 0.710527
3.8 1.151019 -0.518605 -0.290598
4.0 0.464796 -0.662783 -0.514168
4.2 0.186994 -0.587881 -0.487227
4.5 0.074069 -0.468666 -0.398117
4.8 0.028024 -0.359223 -0.308316
5.0 0.009329 -0.271765 -0.234170
5.5 -0.000799 -0.155863 -0.134303
6.0 -0.001838 -0.091985 -0.079003Here are the aug-cc-pVTZ results:
<code | aug-cc-pVTZ>
CP-corrected HF/MP2/CCSD(T)/aug-cc-pVTZ interaction energies
R [Ang] HF MP2 CCSD(T) [kJ/mol]
--------------------------------------------------------------
3.0 17.188613 9.145848 0.000000
3.5 2.869509 -0.396962 0.000000
3.8 1.149645 -0.981602 0.000000
4.0 0.456303 -0.958149 0.000000
4.2 0.180387 -0.775824 0.000000
4.5 0.071272 -0.588018 0.000000
4.8 0.028000 -0.435874 0.000000
5.0 0.010724 -0.322174 0.000000
5.5 0.001081 -0.179785 0.000000
6.0 -0.000352 -0.104699 0.000000And finally the aug-cc-pVQZ results:
<code | aug-cc-pVQZ>
CP-corrected HF/MP2/CCSD(T)/aug-cc-pVQZ interaction energies
R [Ang] HF MP2 CCSD(T) [kJ/mol]
--------------------------------------------------------------
3.0 17.076285 7.893514 0.000000
3.5 2.873486 -0.746453 0.000000
3.8 1.160026 -1.165904 0.000000
4.0 0.463245 -1.062367 0.000000
4.2 0.182901 -0.840188 0.000000
4.5 0.071422 -0.630019 0.000000
4.8 0.027613 -0.463643 0.000000
5.0 0.010560 -0.340454 0.000000
5.5 0.001384 -0.187811 0.000000
6.0 0.000002 -0.108494 0.000000Important
Psi4 uses density-fitting (DF) for all these calculations. This results in an error - a small one, but it does seem to show up in the HF interaction energies in the smaller basis sets. HF does not bind this system, but the interaction energy is negative for the larger separations.
It is also possible that these negative energies result from the CP correction, which will be larger for the smaller basis sets.
Error in CCSD(T) part
Important
This is now fixed. We needed to use the clean() command between calculations. Thanks to Stephen and Charlie.
The CCSD(T) part did not work. If you uncomment the lines in the above file that refer to the CC part of the calculation and run it, it will fail. The first CC calculation would be OK, but the second would fail with the error:
<code | Error message>
*** tstart() called on comanche
*** at Sun Dec 18 20:12:51 2016
An error has occurred python-side. Traceback (most recent call last):
File "<string>", line 57, in <module>
File "/opt/ajm/Psi4/psi4conda/share/psi4/python/driver.py", line 446, in energy
wfn = procedures['energy'][lowername](lowername, molecule=molecule, **kwargs)
File "/opt/ajm/Psi4/psi4conda/share/psi4/python/procedures/proc.py", line 858, in select_ccsd_t_
return func(name, **kwargs)
File "/opt/ajm/Psi4/psi4conda/share/psi4/python/procedures/proc.py", line 1926, in run_ccenergy
psi4.cctransort(ref_wfn)
RuntimeError:
Fatal Error: PSIO Error
Error occurred in file: /scratch/cdsgroup/conda-builds/work/src/lib/libpsio/error.cc on line: 128
The most recent 5 function calls were:
psi::PsiException::PsiException(std::string, char const*, int)
psi::psio_error(unsigned int, unsigned int)
psi::PSIO::write(unsigned int, char const*, char*, unsigned long, psi::psio_address, psi::psio_address*)
psi::PSIO::write_entry(unsigned int, char const*, char*, unsigned long)
psi::cctransort::cctransort(boost::shared_ptr<psi::Wavefunction>, psi::Options&)I have not traced this error, but it seems to be related to one job restarting from another. As the jobs are of different sizes, this may not work. We may need some way of clearing spaces before the next CC calculation, or maybe we need to do the dimers together and then the monomers.
This assumption is correct: I can run multiple CC calculations on the dimer in succession. But it fails with the above error when attempting the monomer CC (the HF and MP2 work).
Supermolecular: DFT (various)
Here is a sample file for supermolecular DFT calculations of the interaction energy:
<code | Ar2_DFT_scan_CPexplicit.psi4>
#! Example potential energy surface scan and CP-correction for Ar2
#Based on a Psi4 example from the wiki
#http://psicode.org/psi4manual/master/tutorial.html#analysis-of-intermolecular-interactions
memory 10 gb
molecule dimer {
Ar
Ar 1 R
}
#Define a monomer with a ghost atom. We need this for the CP correction.
#Only one monomer need be defined as, for this system, the other is
#identical by symmetry.
molecule mono {
Ar
@Ar 1 R
}
Rvals=[3.0, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0]
#Rvals=[3.0, 3.5]
#molecules=[dimer,mono]
# Define the functional
# This did not seem to work!
#set dft_functional wB97X-D
#
# Instead use this to define the functional:
#dft_func = "PBE0-D3MBJ"
#dft_func = "wB97X-D"
#dft_func = "dlDF+D09"
dft_func = "dlDF+D"
basis_set = "aug-cc-pVDZ"
set basis $basis_set
set freeze_core True
# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_dft = {}
e_d_disp = {}
e_m_dft = {}
e_m_disp = {}
for R in Rvals:
dimer.R = R
energy(dft_func,molecule=dimer)
e_d_dft[R] = get_variable('DFT FUNCTIONAL TOTAL ENERGY')
e_d_disp[R] = get_variable('DISPERSION CORRECTION ENERGY')
#
mono.R = R
energy(dft_func,molecule=mono)
e_m_dft[R] = get_variable('DFT FUNCTIONAL TOTAL ENERGY')
e_m_disp[R] = get_variable('DISPERSION CORRECTION ENERGY')
psi4.print_out("\n")
psi4.print_out("CP-corrected %s / %s interaction energies\n\n" % (dft_func, basis_set))
psi4.print_out(" R [Ang] DFT DFT+D [kJ/mol] \n")
psi4.print_out("---------------------------------------------------------\n")
for R in Rvals:
e_dft = (e_d_dft[R] -2.0*(e_m_dft[R])) * psi_hartree2kJmol
e_dftD = (e_d_dft[R]+e_d_disp[R] -2.0*(e_m_dft[R]+e_m_disp[R])) * psi_hartree2kJmol
psi4.print_out(" %3.1f %10.6f %10.6f \n" % (R, e_dft, e_dftD))The table will contain both the pure DFT and DFT+D interaction energies.
Here are some example results:
The wB97XD functional (Psi4 code: wB97X-D) is a range-separated hybrid functional that is dispersion-corrected with the XD dispersion model from Becke and Johnson.
<code | wB97XD aug-cc-pVDZ >
CP-corrected DFT/aug-cc-pVDZ interaction energies
R [Ang] DFT DFT+D [kJ/mol]
---------------------------------------------------------
3.0 12.246337 11.779184
3.5 1.796835 0.952919
3.8 0.742007 -0.148517
4.0 0.299477 -0.506101
4.2 0.078264 -0.578098
4.5 -0.017808 -0.524060
4.8 -0.037940 -0.420009
5.0 -0.029290 -0.316492
5.5 -0.007502 -0.172609
6.0 -0.002706 -0.101213Here is another functional: The PBE0 functional is a hybrid functional. This one has been augmented with Grimme's D3 dispersion correction with the MBJ model (I'm not entirely sure what this is, but the MBJ likely refers to the switching function parameters that joins D3 to PBE0):
<code | PBE0D3MBJ aug-cc-pVDZ>
CP-corrected DFT/aug-cc-pVDZ interaction energies
R [Ang] DFT DFT+D [kJ/mol]
---------------------------------------------------------
3.0 10.986621 8.135669
3.5 0.782656 -0.965454
3.8 -0.081994 -1.377836
4.0 -0.289364 -1.233966
4.2 -0.280621 -0.965640
4.5 -0.213948 -0.711926
4.8 -0.146843 -0.511393
5.0 -0.095005 -0.364408
5.5 -0.036315 -0.188226
6.0 -0.014175 -0.103574Important
How large a basis do we need to converge the DFT energies?
Important
Compare the DFT results with those from MP2 and CCSD(T). Is there any pattern?
You will notice that some DFT models show some binding even without a dispersion correction. But this tends to be erratic.
The functionals often do not perform well for rare-gas dimers. You will have better luck with systems like the water dimer.
Warning
Psi4 is not perfect: some methods/functionals do not work as they should. You have to get into the habit of checking the output file to see if all seems well. Look for keywords, the energies computed, basis sets used,... Developing an ability to spot errors is an art. There is no recipe. You just have to do this for yourself.