Back to the main page of Tatiana Korona

** Contents **

- Theoretical introduction for SAPT electrostatics
- Technical details of calculations with MOLPRO
- Citations
- Sample MOLPRO script for calculation of SAPT electrostatics
- Contact me

where the generalized interaction potential is defined as,

(2) |

(3) |

(4) |

(5) |

Therefore, the electrostatic energy can be easily calculated as any other first-order property of the molecule, provided that densities and are known.

For the analysis of the electrostatic energy in a correlated picture it is
common to divide the total density into the SCF and correlated parts:

(6) |

Then the electrostatic energy (1) can be divided into uncorrelated

and correlated parts:

In the MOLPRO package the correlated densities on the CCSD, MP2, and QCISD level are available. In recent MOLPRO releases those densities can be saved using the command

save,density=recordafter the name of the corresponding method (from 2000.9 for MP2 and QCISD and from a patched 2002.5 for CCSD) and then read in MATROP using a command

load,arbitrary_name,DEN,recordNote that the nonrelaxed CCSD density is saved in the record

record-10The CCSD density is available from the version 2002.1, but is saved in another format to records 5050.2 (the CCSD density with the relaxation contribution) and 5040.2 (nonrelaxed). It can be then read in MATROP using a command

load,arbitrary_name,triang,recordIf the relaxed CCSD density is required (what should usually be the case), no frozen core can be used. The nonrelaxed CCSD density can be calculated with and without frozen core. MP2 and QCISD densities are always relaxed and can be calculated with and without frozen core. Note that the SAPT suite of codes does not allow for the frozen core anyway, so the correlated electrostatics should be calculated without frozen core for consistency with other SAPT corrections.

A simplified algorithm for obtaining the correlated part of the electrostatic energy in the MP2, QCISD, or CCSD approximations is given below:

- Calculate ``effective" potentials for the monomers A and B

(11) - Obtain correlated densities for both monomers
- Calculate
as

(12) - Calculate
as

(13)

Citations required if you use the correlated density from the MOLPRO package:

- CCSD: T. Korona, R. Moszynski, and B. Jeziorski, Mol. Phys.
**100**, 1723 (2002). - MP2: A. ElAzhary, G. Rauhut, P. Pulay, and H.-J. Werner, J. Chem. Phys.
**108**, 5185 (1998). - QCISD: G. Rauhut and H.-J. Werner, Phys. Chem. Chem. Phys.
**3**, 4853 (2001).

Below a sample MOLPRO script is given for calculating the CCSD electrostatics. A care is already taken in the script in order to detect the current MOLPRO version and to read the CCSD densities in MATROP correspondingly. You can also download it directly by clicking . And you find another script, which works with MOLPRO2006.1. The latter script works with MOLPRO2012.1 after one simple change:

symmetry,xinstead of

set zsymel=x

***,He-H2 correlated electrostatics from Molpro memory,4,m ! Thresholds should be set somewhat tighter than defaults ! if you calculate intermolecular interaction energies: gthresh,twoint=1.e-15,energy=1.e-7,coeff=1.e-6 ! just a small test basis basis=avdz ! Set geometry r=4.0 bohr rhh=1.4 bohr geometry={ he,, 0., 0. , r he h1,, -rhh/2, 0. , 0. h h2,, rhh/2, 0. , 0. h} ! Nuclear repulsion term (only necessary for E(100)) inr=2*2/sqrt((rhh/2)**2+r**2) ! Choose the method(s): method=[ccsd] ! available methods: ccsd,qcisd,mp2 ! For the final table: Term=['E(100)','E(1c0)','E(10c)','E(1cc)','E(1)(total)','E(1)(corr)'] do i=1,#method ! A MONOMER A MONOMER A MONOMER A MONOMER A MONOMER A MONOMER A MONOMER ! Do SCF for the A monomer (of course with the ghost functions on the monomer B) scfa=2110.2 ! here SCF results for A will be saved dummy,h1,h2 ! Accuracy of SCF should also be set tighter than default: hf;accu,13;save,scfa ! Use MATROP to calculate the effective potential for the monomer A: matrop load,EPOT ! read the bare potential for A (v_A) save,EPOT,5703.2,triang ! save the bare potential for A in the record 5703.2 ! (it will be overwritten in the calculations for the monomer B) load,DEN,DEN,scfa ! read the SCF density for A coul,J,DEN ! calculate the Coulomb operator with the SCF A density add,WA,1,EPOT,2,J ! calculate the effective potential for A save,WA,5701.2,triang ! and save it in the record 5701.2 ! Calculate the correlated density for the monomer A ! and save it to the record 5051.2 $method(i);core,0;cphf,1 expec save,density=5051.2 if($method(i).eq.CCSD.and.$version.lt.2002.5) then data,copy,5050.2,5051.2 ! with relaxation contribution (better) data,copy,5040.2,5041.2 ! nonrelaxed (worse, but can be ! calculated with frozen core) endif ! B MONOMER B MONOMER B MONOMER B MONOMER B MONOMER B MONOMER B MONOMER dummy,he scfb=2111.2 ! here SCF results for B will be saved hf;accu,13;save,scfb ! Use MATROP to calculate the effective potential for the monomer B: matrop load,EPOT load,DEN,DEN,scfb coul,J,DEN add,WB,1,EPOT,2,J save,WB,5702.2,triang ! Calculate the correlated density for the monomer B ! and save it to the record 5052.2 $method(i);core,0;cphf,1 expec save,density=5052.2 if($method(i).eq.CCSD.and.$version.lt.2002.5) then data,copy,5050.2,5052.2 data,copy,5040.2,5042.2 endif ! CALCULATION of the ELECTROSTATIC ENERGY ! this is the begin of MATROP if($method(i).eq.CCSD.and.$version.lt.2002.5) then matrop load,TOTDENA,triang,5051.2 ! read the correlated density for A load,TOTDENB,triang,5052.2 ! read the correlated density for B load,WA,triang,5701.2 ! read the effective potential for A load,WB,triang,5702.2 ! read the effective potential for B load,DENA,DEN,scfA ! read the SCF density for A load,DENB,DEN,scfB ! read the SCF density for B add,CDENA,1,TOTDENA,-1,DENA ! calculate the correlated part of density for A add,CDENB,1,TOTDENB,-1,DENB ! calculate the correlated part of density for B ! calculate E(100)elst (just for check with SAPT) ! in two ways (see below) load,EPOT ! load bare potential for monomer B trace,e100B,DENA,EPOT trace,e100A,DENB,WA load,EPOT,triang,5703.2 ! load bare potential for monomer A trace,e100A1,DENB,EPOT trace,e100B1,DENA,WB ! (just for check: the trace of the SCF density should be equal to ! the number of electrons and the trace of the correlated part should be 0) load,S ! load the overlap matrix trace,na,DENA,S trace,nb,DENB,S trace,nap,CDENA,S trace,nbp,CDENB,S ! The main point - calculation of the correlated electrostatics: trace,e1n0,CDENA,WB trace,e10n,CDENB,WA coul,J,CDENA trace,e1nn,J,CDENB ! this is the end of MATROP else ! dependent on molpro version matrop load,TOTDENA,den,5051.2 ! read the correlated density for A load,TOTDENB,den,5052.2 ! read the correlated density for B load,WA,triang,5701.2 ! read the effective potential for A load,WB,triang,5702.2 ! read the effective potential for B load,DENA,DEN,scfA ! read the SCF density for A load,DENB,DEN,scfB ! read the SCF density for B add,CDENA,1,TOTDENA,-1,DENA ! calculate the correlated part of density for A add,CDENB,1,TOTDENB,-1,DENB ! calculate the correlated part of density for B ! calculate E(100)elst (just for check with SAPT) ! in two ways (see below) load,EPOT ! load bare potential for monomer B trace,e100B,DENA,EPOT trace,e100A,DENB,WA load,EPOT,triang,5703.2 ! load bare potential for monomer A trace,e100A1,DENB,EPOT trace,e100B1,DENA,WB ! (just for check: the trace of the SCF density should be equal to ! the number of electrons and the trace of the correlated part should be 0) load,S ! load the overlap matrix trace,na,DENA,S trace,nb,DENB,S trace,nap,CDENA,S trace,nbp,CDENB,S ! The main point - calculation of the correlated electrostatics: trace,e1n0,CDENA,WB trace,e10n,CDENB,WA coul,J,CDENA trace,e1nn,J,CDENB ! this is the end of MATROP endif ! dependent on molpro version ! Now extract the results: ! the SCF-correlated part (e10n) and the correlated-SCF part (e1n0) are ready ! the correlated-correlated part: e1nn=2*e1nn ! the SCF-SCF part (calculated in two ways, both should give the same number!) e100=INR+e100A+e100B e100p=INR+e100A1+e100B1 if(abs(e100-e100p).gt.1.e-8) then text, Warning: SCF check failed,$e100,$e100p endif ! the total correlation contribution: e1c=e10n+e1n0+e1nn ! the total electrostatic contribution: e1t=e100+e1c ! hartree - > Kelvin e100k=e100*tokelvin e1n0k=e1n0*tokelvin e10nk=e10n*tokelvin e1nnk=e1nn*tokelvin e1ck=e1c*tokelvin e1tk=e1t*tokelvin ! hartree - > cm-1 e100c=e100*tocm e1n0c=e1n0*tocm e10nc=e10n*tocm e1nnc=e1nn*tocm e1cc=e1c*tocm e1tc=e1t*tocm ! Make the final table: ! final results in microhartrees: n=6 skala=10**n eau(1)=e100*skala eau(2)=e1n0*skala eau(3)=e10n*skala eau(4)=e1nn*skala eau(5)=e1t*skala eau(6)=e1c*skala eke(1)=e100k eke(2)=e1n0k eke(3)=e10nk eke(4)=e1nnk eke(5)=e1tk eke(6)=e1ck ecm(1)=e100c ecm(2)=e1n0c ecm(3)=e10nc ecm(4)=e1nnc ecm(5)=e1tc ecm(6)=e1cc text, FINAL RESULTS table,Term,eau digits,,6 title, Electrostatic interaction energy in the $method approximation title, energies in hartree/10^$n table,Term,eke digits,,6 title, Electrostatic interaction energy in the $method approximation title, energies in Kelvin table,Term,ecm digits,,6 title, Electrostatic interaction energy in the $method approximation title, energies in cm^-1 enddoBack to top of the page

**Bibliography **

1
B. Jeziorski, R. Moszynski, and K. Szalewicz, Chem. Rev. **94**, 1887
(1994).

Tatiana Korona 2007-03-08