
<<TableOfContents()>>

= 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 [[http://www.psicode.org/psi4manual/master/basissets.html#sec-basisuserdefined|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 [[http://www.psicode.org/psi4manual/master/basissets.html#|the Psi4 documentation on basis sets.]]

Here is the file containing the three mid-bond sets we use:

[[attachment:mb-set.gbs]]

{{{
<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):

[[attachment:mb-set-ri.gbs]]

{{{
<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
''''''
}}}

{{{#!wiki important
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:

[[attachment:example of usage]]

{{{
<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''.


{{{#!wiki important
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:

[[attachment:Ar2_DC+_basis.psi4]]

{{{
<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:

[[attachment:Rn__results|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 True
}}}

The 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 True
}}}
In 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.