810 likes | 1.24k Views
MOLPRO. Adem Tekin. Yrd . Doc. Dr. A quantum chemistry package. 15.06.2012 Istanbul. Outline. Introduction Input structure and introductory examples Brief introduction to SCS-MP2 and DFT-SAPT Applications of SAPT: Acetylene-Benzene dimer Towards the Ac n Bz n aggregates
E N D
MOLPRO AdemTekin Yrd. Doc. Dr. A quantum chemistry package 15.06.2012 Istanbul
Outline • Introduction • Input structure and introductory examples • Brief introduction to SCS-MP2 and DFT-SAPT • Applications of SAPT: Acetylene-Benzene dimer • Towards the AcnBzn aggregates • Exercise: Calculation of interaction energies for water dimer at MP2, SCS-MP2, B3LYP-D, DFT-SAPT (PBE0AC), DFT-SAPT (LPBE0AC) and CCSD(T) levels Turbomole 15.06.2012
Molpro • Developed by H. –J. Werner (University of Stuttgart) and P. J. Knowles (Cardiff University) • for more information please visit to: • http://www.molpro.net • A user mailing list is available under http://www.molpro.net/mailman/listinfo/molpro-user/?portal= user&choice=User+mailing+list • In molpro, the emphasis is given to highly accurate computations, with extensive treatment of the electron correlation problem through the multiconfiguration-reference CI and coupled cluster methods. Turbomole 15.06.2012
Running molpro • Molpro is accessed using the molpro command. The syntax is: molpro[options] [datafile] • datafile is the input. Prepare your input with .com extension. • output will be written to a datafile.out file • -d dir1:dir2:… where dir1:dir2:… is a list of directories which may be used for creating scratch files. • There are more options… • Molpro can be run in parallel: -n specifies the number of cpu -N specifies the number of tasks on the node Turbomole 15.06.2012
Submitting molpro jobs on UHEM #!/bin/bash #BSUB -m karadeniz_temp1 #BSUB -a intelmpi # bukismidegistirmeyin !!! #BSUB -J triazine_avtz16 #isinizitanimlayacakbirisim #BSUB -o triazine_avtz16.out # bukismidegistirmeyin !!! #BSUB -e triazine_avtz16.err # bukismidegistirmeyin !!! #BSUB -q workshop#kuyrukismi #BSUB -n 16 # kullanilacakolanislemcisayisi #BSUB –P workshop #BSUB -R "span[ptile=8]" # bukismidegistirmeyin !!! echo 'Starting time:' date molpro-n 16-N 8 --no-xml-output $PWD/triazine_avtz16.com>triazine_avtz16.log echo 'Ending time:' date Turbomole 15.06.2012
Molpro input structure ***,title !title (optional) memory,4,m !memory specification (optional) file,1,name.int !permanent named integral file (optional) file,2,name.wfu !permanent named wavefunction file (optional) gprint,options !global print options (optional) gthresh,options !global thresholds (optional) gdirect[,options] !global direct (optional) gexpec,opnames !global definition of one-electron operators (optional) basis=basisname !basis specification. If not present, cc-pVDZ !is default geometry={...} !geometry specification var1=value,var2=value,... !setting variables for geometry and/or wavefunction {command,options!program or procedure name directive,data,option !directives for command (optional) ... } !end of command block --- !end of input (optional) Turbomole 15.06.2012
Molpro input examples (http://www.be.itu.edu.tr/~adem.tekin/molpro) ! SCF calcultion for water. (ex1.com) ***,h2o !A title r=1.85,theta=104!set geometry parameters geometry={O; !z-matrix geometry input H1,O,r; H2,O,r,H1,theta} hf !closed-shell scf TURBOMOLE usage philosophy The default basis is cc-pVDZ abbreviated as VDZ ! SCF calculation for water using vtz (ex2.com) ***,h2o cc-pVTZ basis !A title r=1.85,theta=104 !set geometry parameters geometry={O; !z-matrix geometry input H1,O,r; H2,O,r,H1,theta} basis=VTZ !use VTZ basis hf !closed-shell scf
Molpro input examples (http://www.be.itu.edu.tr/~adem.tekin/molpro) Turbomole :%s/water-ccsdt-avdz/ex1/g 15.06.2012
Molpro examples ! Geometry optimization at the HF level for water (ex3.com) ***,h2o !A title r=1.85,theta=104 !set geometry parameters geometry={O; !z-matrix geometry input H1,O,r; H2,O,r,H1,theta} basis=6-31g** !use Pople basis set hf !closed-shell scf optg !do scf geometry optimization Turbomole ! Single point CCSD(T) for water (ex4.com) ***,h2o !A title r=1.85,theta=104 !set geometry parameters geometry={O; !z-matrix geometry input H1,O,r; H2,O,r,H1,theta} basis=VTZ !use VTZ basis hf !closed-shell scf ccsd(t) !do ccsd(t) calculation 15.06.2012
Electron correlation methods • Hartree-Fock is a single determinant wave function method where molecular orbitals are varied. • Configuration Interaction (CI) is a multiple determinant wavefunction method where molecular orbitals are not varied. • CI wavefunction is constructed by starting with the HF wavefunction and making new determinants by promoting electrons from the occupied to unoccupied orbitals. • Depending on the number of excitations used to make the determinants, CI methods are called CIS, CISD, CISDT, CISDTQ and full CI. • Multi-configurational self-consistent field (MCSCF) calculations also use multiple determinants. In contrast to the CI, orbitals are also optimized. • An MCSCF calculation in which all combinations of the active space orbitals are included is called a complete active space self-consistent field (CASSCF) calculation. Turbomole 15.06.2012
Electron correlation methods • In a CASSCF wavefunction the occupied orbital space is divided into a set of inactive or closed-shell orbitals and a set of active orbitals.All inactive orbitals are doubly occupied in each Slater determinant. On the other hand, the active orbitals have varying occupations, and all possible Slater determinants (or CSFs) are taken into account which can be generated by distributing the Nact=Nel-2mclosed electrons in all possible ways among the active orbitals, where mclosed is the number of closed-shell (inactive) orbitals, and Nel is the total number of electrons. • It is possible to construct a CI wavefunction starting with an MCSCF calculation rather than starting with a HF wavefunction. This starting wavefunction is called the reference state. These calculations are called multi-reference configuration interaction (MRCI). • Coupled cluster (CC) calculations are similar to CI calculations in that the wavefunction is a linear combination of many determinants. The means for choosing the determinants in CC is more complex than CI. There are variants of CC: CCSD, CCSDT and CCSD(T). Turbomole 15.06.2012
Molpro examples ! Geometry optimization at the HF level for water (ex5.com) ***,h2o !A title r=1.85,theta=104 !set geometry parameters geometry={o; !z-matrix geometry input h1,O,r; h2,O,r,H1,theta} basis=vtz !use VTZ basis hf !closed-shell scf ccsd(t) !do ccsd(t) calculation casscf !do casscf calculation mrci !do mrci calculation Turbomole 15.06.2012
Molpro example – tables ! Put the results into a table (ex6.com) ***,h2o !A title r=1.85,theta=104 !set geometry parameters geometry={o; !z-matrix geometry input h1,O,r; h2,O,r,H1,theta} basis=vtz !use VTZ basis hf !closed-shell scf e(1)=energy !save scf energy in variable e(1) method(1)=program !save the string ’HF’ in variable method(1) ccsd(t) !do ccsd(t) calculation e(2)=energy !save ccsd(t) energy in variable e(2) method(2)=program !save the string ’CCSD(T)’ in variable method(2) casscf !do casscf calculation e(3)=energy !save scf energy in variable e(3) method(3)=program !save the string ’CASSCF’ in variable method(3) mrci !do mrci calculation e(4)=energy !save scf energy in variable e(4) method(4)=program !save the string ’MRCI’ in variable method(4) table,method,e !print a table with results title,Results for H2O, basis=$basis !title for the table
Molpro example – Procedures (ex7.com) ***,h2o !A title proc save_e !define procedure save_e if(#i.eq.0) i=0 !initialize variable i if it does not exist i=i+1 !increment i e(i)=energy !save scf energy in variable e(i) method(i)=program !save the present method in variable method(i) endproc !end of procedure r=1.85,theta=104 !set geometry parameters geometry={o; !z-matrix geometry input h1,O,r; h2,O,r,H1,theta} basis=vtz !use VTZ basis hf !closed-shell scf save_e !call procedure, save results ccsd(t) !do ccsd(t) calculation save_e !call procedure, save results casscf !do casscf calculation save_e !call procedure, save results mrci !do mrci calculation save_e !call procedure, save results table,method,e !print a table with results title,Results for H2O, basis=$basis !title for the table
Molpro example – do loops (ex8.com) ***,H2O potential symmetry,x !use cs symmetry geometry={ o; !z-matrix h1,o,r1(i); h2,o,r2(i),h1,theta(i) } basis=vdz !define basis set angles=[100,104,110] !list of angles distances=[1.6,1.7,1.8,1.9,2.0] !list of distances i=0 !initialize a counter do ith=1,#angles !loop over all angles H1-O-H2 do ir1=1,#distances !loop over distances for O-H1 do ir2=1,ir1 !loop over O-H2 distances(r1.ge.r2) i=i+1 !increment counter r1(i)=distances(ir1) !save r1 for this geometry r2(i)=distances(ir2) !save r2 for this geometry theta(i)=angles(ith) !save theta for this geometry hf; !do SCF calculation escf(i)=energy !save scf energy for this geometry ccsd(t); !do CCSD(T) calculation eccsd(i)=energc!save CCSD energy eccsdt(i)=energy !save CCSD(T) energy enddo !end of do loop ith enddo !end of do loop ir1 enddo !end of do loop ir2 {table,r1,r2,theta,escf,eccsd,eccsdt !produce a table with results head, r1,r2,theta,scf,ccsd,ccsd(t) !modify column headers for table save,h2o.tab !save the table in file h2o.tab title,Results for H2O, basis $basis !title for table sort,3,1,2} !sort table
Program control • ***,text ! Starting a job • --- ! Ending a job • MEMORY,n,scale; ! Allocating dynamic memory • MEMORY,90000 ! allocates 90 000 words of memory • MEMORY,500,K ! allocates 500 000 words of memory • MEMORY,2,M !allocates 2 000 000 words of memory • 1 word = 8 bytes • IF($method.eq.’HF’) then ... ENDIF • GTHRESH,key1=value1,key2=value2,. . . ! Global thresholds • ENERGY ! Convergence threshold for energy (default 1.d-6) • ORBITAL ! Convergence threshold for orbital optimization in the SCF program (default 1.d-5). • gthresh,energy=1d-8,orbital=1d-7 Turbomole 15.06.2012
Molpro example – Expectation values – dipole and quadrupole moments (ex9.com) ***,h2o properties geometry={o;h1,o,r;h2,o,r,h1,theta} !Z-matrix geometry input r=1 ang!bond length theta=104 !bond angle gexpec,dm,sm,qm compute dipole and quarupole moments $methods=[hf,multi,ci] !do hf, casscf, mrci do i=1,#methods !loop over methods $methods(i) !run energy calculation e(i)=energy dip(i)=dmz !save dipole moment in variable dip quadxx(i)=qmxx!save quadrupolemomemts quadyy(i)=qmyy quadzz(i)=qmzz smxx(i)=xx !save second momemts smyy(i)=yy smzz(i)=zz enddo table,methods,dip,smxx,smyy,smzz !print table of first and second moments table,methods,e,quadxx,quadyy,quadzz!print table of quadrupole moments
Geometry specification – using xyz coordinates geomtyp=xyz geometry={ 3 ! number of atoms This is an example of geometry input for water with an XYZ file O ,0.0000000000,0.0000000000,-0.1302052882 H ,1.4891244004,0.0000000000, 1.0332262019 H,-1.4891244004,0.0000000000, 1.0332262019 } hf Turbomole 15.06.2012
Basis set • BASIS,DEFAULT=VTZ,O=AVTZ,H=VDZ • BASIS=VTZ,O=AVTZ,H=VDZ • BASIS DEFAULT=VTZ O=AVTZ H=VDZ END Turbomole 15.06.2012
Density fitting • Density fitting can be used to approximate the integrals in spin restricted Hartree-Fock (HF), density functional theory(KS), second-order Møller-Plesset perturbation theory (MP2 and RMP2), explicitly correlated MP2 (MP2-F12), and all levels of closed-shell local correlation methods (LCC2, LMP2-LMP4, LQCISD(T), LCCSD(T)). • Density fitting is invoked by adding the prefix DF- to the command name, e.g. DF-HF, DF-KS, DF-MP2 and so on. • Symmetry is not implemented for density fitting programs. • By default, a fitting basis set will be chosen automatically that corresponds to the current orbital basis set and is appropriate for the method. For instance, if the orbital basis set is VTZ, the default fitting basis is VTZ/JKFIT for DF-HF or DF-KS, and VTZ/MP2FIT for DF-MP2. • Examples: BASIS=VTZ !use VTZ orbital basis DF-HF,DF_BASIS=VQZ !use VQZ/JKFIT fitting basis DF-MP2,DF_BASIS=VQZ !use VQZ/MP2FIT fitting basis Turbomole 15.06.2012
Calculating the intermolecular interactions in water dimer • MP2 • SCS-MP2 • B3LYP-D • DFT-SAPT • CCSD(T) Turbomole 15.06.2012
Spin component scaled (SCS) MP2 • SCS-MP2 is based on a separate scaling of the second-order parallel (αα+ββ) and antiparallel-spin (αβ) pair correlation energies e(2). • The SCS-MP2 approximation to the correlation energy Ecorr is given by • Here, pTand pS are empirical parameters determined to be 1/3 and 6/5. • If you perform an MP2 calculation, you already have the SCS-MP2 energy. Turbomole S. Grimme, J. Chem. Phys., 2003, 118, 9095.
The iterative formulation of Rayleigh-Schrödinger perturbation theory Turbomole
, , • Intermolecular perturbation theory Turbomole
Symmetry adapted perturbation theory • The total interaction energy : Turbomole • Contributions of third and higher orders : 15.06.2012
A more simple way to calculate intermolecular interaction energy • Supermolecular approach (SMA) Turbomole • Basis set superposition error (BSSE) • Employ Counterpoise (CP) technique to avoid BSSE in SMA • SAPT does not contain BSSE 15.06.2012
The analysis of different types of intermolecular interactions Turbomole 15.06.2012
T S SS • Applications of SAPT: Acetylene - Benzene Turbomole -2.9 -5.8 -10.8 15.06.2012
Comparison of SAPT with supermolecular calculations: Acetylene - Benzene Turbomole Aug-cc-pVTZ
Comparison of SAPT with supermolecular calculations: Acetylene - Furan Turbomole 15.06.2012 Aug-cc-pVTZ
The intermolecular PES: extrapolation to the cbs limit for Acetylene - Benzene Turbomole • 42 test geometries (including global minimum region) • Only • extrapolated • Extrapolation with double- and triple- zeta basis sets sufficient here formula of Halkier et al. [CPL 286 (1998) 243] 15.06.2012
Intermolecular potential energy surface of the Acetylene-Benzene dimer: Choice of geometries C. A. Hunter et al. J. Am. Chem. Soc. 112 (1990), 5525 J. Gauss et al. J. Phys. Chem. 104 (2000),2865 J. L. M. Martin et al. J. Chem. Phys. 108 (1998), 676
Asymptotic interactions: Electrostatic energy via multipole expansion and Dispersion energy Turbomole Quadrupole-Quadrupole interaction Quadrupole-Octapole interaction Octapole-Octapole interaction • Quadrupole moments and Dispersion coefficients are calculated only once. • Rotated for other orientations using Euler`s rotation matrices. • 10.0 ≤ R ≤ 49.0 Å, R increased by 3 Å • 10192 orientations produced in this manner. • 10192+ 693= 10885 (Total number of points used in fit procedure) 15.06.2012
Site-site fit model: Dispersion & Induction Repulsion Electrostatics • Fitting the intermolecular PES: Fit of extrapolated total interaction energy • Damping function: • Merit function (Chi-Square): C. Leforestier, A. Tekin, G. Jansen , M. Herman, J. Chem. Phys., 135, 234306, (2011)
For 693 points SD [kJ/mol]: 3.63 0.18 6.9 (15.5 %) For 10885 points SD [kJ/mol]: 0.92 -唴 唴 0.04 6.9 (15.5 %) Eint [kJ/mol] 2.6 -唴 唴 Eint [kJ/mol] 2.6 • Fitting the intermolecular PES: Fit quality I Turbomole 15.06.2012
Fitting the intermolecular PES: Fit quality II Turbomole 15.06.2012
Fitting the intermolecular PES: Fit quality III Turbomole 15.06.2012
Attention!! Turbomole 15.06.2012
(SS1) • Fitting the intermolecular PES: Comparison to Experiment Turbomole They “confirmed” the existence of 1B employing HF, MP2 and B3LYP with 6-31G++(d,p) geometry optimizations. 15.06.2012
Intermolecular potential energy surface of the Acetylene - Acetylene dimer: Choice of geometries Turbomole 15.06.2012 C. Leforestier, A. Tekin, G. Jansen , M. Herman, J. Chem. Phys., 135, 234306, (2011)
Intermolecular potential energy surface of the Benzene - Benzene dimer: Choice of geometries Turbomole 15.06.2012
D6h Ring C2V • AcnBz cluster geometry optimizations: Test system Ac2Bz Turbomole 15.06.2012 a Fujii et al. J. Phys. Chem. A 108 (2004) 2652.
B A C A D C G E B F H • Structural comparison of Ac2Bz geometries Turbomole 15.06.2012
-10.79 -7.11 -10.73 -14.46 -19.30 -14.82 • Acn (n=3-6) clusters Turbomole 15.06.2012 (Energies are in mH)
-13.49 -13.51 -19.11 -24.43 -29.99 -30.52 • AcnBz (n=3-6) clusters Turbomole 15.06.2012 (Energies are in mH)
-10.53 -8.74 -8.31 • Ac1Bz2 clusters Turbomole A Angew. Chem. Int. Ed. 47 (2008) 10094 (Energies are in mH) 15.06.2012
-16.25 -13.36 -15.23 -15.39 -12.97 -15.37 -14.97 -15.49 -15.04 • Ac2Bz2 clusters Turbomole 15.06.2012 (Energies are in mH)
Intermolecular interactions [in kJ/mol] in water dimer Turbomole • Molpro energy results are in hartree. Convert them to kJ/mol by multiplying 2625.5. • DFT-SAPT results are in mH (miliHartree). Multiply them by 2.6255. • HF, MP2, SCS-MP2, CCSD and CCSD(T) results are in CCSD(T) outputs. • Only the SCS-MP2 interaction energy is not written in the output. To calculate it grep SCS-MP2 energies in the ccsd(t) outputs: • grep"SCS-MP2 total energy"water-ccsdt-avdz.out 15.06.2012 ESCS-MP2=(-152.52056354-(-76.25674577-76.25727940)) in hartree
Typical CCSD(T) output !RHF STATE 1.1 Energy MP2 total energy SCS-MP2 total energy !CCSD(T) total energy SETTING E_MB_CP_MP2 = -76.26175127 AU SETTING E_MB_CP_CCSD = -76.26944486 AU SETTING E_MB_CP_CCSDT = -76.27482665 AU SETTING E_INT_HF = -0.00605271 AU SETTING E_INT_MP2 = -0.00716194 AU SETTING E_INT_CCSD = -0.00685369 AU SETTING E_INT_CCSDT = -0.00710741 AU Turbomole 15.06.2012
Typical CCSD(T) output cp water-ccsdt-avdz.com water-ccsdt-avtz.com cp water-ccsdt-avtz.com water-ccsdt-avqz.com Turbomole 15.06.2012