[molpro-user] dipole-hexadecapole polarizability

[Gerald Knizia]

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}