How to use MOLPRO CCSD one electron density matrix for SAPT electrostatics

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

Theoretical introduction for SAPT electrostatics

As it is well known [1] the electrostatic energy $E^{(1)}_{\rm elst}$ can be written in terms of the electron densities of the isolated monomers A and B, $\rho_{\rm A}(\mbox{r})$ and $\rho_{\rm B}(\mbox{r})$ ,

$\begin{displaymath} E^{(1)}_{\rm elst} = \int\int \rho_{\rm A}(\mbox{r}_1) \rho_... ...}(\mbox{r}_2)v(\mbox{r}_1,\mbox{r}_2) d\mbox{r}_1 d\mbox{r}_2, \end{displaymath}$

(1)

where the generalized interaction potential $v(\mbox{r}_1,\mbox{r}_2)$ is defined as,

$\begin{displaymath} v(\mbox{r}_1,\mbox{r}_2)=\frac{1}{r_{12}} +\frac{1}{N_A}v_A(\mbox{r}_2) +\frac{1}{N_B}v_B(\mbox{r}_1) +\frac{1}{N_A N_B} n.r. \end{displaymath}$

(2)

where $v_A(\mbox{r}_2)$ describes the attraction of the electron 2 by the nuclei of the monomer A,

$\begin{displaymath} v_A(\mbox{r}_2)= -\sum_\alpha \frac{Z_\alpha}{r_{\alpha 2}} \end{displaymath}$

(3)

analogically, $v_A(\mbox{r}_1)$ describes the attraction of the electron 1 by the nuclei of the monomer B,

$\begin{displaymath} v_B(\mbox{r}_1)= -\sum_\beta \frac{Z_\beta}{r_{\beta 1}} \end{displaymath}$

(4)

and ``n.r." denotes the repulsion of nuclei belonging to the monomer A and B,

$\begin{displaymath} n.r.= \sum_\alpha \sum_\beta \frac{Z_\alpha Z_\beta}{R_{\alpha\beta}}. \end{displaymath}$

(5)

Therefore, the electrostatic energy can be easily calculated as any other first-order property of the molecule, provided that densities $\rho_{\rm A}(\mbox{r})$ and $\rho_{\rm B}(\mbox{r})$ 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:

$\begin{displaymath}\rho_{\rm X}(\mbox{r})=\rho^{\rm SCF}_{\rm X}(\mbox{r})+ \rho^{\rm corr}_{\rm X}(\mbox{r}), \end{displaymath}$

(6)

where

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

$\begin{displaymath} E^{(100)}_{\rm elst} = \int\int \rho^{\rm SCF}_{\rm A}(\mbox... ...}(\mbox{r}_2)v(\mbox{r}_1,\mbox{r}_2) d\mbox{r}_1 d\mbox{r}_2, \end{displaymath}$

(7)

and correlated parts:

$\begin{displaymath} E^{(10c)}_{\rm elst} = \int\int \rho^{\rm SCF}_{\rm A}(\mbox... ...}(\mbox{r}_2)v(\mbox{r}_1,\mbox{r}_2) d\mbox{r}_1 d\mbox{r}_2, \end{displaymath}$

(8)

$\begin{displaymath} E^{(1c0)}_{\rm elst} = \int\int \rho^{\rm corr}_{\rm A}(\mbo... ...}(\mbox{r}_2)v(\mbox{r}_1,\mbox{r}_2) d\mbox{r}_1 d\mbox{r}_2, \end{displaymath}$

(9)

$\begin{displaymath} E^{(1cc)}_{\rm elst} = \int\int \rho^{\rm corr}_{\rm A}(\mbo... ...}(\mbox{r}_2)v(\mbox{r}_1,\mbox{r}_2) d\mbox{r}_1 d\mbox{r}_2, \end{displaymath}$

(10)

Technical details of calculations with MOLPRO

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=record

after 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,record

Note that the nonrelaxed CCSD density is saved in the record

record-10

The 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,record

If 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

$\begin{displaymath}v_X^{\rm eff}= v_X+2J(\rho_X^{\rm SCF}) \end{displaymath}$ (11)
Obtain correlated densities for both monomers
Calculate $E^{(10c)}_{\rm elst}$ as

$\begin{displaymath}E^{(10c)}_{\rm elst} = \int \rho^{\rm corr}_{\rm B}(\mbox{r}_2)v_A^{\rm eff}(\mbox{r}_2) d\mbox{r}_2 \end{displaymath}$ (12)

and analogically $E^{(1c0)}_{\rm elst}$ ,
Calculate $E^{(1cc)}_{\rm elst}$ as

$\begin{displaymath}\int \rho^{\rm corr}_{\rm B}(\mbox{r}_2)J(\rho^{\rm corr}_{\rm A})(\mbox{r}_2) d\mbox{r}_2 \end{displaymath}$ (13)

Citations

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).

Sample MOLPRO script for calculation of SAPT electrostatics

The very latest version of the script is presented here (works for Molpro2020) here .

Below are archive versions (used in the explanation of what is done - versions of scripts differ mainly in the set number, where the needed density is stored in the record) 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 here . And here you find another script, which works with MOLPRO2006.1. The latter script works with MOLPRO2012.1 after one simple change:

symmetry,x

instead 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

enddo

Contact me

If you experience any problems with this script do not hesitate to contact me, my e-mail: tatiana.korona@tiger.chem.uw.edu.pl

Bibliography

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

Tatiana Korona 2007-03-08