Gerald Knizia knizia at theochem.uni-stuttgart.de
Wed Feb 24 07:59:57 GMT 2010

```On Wednesday 24 February 2010 00:10, Yulia Kalugina wrote:
> Do anyone know how to calculate the dipole-hexadecapole polarizability at
> CCSD(t) level? what keywords should be used?

Unfortunately, Molpro is still somewhat lacking in its facilities to calculate
properties. The highest level at which Molpro can calculate first order
properties analytically (see ``PROPERTIES AND EXPECTATION VALUES'' in the
manual) are closed-shell CCSD (without T) and MRCI. For second order
properties (like polarizabilities) there is nothing.

You can, however, use the finite field approach to evaluate these properties
with arbitrary electronic structure methods as numerical second derivatives.
If correct field strengths are chosen, this approach is very accurate. If it
is sufficient for hexadecapole polarizabilities is a different matter of
course... (should work fine for dipole and quadrupole, tho).
(I attached a file which I used to find optimal field strengths for the dipole
polarizability of the carbon atom.)
Similarly, to get the polarizabilities when you have first-order moments, a
first order numerical derivative with respect to the finite field strength
would have to be evaluated with a similar approach.
--
Gerald Knizia
-------------- next part --------------
memory,64,m;

gthresh,energy=1e-14,coeff=1e-7
gthresh,twoint=1e-15

basis=davdz

! geometry={B};
! {rhf; accu,20; wf,sym=2,spin=1; save,2110.2}
geometry={C};
{rhf; accu,20; wf,sym=4,spin=2; open,1.2,1.3; save,2110.2}
! geometry={N};
! {rhf; accu,20; wf,sym=8,spin=3; save,2110.2}
! geometry={O};
! {rhf; accu,20; wf,sym=4,spin=2; save,2110.2}
! geometry={F};
! {rhf; accu,20; wf,sym=2,spin=1; save,2110.2}
! geometry={B};
! {rhf; accu,20; wf,sym=2,spin=1; start,atden}

\$basislist = ['DAVDZ','DAVTZ','DAVQZ','DAV5Z']

_CC_NORM_MAX=2.0

fss = [.032,.01,.0032,.001,.00032,.0001,.000032,.00001,.0000032,.000001,.00000032,.0000001]

do i = 1,#fss
fs = fss(i)
!    bases(i)=\$basislist(i)
bases(i)='DAVTZ'
basis=\$bases(i)

field
set,zsymel=[x,y,z]
{rhf; accu,20; wf,sym=4,spin=2; open,1.2,1.3; start,atden}
{uccsd(t); maxit,200; diis,1,1,8}
EShift = 0
EShift = ENERGY
En(i) = ENERGY - EShift

set,zsymel=[x,y]
dip,0.00,0.00,fs
{rhf; accu,20; wf,spin=2,sym=4; occ,2,1,1,0; open,1.2,1.3; start,atden}
{uccsd(t); maxit,200; diis,1,1,8}
EpP(i) = ENERGY - EShift
EmP(i) = ENERGY - EShift

dip,0.00,0.00,2*fs
{rhf; accu,20; wf,spin=2,sym=4; occ,2,1,1,0; open,1.2,1.3; start,atden}
{uccsd(t); maxit,200; diis,1,1,8}
EpP2(i) = ENERGY - EShift
EmP2(i) = ENERGY - EShift

set,zsymel=[x,y]
dip,0.00,0.00,fs
{rhf; accu,20; wf,spin=2,sym=2; occ,3,1,0,0; open,3.1,1.2; start,atden}
{uccsd(t); maxit,200; diis,1,1,8}
EpX(i) = ENERGY - EShift
EmX(i) = ENERGY - EShift

dip,0.00,0.00,2*fs
{rhf; accu,20; wf,spin=2,sym=2; occ,3,1,0,0; open,3.1,1.2; start,atden}
{uccsd(t); maxit,200; diis,1,1,8}
EpX2(i) = ENERGY - EShift
EmX2(i) = ENERGY - EShift

table,bases,En,EpP,EpP2,EpX,EpX2

! first order formulas (need two calculations less per point, but need lower fs value)
!    PolP(i) = (2*En(i)-EpP(i)-EmP(i))/(fs**2)
!    PolX(i) = (2*En(i)-EpX(i)-EmX(i))/(fs**2)

! second order formulas (more stable with respect to variations in fs)
PolP(i) = -((-2*1.0*EpP2(i)+2*16.0*EpP(i))-30.0*En(i) )/(12.0*fs**2)
PolX(i) = -((-2*1.0*EpX2(i)+2*16.0*EpX(i))-30.0*En(i) )/(12.0*fs**2)
PolM(i) = (2.0/3.0) * PolX(i) + (1.0/3.0) * PolP(i)
PolA(i) = PolX(i) - PolP(i)

text,PolP: atom/field polarization parallel
text,PolX: atom/field polarization crossed
text,PolM: mean dipole polarization
text,PolA: dipole polarization anisotropy
table,fss,PolP,PolX,PolM,PolA
end do
-------------- next part --------------
memory,64,m;

gthresh,energy=1e-14,coeff=1e-7
gthresh,twoint=1e-15

basis=davtz
geometry={C};
{rhf; accu,20; wf,sym=4,spin=2; open,1.2,1.3; save,2110.2}
mcscf
{ci; expec,dm,qm,MLTP6}
```