Contents
Using mid-bonds and user-defined basis sets with Psi4
To use mid-bond functions we need to provide Psi4 with a user-defined basis. Instructions for this can be found in the Psi4 User-defined basis documentation. Here I will show you how this can be done for the mid-bond sets we use in most SAPT/SAPT(DFT) calculations. Of course, these sets can also be used in supermolecular calculations of the interaction energy, as the examples below will demonstrate.
User-defined basis sets
Create these in a file that takes the form: BASIS-NAME.gbs. There is a naming convention used for this file name. For details see the Psi4 documentation on basis sets.
Here is the file containing the three mid-bond sets we use:
<code | mb-set.gbs > spherical '''''' Kr 0 S 1 1.00 0.553063 1.00000000 S 1 1.00 0.250866 1.00000000 S 1 1.00 0.117111 1.00000000 P 1 1.00 0.392 1.00000000 P 1 1.00 0.142 1.00000000 D 1 1.00 0.328 1.00000000 '''''' Xe 0 S 1 1.00 0.393863 1.00000000 S 1 1.00 0.195206 1.00000000 S 1 1.00 0.800468 1.00000000 P 1 1.00 0.8025 1.00000000 P 1 1.00 0.3352 1.00000000 D 1 1.00 0.5 1.00000000 F 1 1.00 0.5 1.00000000 '''''' Rn 0 S 1 1.00 0.9 1.00000000 S 1 1.00 0.3 1.00000000 S 1 1.00 0.1 1.00000000 P 1 1.00 0.9 1.00000000 P 1 1.00 0.3 1.00000000 P 1 1.00 0.1 1.00000000 D 1 1.00 0.6 1.00000000 D 1 1.00 0.2 1.00000000 F 1 1.00 0.6 1.00000000 F 1 1.00 0.2 1.00000000 ''''''
Download this file and place it in your Psi4 basis library.
You will also need the following RI basis for the calculations using density-fitting (DF-MP2 etc):
<code | mb-set-ri.gbs> # RI basis for mb-set : AJ Misquitta : Podeszwa et al. J. Phys. Chem. A, Vol. 110, No. 34, 2006 spherical '''''' Kr 0 S 1 1.00 1.1061 1.0 S 1 1.00 0.5017 1.0 S 1 1.00 0.2342 1.0 P 1 1.00 0.94 1.0 P 1 1.00 0.5 1.0 P 1 1.00 0.25 1.0 D 1 1.00 0.90 1.0 D 1 1.00 0.60 1.0 D 1 1.00 0.30 1.0 F 1 1.00 0.7 1.0 F 1 1.00 0.4 1.0 G 1 1.00 0.65 1.0 '''''' Xe 0 S 1 1.00 1.8 1.0 S 1 1.00 1.2 1.0 S 1 1.00 0.6 1.0 S 1 1.00 0.4 1.0 S 1 1.00 0.2 1.0 P 1 1.00 1.8 1.0 P 1 1.00 1.2 1.0 P 1 1.00 0.6 1.0 P 1 1.00 0.4 1.0 P 1 1.00 0.2 1.0 D 1 1.00 1.8 1.0 D 1 1.00 1.2 1.0 D 1 1.00 0.6 1.0 D 1 1.00 0.4 1.0 D 1 1.00 0.2 1.0 F 1 1.00 1.5 1.0 F 1 1.00 0.9 1.0 F 1 1.00 0.5 1.0 F 1 1.00 0.3 1.0 G 1 1.00 1.5 1.0 G 1 1.00 0.9 1.0 G 1 1.00 0.3 1.0 '''''' Rn 0 S 1 1.00 1.8 1.0 S 1 1.00 1.2 1.0 S 1 1.00 0.6 1.0 S 1 1.00 0.4 1.0 S 1 1.00 0.2 1.0 P 1 1.00 1.8 1.0 P 1 1.00 1.2 1.0 P 1 1.00 0.6 1.0 P 1 1.00 0.4 1.0 P 1 1.00 0.2 1.0 D 1 1.00 1.8 1.0 D 1 1.00 1.2 1.0 D 1 1.00 0.6 1.0 D 1 1.00 0.4 1.0 D 1 1.00 0.2 1.0 F 1 1.00 1.5 1.0 F 1 1.00 0.9 1.0 F 1 1.00 0.5 1.0 F 1 1.00 0.3 1.0 G 1 1.00 1.5 1.0 G 1 1.00 0.9 1.0 G 1 1.00 0.3 1.0 ''''''
Important
Your Psi4 basis library will be in
- psi4conda/share/psi4/basis/
Once this is done the basis becomes accessible under the name mb-set as follows:
<code | example of usage>
basis {
assign aug-cc-pVDZ
assign Rn mb-set # gets the basis set for the ghost Rn atoms from mb-set.gbs
}In this example, all atoms get the default aug-cc-pVDZ basis and Rn gets mb-set.
Important
Why have I used the heavier rare-gas atoms for the mid-bond set?
This is because NWChem needs all sites, even the dummy ones, to have recognised atomic symbols. So it will not do to label a mid-bond site as Mb. As it is unlikely that we will use the heavier rare-gas atoms in a calculation, I have used these names. And, as you can see from above, the sizes of the mid-bond sets increases from Kr, through Xe, to Rn. Use Rn for the largest set.
In the conventional notation, Kr is Mb1, Xe is Mb2, and Rn is Mb3.
Ar$_2$ interaction energy with mid-bonds
Here's the Psi4 command file:
<code | Ar2_DC+_basis.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
@Rn 1 RMb
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
@Rn 1 RMb
@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, 3.75, 4.0 ]
molecules=[dimer,mono]
basis {
assign aug-cc-pVDZ
assign Rn mb-set # gets the basis set for the ghost Rn atoms from mb-set.gbs
}
# set scf_type df
# set cc_type df
set freeze_core True
# Initialize a blank dictionary of counterpoise corrected energies
# (Need this for the syntax below to work)
e_d_hf = {}
e_d_mp2 = {}
e_d_cc = {}
e_m_hf = {}
e_m_mp2 = {}
e_m_cc = {}
for R in Rvals:
dimer.R = R
dimer.RMb = R/2.0 # position of the midbond
#energy('mp2',molecule=dimer)
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
mono.RMb = R/2.0 # position of the midbond
#energy('mp2',molecule=mono)
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')
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
#e_cc = 0.0
psi4.print_out(" %3.1f %10.6f %10.6f %10.6f \n" % (R, e_hf,e_mp2,e_cc))This file should be easy to understand. The mid-bond set is placed (in this case) at the mid-point of the COMs of the two systems. So its location is easy to compute. This should be the result of this calculation:
Ar2 aug-cc-pVDZ DC+ (Mb3/Rn) results
<code | Ar2 aug-cc-pVDZ DC+ (Mb3/Rn) results>
CP-corrected HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies
R [Ang] HF MP2 CCSD(T) [kJ/mol]
--------------------------------------------------------------
3.0 17.147468 8.304114 8.842833
3.5 2.850430 -0.801159 -0.673696
3.8 1.144161 -1.226052 -1.155367
4.0 0.456287 -1.102277 -1.056947
'''* Psi4 exiting successfully. Buy a developer a beer!
Mid-bonds and DF-MP2, DF-SCF
These modules need the RI basis sets too. Without any RI basis specified I get the following error:
$ Traceback (most recent call last):
File "/home/alston/Psi4/psi4conda/share/psi4/python/qcdb/libmintsbasisset.py", line 599, in pyconstruct
bs, msg = BasisSet.construct(parser, mol, fitrole, None if fitrole == 'BASIS' else fitrole, basstrings[fitrole])
File "/home/alston/Psi4/psi4conda/share/psi4/python/qcdb/libmintsbasisset.py", line 748, in construct
(at + 1, role, text2))
qcdb.exceptions.BasisSetNotFound: BasisSet::construct: Unable to find a basis set for atom 2 for role JKFIT among:
Shell Entries: ['RN']
Basis Sets: ['def2-qzvpp-jkfit']
File Path: , /home/alston/Psi4/psi4conda/share/psi4/basis
Input Blocks: This is with commands:
basis {
assign aug-cc-pVDZ
assign Rn mb-set # gets the basis set for the ghost Rn atoms from mb-set.gbs
}
set scf_type df
set cc_type df
set freeze_core TrueThe Psi4 page http://www.psicode.org/psi4manual/master/basissets.html#sec-basisuserdefined says that we should define the DF-MP2 basis using
basis {
assign cc-pVDZ-RI df_basis_mp2
}
or to make an assignment for a specific atom:
basis {
assign C aug-cc-pVDZ-RI df_basis_mp2
}As I got a JKFIT error, I used df_basis_jkfit to specify the basis and got this error:
$ Traceback (most recent call last):
File "/home/alston/Psi4/psi4conda/share/psi4/python/inputparser.py", line 732, in process_input
temp = re.sub(basis_block, process_basis_block, temp)
File "/home/alston/Psi4/psi4conda/lib/python2.7/re.py", line 155, in sub
return _compile(pattern, flags).sub(repl, string, count)
File "/home/alston/Psi4/psi4conda/share/psi4/python/inputparser.py", line 356, in process_basis_block
raise TestComparisonError(message)
p4util.exceptions.TestComparisonError: Conflicting basis set specification: assign lines present but shells have no [basname] label.The correct syntax to select an auxiliary basis for all atoms is
assign df_basis_scf <NAME> assign df_basis_mp2 <NAME>
This is not what is stated on the Psi4 wiki.
Psi4 will work correctly if I use the following commands:
set {
Print 3 # will print information about the basis sets in the code
}
basis {
assign aug-cc-pVDZ
assign df_basis_scf aug-cc-pVDZ-JKFIT
assign Rn mb-set # gets the basis set for the ghost Rn atoms from mb-set.gbs
# assign Rn df_basis_scf mb-set-ri # this does not work
}
set scf_type df
set cc_type df
set freeze_core TrueIn addition to the above, I needed to inlcude the auxiliary basis for "Rn" in def2-qzvpp-jkfit.gbs. When this was done Psi4 worked correctly.
Interim solution
- Use only the 'Rn' mid-bond set.
- Edit file def2-qzvpp-jkfit.gbs. Right at the bottom include the Rn auxiliary-basis set from mb-det-ri.gbs (shown above). Only this set! Do the same for def2-qzvpp-ri.gbs
- Define the mid-bond site and basis sets as follows. Also include a specific definition of the DF basis sets:
molecule dimer {
Ar
@Rn 1 RMb
Ar 1 R
}
set {
Print 3
}
basis {
assign aug-cc-pVTZ # This is the main basis.
assign df_basis_scf aug-cc-pVTZ-JKFIT # This and the next DF basis set should be
assign df_basis_mp2 aug-cc-pVTZ-RI # consistent with the main basis
assign Rn mb-set # Basis set for the ghost Rn atoms from mb-set.gbs
}
set scf_type df
set mp2_type df
set cc_type df
Psi4 : Proper solution
Here's the way to solve the problem properly.
From https://github.com/psi4/psi4/issues/626:
basis {
assign aug-cc-pVDZ
assign Rn mb-set
}
df_basis_scf {
assign aug-cc-pVDZ-jkfit
assign Rn mb-set-ri
}
df_basis_cc {
assign aug-cc-pVDZ-RI
assign Rn mb-set-ri
}If you leave your mb-set.gbs and mb-set-ri.gbs alongside the input file, The parser will catch them. But you can also add directories to the envvar PSIPATH for it to look for basis set files. Please let me know if basis sets continue to give problems or confusion.