## page was renamed from AJMGrpOnly/teaching/electronic-structure/psi4-argon2-pes-hf-mp2 ## page was renamed from ajm/teaching/electronic-structure/psi4-argon2-pes-hf-mp2 <> = Psi4: Ar$_2$ interaction energies = In these series of examples we will explore how to best calculate the interaction energies of the Ar dimer system. We will use only MP2 (Moller-Plesset perturbation theory at second-order) and Hartree-Fock (HF) theory. The latter does not lead to binding of rare-gas complexes as these are bound entirely by the van der Waals interaction (dispersion energy) with is the result of electron correlation. MP2 is not a perfect theory but we use it here as it is computationally efficient enough for our purposes. We will use the Psi4 code here. This is a cutting-edge research-level electronic structure program that is fast, versatile, and which contains a large number of techniques. You will find extensive documentation on [[http://www.psicode.org/|the Psi4 website.]] The code is written in C++, F90 and Python. The input file can use a lot of Python commands - something we will make use of. == Setting up Psi4 == In addition to the [[ajm/teaching/electronic-structure/getting-started|NWChem set up instructions]], do the following. Include the lines in your '''.bashrc''' file: {{{ # Psi4 paths export PSI4_HOME=/opt/ajm/Psi4/current export PATH="${PSI4_HOME}/bin:$PATH" # Make sure that this command comes after you have defined the SCRATCH variable: export PSI_SCRATCH=${SCRATCH} }}} Then execute it using {{{ $ source .bashrc }}} If done correctly, you will be able to find the Psi4 executable using the '''which''' command: {{{ $ which psi4 /opt/ajm/Psi4/current/bin/psi4 }}} If you get //command not found// then something is wrong and ask your demonstrators for help. == Running Psi4 == Psi4 is a parallel code and you need to specify the number of processors it will use using the invocation: {{{ $ psi4 -n 1 job.psi4 job_basis.out & }}} where '''job.psi4''' is the Psi4 command file and it will typically use a certain '''basis'''. I tend to label my output files using both the job name and the basis like '''job_basis.out'''. This makes it clear to me which Psi4 command file I used, and with which basis. In this class you will always use only 1 processor. == Basis sets == The basis sets we will use are the so called Dunning correlation-consistent basis set. You can see this from the [[https://bse.pnl.gov/bse/portal| EMSL Basis Set Exchange]]. Go there and see if you can find the basis sets for Argon. If you want to learn more on basis sets see [[https://www.google.co.uk/url?sa=t&rct=j&q=&esrc=s&source=web&cd=3&cad=rja&uact=8&ved=0ahUKEwi3lJvJ1NLXAhXJJlAKHUV1B1kQFgg9MAI&url=http%3A%2F%2Fwww.theochem.unito.it%2Fcrystal_tuto%2Fmssc2013_cd%2Ftutorials%2Fbasis_set%2Fbasis_sets_2000.pdf&usg=AOvVaw1F9TlyeknmIF7RGtED543U| this introductory article]] by a colleague of mine, Dr Mike Towler. In the examples below you will need to modify the files to use these basis sets: * cc-pVDZ * cc-pVTZ * aug-cc-pVDZ * aug-cc-pVTZ These are all Dunning basis sets. The important parts are the 'D' and 'T': these give you an idea of the size of the basis. A double-zeta ('D') basis is smaller than a triple-zeta ('T') basis. Also, the 'aug' means that the basis has been augmented with basis functions that are more slowly decaying (the //diffuse// functions). These tend to be better for van der Waals systems. == No Counterpoise (CP) correction == First of all let's see how the interaction energy can be calculated for a pair of argon atoms at a fixed separation. By definition (see my notes on [[ajm/teaching/intermolecular-interactions|Intermolecular Interactions]]) the interaction energy is given by $$ E_{\rm int} = E[A..B] - E[A] - E[B]. $$ Here $E[A..B]$ is the energy of the complex A..B, and $E[A]$ and $E[B]$ are the energies of the isolated subsystems. Here we have two argon atoms separated by distance $R$, so this becomes $$ E_{\rm int}(R) = E[{\rm Ar..Ar}] - 2 E[{\rm Ar}] $$ Therefore all we need to do is calculate the energy of two argon atoms separated by $R$ and subtract twice the energy of one argon. This can be done using the following command file for Psi4. See the comments in the file for explanations of what's going on. [[attachment:Ar2_MP2_single_noCP.psi4]] {{{ # Ar2 energy single geometry : MP2 : no counterpoise correction memory 8 gb molecule dimer { Ar Ar 1 3.75 } #Define a monomer with a ghost atom. #Only one monomer need be defined as, for this system, the other is #identical by symmetry. molecule mono { Ar } molecules=[dimer,mono] set { Print 3 } # Use basis sets: # cc-pVDZ # cc-pVTZ # aug-cc-pVDZ # aug-cc-pVTZ # basis { assign aug-cc-pVTZ } # This option makes the calculations faster without degrading the # energies by a significant amount. It causes the inner electrons # to be 'frozen', i.e., not considered for the more computationally # demanding parts of the calculation. set freeze_core True # First perform the monomer energy calculation. energy('mp2',molecule=mono) # Psi4 has variables that store important energies. We extract the # HF (==SCF) and MP2 energies using these and store them in our own (Python) # variables: e_m_hf = get_variable('SCF TOTAL ENERGY') e_m_mp2 = get_variable('MP2 TOTAL ENERGY') # We need this statement to get Psi4 to clear its settings. # This is needed as it currently has various arrays set to the size # of the previous calculation, but when we perform the Ar2 calculation, # the system and basis size will change. This will cause the code to # fail with a very cryptic error message. If you want to see it, # comment out this line and run this code. clean() # Now the same for the dimer: energy('mp2',molecule=dimer) # Once again we store the HF and MP2 energies in our variables: e_d_hf = get_variable('SCF TOTAL ENERGY') e_d_mp2 = get_variable('MP2 TOTAL ENERGY') #Remove tmp files etc or else the mono jobs will fail. clean() # Psi4 conversion factors # These constants can be used to convert energy units: # psi_hartree2kJmol # psi_hartree2wavenumbers # # Now define the interaction energy as # Eint = E[Ar..Ar] - E[Ar] - E[Ar] # = E[Ar..Ar] - 2*E[Ar] # and convert units from Hartree (atomic units) to wavenumbers (cm^{-1}) # We do this as the interaction energies of rare gas dimers is usually # expressed in these units. psi4.print_out("\n") psi4.print_out("no-CP correction HF/MP2 interaction energies\n\n") psi4.print_out(" R [Bohr] HF MP2 [cm^{-1}] \n") psi4.print_out("-------------------------------------------------------\n") e_hf = (e_d_hf -2.0*e_m_hf ) * psi_hartree2wavenumbers e_mp2 = (e_d_mp2 -2.0*e_m_mp2) * psi_hartree2wavenumbers psi4.print_out(" %6.3f %12.6f %12.6f \n" % (R, e_hf,e_mp2)) }}} It looks like a lot, but examine it carefully and you will see that most of it makes sense. We have used the //aug-cc-pVTZ// basis here. Run this using {{{ $ psi4 -n 1 Ar2_MP2_single_noCP.psi4 Ar2_MP2_single_noCP_aTZ_R3.75.out & }}} Notice I have included the basis set as well as the separation $R=3.75$ Angstrom in the output file name. The output file will be quite large. Have a look at it and see if you can spot the SCF (same as Hartree-Fock) energies and MP2 energies. At the very end of the file you should see the following: {{{ no-CP correction HF/MP2 interaction energies R [Bohr] HF MP2 [cm^{-1}] ------------------------------------------------------ 8.3 89.551891 -110.712792 '''* Psi4 exiting successfully. Buy a developer a beer! }}} The interaction energy at the HF level of theory is $89.55 ~ cm^{-1}$ and at the MP2 level it is $-110.71 ~ cm^{-1}$. This is using the //aug-cc-pVTZ// basis. There are a few things to note: * With such an input file we get the interaction energy for a single geometry only. * Also, we have used a single basis. Try out this calculation with each of the four basis sets listed above. * There is an error in this calculation that originated from the fact that the basis sets we have used are finite in extent in space. This is called the '''basis set superposition error, or BSSE.''' For now ignore this error. {{{#!wiki important Important Before proceeding to the next step make a note of the interaction energies as a function of the basis sets listed above. Do you see any pattern? Does the energy appear to converge? Reference $E_{\rm int}$ at a separation of $R=3.75$ Angstrom is: -99.47 cm-1 (Ref: [[attachment:ar2_accurate_potential_patkowskimfs05_molp103.pdf|Patkowski, Murdachaew, Fou, and Szalewicz, Mol. Phys., 103, 10-20 (2005)]] Have a look at this paper as it contains a lot more data.) }}} === Energy scan: no CP correction : Dimer basis === In the previous example we calculated the interaction energy for a single separation $R$. This is not sufficient and we can gain far more information when we use multiple separations. In this example you will calculate the interaction energy without the counterpoise correction for a range of $R$ values. Here we use a standard basis for the dimer (Ar$_2$) calculation and a monomer basis for Ar. Compute the interaction energies without the CP correction for the set of basis sets listed above. [[attachment:Ar2_MP2_scan_noCP.psi4]] {{{ # Ar2 energy scan : MP2 : no counterpoise correction memory 8 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 } # This array contains the list of separations at which the calculations will be performed. # Use the longer array only when you are sure the code is working correctly. # Rvals=[3.0, 3.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0] Rvals=[3.0, 3.5, 3.75 ] molecules=[dimer,mono] set { Print 3 } # Use basis sets: # cc-pVDZ # cc-pVTZ # aug-cc-pVDZ # aug-cc-pVTZ # basis { assign aug-cc-pVTZ } # This option makes the calculations faster without degrading the # energies by a significant amount. It causes the inner electrons # to be 'frozen', i.e., not considered for the more computationally # demanding parts of the calculation. 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 = {} # First perform the monomer energy calculation. As we are not using # the CP correction, this needs to be done only once. # Also, we do not need to put the energies in arrays. # energy('mp2',molecule=mono) # energy('ccsd(t)',molecule=mono) e_m_hf = get_variable('SCF TOTAL ENERGY') e_m_mp2 = get_variable('MP2 TOTAL ENERGY') # e_m_cc = get_variable('CCSD(T) TOTAL ENERGY') # We need this statement to get Psi4 to clear its settings. # This is needed as it currently has various arrays set to the size # of the previous calculation, but when we perform the Ar2 calculation, # the system and basis size will change. This will cause the code to # fail with a very cryptic error message. If you want to see it, # comment out this line and run this code. clean() for R in Rvals: dimer.R = R 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() # Psi4 conversion factors # psi_hartree2kJmol # psi_hartree2wavenumbers psi4.print_out("\n") psi4.print_out("no-CP correction HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies\n\n") psi4.print_out(" R [Ang] HF MP2 CCSD(T) [cm^{-1}] \n") psi4.print_out("---------------------------------------------------------------\n") for R in Rvals: e_hf = (e_d_hf[R] -2.0*e_m_hf ) * psi_hartree2wavenumbers e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2) * psi_hartree2wavenumbers # e_cc = (e_d_cc[R] -2.0*e_m_cc ) * psi_hartree2wavenumbers e_cc = 0.0 psi4.print_out(" %6.3f %12.6f %12.6f %12.6f \n" % (R, e_hf,e_mp2,e_cc)) }}} Run it using {{{ $ psi4 -n 1 Ar2_MP2_scan_noCP.psi4 Ar2_MP2_scan_noCP_aTZ.out & }}} This may take a few minutes. The result summary will be at the end of the output and should look like this {{{ no-CP correction HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies R [Ang] HF MP2 CCSD(T) [cm^{-1}] --------------------------------------------------------------- 3.0 1421.719436 683.451464 0.000000 3.5 229.679934 -75.276006 0.000000 3.8 89.551891 -110.712792 0.000000 '''* Psi4 exiting successfully. Buy a developer a beer! }}} Ignore the last column. It would normally contant energies calculated using a very high level of theory (CCSD(T)) but we will not use this method as it is quite computationally intensive. {{{#!wiki important Important Edit the command file to activate the long form of the Rvals array and re-run the job. Use the four basis sets. You will now be able to plot four graphs as a function of $R$. If you wish, include more separation values. }}} === Dimer basis with mid-bonds === {{{#!wiki warning Warning Skip this! }}} Mid-bond functions are basis functions placed between the interacting systems. In our case, this would be mid-way between the two argon atoms. These basis functions are known to significantly accelerate the convergence of the interaction energies. I.e., with the mid-bonds present you can use smaller basis sets and yet achieve results close to those from large basis calculations obtained without the mid-bonds. This is great as the computational cost thereby decreases. However mid-bond functions pose a massive problem when the counterpoise correction is not applied. And the sole purpose of this part of the exercise is to show you how bad a problem this can be. [[attachment:Ar2_MP2_scan_noCP_Mb.psi4]] {{{ # Ar2 energy scan : MP2 : no counterpoise correction # using Mid-Bonds (MB) # # The mid-bond functions accelerate convergence with basis # size. But here, when the CP correction is not applied, # they wreck havoc with the interaction energies! # The situation improves somewhat when we move to # larger basis sets, so you should see that the aug-cc-pVDZ # basis results in nonsense for the interaction energy, but # aug-cc-pVTZ is slightly better and aug-cc-pVQZ still better. # But none of these result in interaction energies we can # use. # memory 8 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 } # Rvals=[3.0, 3.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0] Rvals=[3.0, 3.5, 3.75 ] molecules=[dimer,mono] set { Print 3 } basis { assign aug-cc-pVDZ assign Rn mb-set } 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 = {} # First perform the monomer energy calculation. As we are not using # the CP correction, this needs to be done only once. # Also, we do not need to put the energies in arrays. # energy('mp2',molecule=mono) # energy('ccsd(t)',molecule=mono) e_m_hf = get_variable('SCF TOTAL ENERGY') e_m_mp2 = get_variable('MP2 TOTAL ENERGY') # e_m_cc = get_variable('CCSD(T) TOTAL ENERGY') # We need this statement to get Psi4 to clear its settings. # This is needed as it currently has various arrays set to the size # of the previous calculation, but when we perform the Ar2 calculation, # the system and basis size will change. This will cause the code to # fail with a very cryptic error message. If you want to see it, # comment out this line and run this code. clean() 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() # Psi4 conversion factors # psi_hartree2kJmol # psi_hartree2wavenumbers psi4.print_out("\n") psi4.print_out("no-CP correction HF/MP2/CCSD(T)/aug-cc-pVDZ interaction energies \n\n") psi4.print_out(" R [Ang] HF MP2 CCSD(T) [cm^{-1}] \n") psi4.print_out("---------------------------------------------------------------\n") for R in Rvals: e_hf = (e_d_hf[R] -2.0*e_m_hf ) * psi_hartree2wavenumbers e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2) * psi_hartree2wavenumbers # e_cc = (e_d_cc[R] -2.0*e_m_cc ) * psi_hartree2wavenumbers e_cc = 0.0 psi4.print_out(" %6.3f %12.6f %12.6f %12.6f \n" % (R, e_hf,e_mp2,e_cc)) }}} == Counterpoise Corrected interaction energies == To correct for the errors the result with finite basis sets we use what is known as the Boys and Bernadi //counterpoise correction//. The idea is very simple: calculate the energies of the dimer and monomers in the same basis set. That way we avoid, or suppress, any tendency for one of these energies to be inconsistent with the others. So here we do the following: * For separation //R//, * compute the dimer energy as before: Ar...Ar * But compute the monomer energy with basis functions located both on the Argon and //at the position where its partner would be//: Ar...Ar-Ghost. I have indicate the partner with //Ar-Ghost// as it is, in a way, a ghost atom. * The complication is that now the monomer energy depends on //R// as the Ar-Ghost position varies with //R//. Here is a command file to perform a CP-corrected calculation of the interaction energy for a fixed geometry. The only change from the uncorrected file (first one on this page) is that the monomer is defined with a ghost atom denoted as //@Ar//: [[attachment:Ar2_MP2_single_CP.psi4]] {{{ # Ar2 energy single geometry : MP2 : no counterpoise correction memory 8 gb R=3.75 molecule dimer { Ar Ar 1 3.75 } #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 3.75 } molecules=[dimer,mono] set { Print 3 } # Use basis sets: # cc-pVDZ # cc-pVTZ # aug-cc-pVDZ # aug-cc-pVTZ # basis { assign aug-cc-pVTZ } # This option makes the calculations faster without degrading the # energies by a significant amount. It causes the inner electrons # to be 'frozen', i.e., not considered for the more computationally # demanding parts of the calculation. set freeze_core True # First perform the monomer energy calculation. energy('mp2',molecule=mono) # Psi4 has variables that store important energies. We extract the # HF (==SCF) and MP2 energies using these and store them in our own (Python) # variables: e_m_hf = get_variable('SCF TOTAL ENERGY') e_m_mp2 = get_variable('MP2 TOTAL ENERGY') # We need this statement to get Psi4 to clear its settings. # This is needed as it currently has various arrays set to the size # of the previous calculation, but when we perform the Ar2 calculation, # the system and basis size will change. This will cause the code to # fail with a very cryptic error message. If you want to see it, # comment out this line and run this code. clean() # Now the same for the dimer: energy('mp2',molecule=dimer) # Once again we store the HF and MP2 energies in our variables: e_d_hf = get_variable('SCF TOTAL ENERGY') e_d_mp2 = get_variable('MP2 TOTAL ENERGY') #Remove tmp files etc or else the mono jobs will fail. clean() # Psi4 conversion factors # These constants can be used to convert energy units: # psi_hartree2kJmol # psi_hartree2wavenumbers # # Now define the interaction energy as # Eint = E[Ar..Ar] - E[Ar] - E[Ar] # = E[Ar..Ar] - 2*E[Ar] # and convert units from Hartree (atomic units) to wavenumbers (cm^{-1}) # We do this as the interaction energies of rare gas dimers is usually # expressed in these units. psi4.print_out("\n") psi4.print_out("CP corrected HF/MP2 interaction energies\n\n") psi4.print_out(" R [Bohr] HF MP2 [cm^{-1}] \n") psi4.print_out("-------------------------------------------------------\n") e_hf = (e_d_hf -2.0*e_m_hf ) * psi_hartree2wavenumbers e_mp2 = (e_d_mp2 -2.0*e_m_mp2) * psi_hartree2wavenumbers psi4.print_out(" %6.3f %12.6f %12.6f \n" % (R, e_hf,e_mp2)) }}} As before, run this using {{{ $ psi4 -n 1 Ar2_MP2_single_CP.psi4 Ar2_MP2_single_CP_aTZ_R3.75.out & }}} The output will contain {{{ CP corrected HF/MP2 interaction energies R [Bohr] HF MP2 [cm^{-1}] ------------------------------------------------------- 8.3 96.102802 -82.055480 }}} {{{#!wiki important Important Calculate the CP-corrected energies with the four basis sets and compare them with the uncorrected energies. What do you see? }}} === Energy scan with the CP correction : dimer basis === Here is a Psi4 command file that is very similar to the ones above, but this one performs an energy scan with the counterpoise correction applied. [[attachment:Ar2_MP2_scan_CP.psi4]] {{{ # Ar2 : MP2 energy scan : Counterpoise correction applied # Using mid-bonds (MB) and density-fitting (DF) #! 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.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0] Rvals=[3.0, 3.5, 3.75 ] molecules=[dimer,mono] # This is one way of defining a basis set: # In the other examples we use a different, and more flexible way. set basis aug-cc-pVTZ 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 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 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 conversion factors # psi_hartree2kJmol # psi_hartree2wavenumbers 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) [cm^{-1}] \n") psi4.print_out("---------------------------------------------------------------\n") for R in Rvals: e_hf = (e_d_hf[R] -2.0*e_m_hf[R] ) * psi_hartree2wavenumbers e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2[R]) * psi_hartree2wavenumbers # e_cc = (e_d_cc[R] -2.0*e_m_cc[R] ) * psi_hartree2wavenumbers e_cc = 0.0 psi4.print_out(" %6.3f %12.6f %12.6f %12.6f \n" % (R, e_hf,e_mp2,e_cc)) }}} === Energy scan with CP correction : dimer basis with mid-bonds === The above Psi4 command file works fine, but you may have noticed that there was no clear indication that the interaction energies had converged with basis set. By converged we mean that a large basis would not change the energies by much (or at all). Convergence can be vastly improved if we include the so-called //mid-bond basis sets//. These are basis sets placed in the middle of the interacting species. So in our case we'd have the following {{{ Ar....Mb....Ar }}} where '''Mb''' indicates the mid-bond set. In the following example we use mid-bond functions. But we also use a numerical device called //density-fitting// to significantly speed up our calculations. This can allow us to use even larger basis sets. The complication is that we now have to define two kinds of basis sets: {{{ basis { assign aug-cc-pVDZ # Do not alter this one: assign Rn mb-set } df_basis_scf { assign aug-cc-pVDZ-jkfit # Do not alter this one: assign Rn mb-set-ri } }}} Here the first basis (indicated simply as //basis//) is the main basis set. For technical reasons the '''Rn''' is not radon but is actually used for the mid-bond basis. The second basis is //df_basis_scf// and is the density-fitting basis. This needs to be matched with the //basis// definition. So if you alter the main basis to be '''aug-cc-pVTZ''' then also change the //df_basis_scf// to use '''aug-cc-pVTZ-jkfit'''. Here is the full command file: [[attachment:Ar2_MP2_scan_CP_DF_Mb.psi4]] {{{ # Ar2 : CCSD(T) energy scan : Counterpoise correction applied # Using mid-bonds (MB) and density-fitting (DF) # # Mid-Bonds # --------- # These are a small set of basis functions placed between the # interacting atoms. Here they are at the site labeled 'Rn'. # # Mid-bonds significantly accelerate the convergence of energies # with basis set. So you should find interaction energies calculated # using the mid-bond sets to be closer to the reference energies. # # Density-Fitting # --------------- # This is a numerical device used to significantly lessen the # computational cost of a calculation. It does not usually incurr # a significant loss in accuracy. memory 8 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.25, 3.5, 3.65, 3.75, 3.85, 4.0, 4.25, 4.5, 4.75, 5.0, 5.5, 6.0] Rvals=[3.0, 3.5, 3.75 ] molecules=[dimer,mono] set { Print 3 } # We have more basis sets defined here as we have mid-bonds and # are using density-fitting. # Alter the basis sets consistently. So if you change the first # basis to aug-cc-pVTZ, then also change the df_basis_set to # aug-cc-pVTZ-jkfit basis { assign aug-cc-pVDZ # Do not alter this one: assign Rn mb-set } df_basis_scf { assign aug-cc-pVDZ-jkfit # Do not alter this one: assign Rn mb-set-ri } #df_basis_cc { # assign aug-cc-pVDZ-RI # # Do not alter this one: # assign Rn mb-set-ri #} set scf_type df set mp2_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 conversion factors # psi_hartree2kJmol # psi_hartree2wavenumbers 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) [cm^{-1}] \n") psi4.print_out("---------------------------------------------------------------\n") for R in Rvals: e_hf = (e_d_hf[R] -2.0*e_m_hf[R] ) * psi_hartree2wavenumbers e_mp2 = (e_d_mp2[R] -2.0*e_m_mp2[R]) * psi_hartree2wavenumbers # e_cc = (e_d_cc[R] -2.0*e_m_cc[R] ) * psi_hartree2wavenumbers e_cc = 0.0 psi4.print_out(" %6.3f %12.6f %12.6f %12.6f \n" % (R, e_hf,e_mp2,e_cc)) }}} {{{#!wiki important Important Run this with the four basis sets. Make sure to change both //basis// and //df_basis_scf//. These calculations must be compared with the results in the paper by [[attachment:ar2_accurate_potential_patkowskimfs05_molp103.pdf|Patkowski, Murdachaew, Fou, and Szalewicz, Mol. Phys., 103, 10-20 (2005)]] Have a look at this paper as it contains a lot more data.) Graph the data for all basis sets. Do you now see a convergence with basis? How well do your results agree with those from Patkowski et al.? }}} == Final results of the practical == At the end of this lab session you should have consolidated your data into a single table saved as a text file. The table may take the form: [[attachment:Ar2_MP2.txt]] {{{ # Ar2 MP2 interaction energies # Units: R : Angstrom # E : cm^{-1} # # --------Using mid-bonds--- ----No mid-bonds --- # R aDZ-noCP aTZ-noCP DZ-noCP TZ-noCP aDZ-CP-Mb aTZ-CP-Mb DZ-CP-Mb TZ-CP-Mb aDZ-CP aTZ-CP ... 3.000 849.203 683.451 ... ... 693.597 646.098 ... ... ... 3.250 205.203 106.749 114.204 94.645 3.500 -14.944 -75.276 -66.813 -76.199 3.650 -63.415 -105.688 -97.067 -104.068 ... ... }}} This table will contain at least 12 data sets: * no-CP: 4 basis sets * CP: 4 basis sets * CP with mid-bonds : 4 data sets Use column headers that are short and descriptive. Only the MP2 data are important as HF yields no binding. {{{#!wiki important Important Also include the reference data from the paper supplied above. And make sure you include this data in your plots. }}} You should visualise the data by plotting it. I suggest you choose $x$- and $y$-axis limits to emphasise the bowl of the potential well. I used all $R$ values, but restricted the energy range to $-150$ to $100$ cm$^{-1}$. I suggest you use '''Gnuplot''' to graph the data. This code is available for all platforms so you will be able to install it on your laptop. It is not present on our workstations. {{{#!wiki important Important You will find that you cannot include all the data sets on one graph as you will not be able to make sense of it. One possibility is that you make three sets of graphs, each containing: * the reference data (in bold, black) * the four basis sets for one method (noCP, CP, CP+Mid-bonds) This way you will be able to see what's going on very clearly. Use either the last or second-last column of energies from Table 1 of the Patkowski paper as your reference. These are both very high levels of theory. Here is a sample graph of part of the no-CP data: [[attachment:ar2_mp2_nocp.pdf]] }}} {{{#!wiki important Important What you will find is that your best MP2 results (CP+Mid-bonds in the aTZ basis) will result in an over-binding compared with the Patkowski et al. results. This is due to a short-coming of the MP2 method. Now suppose you used this best MP2 data to create a potential for the argon gas or liquid. What effect would the overbinding have on the thermodynamic properties of the gas/liquid? The only way to fix this problem is to use the CCSD(T) method. But this has a computational scaling of $N^7$, that is, if you double the basis set or system size, it will take $128$ times longer to run. So you see that we cannot easily improve the MP2 result. Well, we can for this system, but the CCSD(T) method cannot be applied to larger systems which not only consist of many more electrons, but for which we typically need to evaluate the interaction energy on many thousands of dimer configurations. We need an alternative... }}}