<div><font size="4" face="verdana,sans-serif">Dear Ms/Mr,</font></div><div><font size="4" face="verdana,sans-serif">         Now, I am obsessed by one question, can you help me solve the problem? Thank you very much!</font></div>
<div><font size="4" face="verdana,sans-serif">         It is about calculating the polarizability of one molecule by differential method. Based on a case (H2O molecule) in Molpro manual, the input file is displayed as following: </font></div>
<div><font size="4" face="verdana,sans-serif">"    </font></div><div><pre><font size="4" face="verdana,sans-serif"><font face="simsun,serif">***,H2O finite field calculations

r=1.85,theta=104                   !set geometry parameters
geometry={O;                       !z-matrix input
          H1,O,r;
          H2,O,r,H1,theta}
basis=avtz                         !define default basis
field=[0,0.005,-0.005]             !define finite field strengths
$method=[hf,mp4,ccsd(t),casscf,mrci]

k=0
do i=1,#field                      !loop over fields
  dip,,,field(i)                   !add finite field to H
  do m=1,#method                   !loop over methods
    k=k+1
    $method(m)                     !calculate energy
    e(k)=energy                    !save energy</font>
<font face="simsun,serif">  enddo
enddo

k=0
n=#method
do m=1,#method
  k=k+1
  energ(m)=e(k)
  dipmz(m)=(e(k+n)-e(k+2*n))/(field(2)-field(3))   !dipole moment as first energy derivative
  dpolz(m)=(e(k+n)+e(k+2*n)-2*e(k))/((field(2)-field(1))*(field(3)-field(1)))  !polarizability as second der.
</font><font face="simsun,serif">enddo

table,method,energ,dipmz,dpolz
title,results for H2O, r=$R, theta=$theta, basis=$basis
---</font>
</font></pre><pre><font size="4" face="verdana,sans-serif">"</font></pre></div><div><font size="4" face="verdana,sans-serif"> The result is as following:</font></div><div><font size="4" face="Verdana">"</font></div>
<div><font size="4" face="simsun,serif"> METHOD        ENERG        DIPMZ        DPOLZ<br> HF        -76.05828804   0.79028512   8.66914177<br> MP4       -76.34319951   0.72435361   9.84241227<br> CCSD(T)   -76.34178065   0.72957181   9.67032932<br>
 CASSCF    -76.11356760   0.72500885   8.89662369<br> MRCI      -76.32839170   0.70823783   9.20622885</font></div><div><font size="4" face="Verdana">"<br></font></div><div><font size="4" face="Verdana">As showing above, the value of polarizability in z direction is positive. But when I calculate the z direction component of the polarizability by adding a finite electronic field in z direction, and then compute the derivative of dipole with respect to the field, the value of this derivative is found be negative but have the same magnitude by and large. The input file of normal calculation is:</font></div>
<font size="4" face="Verdana"><div><div><font size="4" face="Verdana">"</font></div><div><font size="4" face="simsun,serif">***,H2O finite field calculations</font></div><div><font size="4" face="simsun,serif">r=1.85,theta=104                                !-->set geometry parameters<br>
geometry={O;                   !-->z-matrix input<br>          H1,O,r;<br>          H2,O,r,H1,theta}<br>basis=avtz                                            !-->define default basis</font></div><div><font size="4" face="simsun,serif">{hf}<br>
{mrci;dm;<br>natorb,2140.2}</font></div></div></font><div></div><div><font size="4" face="Verdana">"</font></div><div><font size="4" face="Verdana">The result is:</font></div><div><font size="4" face="Verdana">"</font></div>
<div><font size="4" face="simsun,serif"> !RHF STATE 1.1 Dipole moment           0.00000000     0.00000000     0.79032209</font></div><div><font size="4" face="simsun,serif"> !CI(SD) STATE 1.1 Dipole moment        0.00000000     0.00000000     0.74610325</font></div>
<div><font size="4" face="Verdana">"</font></div><div><font size="4" face="Verdana">Add the field (0.001):</font></div><div><font size="4" face="Verdana">"</font></div><div><font size="4" face="Verdana"><div><font size="4" face="simsun,serif">***,H2O finite field calculations</font></div>
<div><font size="4" face="simsun,serif">r=1.85,theta=104                                                 !-->set geometry parameters<br>geometry={O;                                        !-->z-matrix input<br>          H1,O,r;<br>
          H2,O,r,H1,theta}</font></div><div><font size="4" face="simsun,serif">basis=avtz                         !-->define default basis</font></div><div><font size="4" face="simsun,serif">field=0.001                        !-->define finite field strengths<br>
dip,,,field</font></div><div><font size="4" face="simsun,serif">{hf} <br>{mrci;dm;<br>natorb,2140.2</font></div></font></div><div><font size="4" face="Verdana">"</font></div><div><font size="4" face="Verdana">The result is:</font></div>
<div><font size="4" face="Verdana">"</font></div><div><font size="4" face="simsun,serif">!RHF STATE 1.1 Dipole moment                    0.00000000     0.00000000     0.78165042</font></div><div><font size="4" face="simsun,serif">!CI(SD) STATE 1.1 Dipole moment              0.00000000     0.00000000     0.73678704</font></div>
<div> </div><div><font size="4" face="Verdana">"</font></div><div><font size="4" face="Verdana">So, for the HF method, polzz<font face="verdana,sans-serif">=(0.78165042-0.79032209)/0.001=-8.67167</font></font></div><div>
<font size="4" face="Verdana">and, for the CD(SD),    <font face="verdana,sans-serif">polzz=(0.73678704-0.74610325)/0.001=-9.31621</font></font></div><div><font size="4" face="Verdana">Actually, this derivative is also the value of polarizability [I also found this definition in some references pol_zz=d(dip_z)/d(E_z) ]. So why the result of this method is negative? Or there're some details in Molpro calculations which I neglected? Thanks a lot!</font></div>
<div><font size="4" face="Verdana"></font> </div><div><font size="4" face="Verdana">Best regards,</font></div><div><font size="4" face="Verdana">Changjian     </font></div><font size="4" face="verdana,sans-serif">
</font>