Introduction+to+Psi4

=Installing Psi4 (default) and setting up for solving problems with larger basis sets (angular momentum functions like H, I, J, or K)=

On Ruby... ./setup --max-am-eri 6 --cc icc --cxx icpc --fc ifort --mkl=parallel --prefix=`pwd` --type=release --pcmsolver=off --ambit=off --chemps2=off --gdma=off --cxx11=on objdir

DEFAULT
** make ** --- if building boost module this takes awhile (will also increase when we increase the angular momentum fxn limits) ** make tests ** --- these take quite a bit of time ** make install ** [if the tests passed obviously] ** make clean ** [optional for saving space, do ONLY in the obj directory]  Thanks to Brian at OSC, we have Psi4 on the supercomputer! I use the following configure command with an additional call to load the MKL module:
 * cd PATH/TO/PSI4 **
 * mkdir obj **
 * cd obj **
 * ../configure --prefix=/where/to/install **
 * ../configure --prefix=/nfs/04/ucn1093/quantum/psi4.0b5_mod --with-blas='-mkl' --with-cc=icc --with-cxx=icpc --with-fc=ifort --with-opt='-O2 -static -no-prec-div' --with-incdirs=-mkl CFLAGS="$MKL_CFLAGS" FFLAGS="$MKL_FFLAGS" LDFLAGS="$MKL_LIBS" --with-max-am-eri=8 --with-max-am-deriv1=6 --with-max-am-deriv2=6 **
 * make -j4 **
 * make tests **
 * make install **
 * make clean **

Add the executables and PSIDATADIR=/PATH/TO/PSI4/share/psi variables to your .bashrc and away you go. My compilation is modded to include source code and Boost package from the public master branch (essentially a pre-beta 6 release). A number of additional basis sets are included as well.  =** Errors and how to resolve them... ** =

vector3.i not found: (from obj directory) 

=Adjust frozen core options (for alkali Earths) =

Default scheme is to implement nearest Noble gas frozen core approximation which is good for anions, awful for mono- and divalent cations. Here's the tweak to the code to rename alkali Earth metals to use the (N-1)th Noble gas configuration for frozen core (up to K+). A similar approach can be used to adjust any atom. The number you define corresponds to...

File: Path/To/Psi4/src/lib/libmints/molecule.cc ==

== =Faux-JKfit basis sets for RI-SCF methods=

If you study cations in PSI4, you have a problem -- none of them are in the JKfit series for the Dunning type basis sets. This forces the use of the PK algorithm which I'm not sure works 100% correctly as it spends ages writing and sorting 4-center integrals. If, however, you supply a reasonable guess at a JKfit basis set you can take advantage of the fact that all other algorithms will attempt a density fit calculation before proceeding to solve the SCF equations exactly. A common rule is to uncontract the basis set of (zeta+1) quality to the parent basis set and call that one the JKfit one. Here's an example where I took the Na+ cc-pVQZ basis set, uncontracted it, and supplied it in a separate file for use as a 'JKfit' basis set. Since these basis sets are not optimized, your exact algorithm will have to work a bit harder, but it does knock off a handful of iterations and gets rid of the need to perform the awful PK method. Just use scf_type out_of_core.

When paired with the aug-cc-pVDZ parent basis set, we've found the 6-311++G(2d,2p) basis set when uncontracted to give far better energy convergence to begin the Direct or PK algorithms with (at least for sodium). You must use the TZ Pople basis sets here to ensure all functions are optimized for spherical coordinates (5D, 7F). Polarization functions in DZ Pople basis sets are Cartesian optimized (6D, 10F).



=Corrections to the 2nd and 3rd order exchange-induction coupled terms... =

By default, SAPT in Psi4 may have issues for strongly interacting species when modeling the repulsive wall for R < R_eq. John Herbert of OSU recently published on the weakened, unquenched exch-ind terms to 3rd order due to the use of the single-exchange approximation (SEA). To account for this he furthers an idea from a previous study (author eludes me for the moment), which scales these terms by the ratio of [exch(10)]/[exch(10)(S^2)] where the numerator is the exact 1st order exchange term, the denominator being the single-exchange approximation of this very same 1st order contribution. He looks at adding a power to this scaling term and attempts to recover the CCSD(T) potential energy curve for fluoride with water (a definite trouble case -- from personal experience). In effect, this leaves us with the scaling term in the form... {[exch(10)]/[exch(10)(S^2)]}^a where a = 1 in the study Herbert expands upon. Herbert uses a = 1, 2, 3 but cannot quite recover the exact shape of the PEC for R < R_eq while you get very good agreement at R_eq and beyond. Still these corrections are better than nothing and without even correcting every single term that depends upon the values of exch-ind(20) and exch-ind(30), I recover a ~2 kcal/mol destabilizing contribution that pushes my SAPT2+3 interaction energy within 0.5 kcal/mol of the CCSD(T) supermolecular interaction energy.

So... in the code, I just append the corrections in the print lines (so instead of printing say exch-ind(20), it now prints exch-ind(20) (SCALED) and I append the ratio to the variable being printed. I do keep the original code intact -- just copy/paste the line, correct it, and add // to denote a change from the original at the tail of the line

I no longer correct for this in the code, I just use AWK to scale my answer at the end. Either way, the scaling is used as a fudge factor only.

JEAN CLAUDE VAN D-AM'ing your Psi4 install to run when using basis sets with very high angular momentum fxns
 Deriv1 and Deriv2 are derivative modules, set them to the lowest value you need <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">r12 is used in correlated methods, but I think the documentation says this value isn't used in the BETA version <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">ERI must be set to at least N+2, if N is the integer set for the derivatives... this is for electron repulsion integrals <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;"> <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">You'll find these values multiplied by 2 in the libint and libderiv source directories in your obj directory (the README in these folders explains). <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;"> <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">At this point you just follow the same make, test, and install procedure as spelled out above. However, before carrying this out, we need to change two things... <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;"> <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">** cd .. ** (should be the psi4 home directory) <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;"> <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">The changes I made are surrounded by, namely I added an additional zero to the MAXNODE and MAX_NUM_TARGET_VRR_NODES allowing me to compile the source code required for such high angular momentum functions. Once altered, this compiled for me in full without error. All test cases achieving a 'PASSED' mark. code format="c" /*-- Builds a library of functions-calls for applying HRR and VRR --*/
 * cd PATH/TO/PSI4 **
 * mkdir obj **
 * cd obj **
 * ../configure --prefix=/where/to/install --with-max-am-eri=8 --with-max-am-deriv1=6 --with-max-am-deriv2=6 **
 * cd src/lib/libderiv **
 * vim emit_deriv12_managers.c **


 * 1) include <stdio.h>
 * 2) include <string.h>
 * 3) include <stdlib.h>
 * 4) include <libint/libint.h>
 * 5) include "mem_man.h"
 * 6) include "build_libderiv.h"
 * 7) define MAXNODE 10000(0)
 * 8) define MAX_NUM_TARGET_VRR_NODES 2000(0)
 * 9) define NONODE -1000000
 * 10) define NUMPARENTS 40

code <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;"> <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif;">Examples using Psi4
code format="python"
 * 1) ! Hey, this is a Psi4 comment!
 * 2) ! We'll run an energy optimization of water using BLYP-D

molecule water_dimer { 0 1 #! charge, spin multiplicity O x y z H x y z H x y z -- 0 1 #! repeat for second water O x y z H x y z H x y z units angstrom symmetry c1 no_reorient #! deadly important in SAPT } #! close molecule
 * 1) ! First we'll define the molecule and split the
 * 2) ! non-covalently bound waters with the -- tag

set globals { basis cc-pvdz dft_functional blyp guess sad scf_type df puream true geom_maxiter 100 dft_basis_tolerance 1.0E-10 ints_tolerance 1.0E-8 g_convergence gau_tight full_hess_every 10 }
 * 1) ! Psi4 beta release 5 supports only JK-fit gradients in optimization
 * 2) ! Reduced tolerances here to speed up calculation and we ask for
 * 3) ! coordinates of Gaussian's 'tight' keyword quality
 * 4) ! NOTE: You can compute the Hessian every N-steps; need a for loop
 * 5) ! if you'd like to compute the Hessian at SCF-level to use in MP2, etc.

optimize('blyp-d') water_dimer.update_geometry e_dimer=get_variable("CURRENT ENERGY") energy('sapt2+') e_sapt=get_variable("SAPT ENERGY")
 * 1) ! Optimize coordinates and update $molecule

w_gho=water_dimer.extract_subsets(1,2) #! E(A) gho_w=water_dimer.extract_subsets(2,1) #! E(B) self1=water_dimer.extract_subsets(1) #! E(a) self2=water_dimer.extract_subsets(2) #! E(b)
 * 1) ! Binding energy is...
 * 2) ! E(AB)-E(A)-E(B)+[E(am)-E(a)]+[E(bm)-E(b)]
 * 3) ! Capital letter A, B means AB basis set
 * 4) ! Lowercase letter a, b means monomer basis
 * 5) ! Lowercase m means monomer geometry; else dimer geometry

activate(w_gho) set globals { basis cc-pvdz dft_functional blyp guess sad scf_type df puream true dft_basis_tolerance 1.0E-10 ints_tolerance 1.0E-8 }

energy('blyp-d') e_w_gho=get_variable("CURRENT ENERGY")

energy('blyp-d') e_self1=get_variable("CURRENT ENERGY") optimize('blyp-d') e_self1_vac=get_variable("CURRENT ENERGY") #! E(am/bm)
 * 1) ! Repeat for gho_w, self1, and self2
 * 2) ! Optimize self1 and self2 at BLYP-D using the same procedure as above
 * 3) ! So the calculations you'd perform for self1 and self2 are...

delE_uncorr=(e_dimer-e_w_gho-e_gho_w)*627.5095 delE_corr=(e_dimer-e_w_gho-e_gho_w+e_self1_vac-e_self1+e_self2_vac-e_self2)*627.5095 delE_sapt_uncorr=e_sapt*627.5095 delE_sapt_corr=(e_sapt+e_self1_vac-e_self1+e_self2_vac-e_self2)*627.5095
 * 1) ! Now we assemble the binding energy (corrected and uncorrected for
 * 2) ! monomer relaxation effects) and convert to kcal/mol

table=Table(rows=["unit: kcal/mol"],cols=["delE(corr)"."delE(SAPT-corr)","delE(uncorr)","delE(SAPT-uncorr)"] table[1]=[delE_corr,delE_sapt_corr,delE_uncorr,delE_sapt_uncorr] print table code

<span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">

<span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;"> <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;"> <span style="background-color: #ffffff; color: #222222; display: block; font-family: arial,sans-serif;">Last Updated: TPP 22/10/2013