How to use ccsd to get analytical polarizability

How To
#1

I know analytic gradients are available for the CCSD, and the coupled-cluster package in Q-CHEM can calculate properties of target EOM states including permanent dipoles, angular momentum projections, static polarizabilities. So I want to know if I can get analytical polarizability from CCSD? I have tried to use the key word

$rem
   JOBTYPE                      POLARIZABILITY
   METHOD                       CCSD
   BASIS                            aug-cc-pVDZ
   MEM_TOTAL                120000
   RESPONSE_POLAR    0
$end

but I find the output polarizability is actually the analytical polarizability from HF method. So can you give me an example of the input file if Q-chem can give analytical polarizability of CCSD?

#2

See Section 7.10.20.8 of the manual. Examples are in Section 7.10.21.
https://manual.q-chem.com/latest/sec_EOMGRAD.html

#3

I have tried to get the analytical polarizability using two methods. However, the results are different and neither is the same as the numerical results. So I am puzzled to use which method. Here are my input,

$rem
   METHOD            CCSD
   BASIS             aug-cc-pVDZ
   MEM_TOTAL         200000
   CC_REF_PROP       1
   CC_POL            2
   CC_DIIS_SIZE      15
   CC_FULLRESPONSE   1
   SYM_IGNORE        TRUE
   no_reorient       true
$end   

$rem
   METHOD            CCSD
   BASIS             aug-cc-pVDZ
   MEM_TOTAL         200000
   CC_REF_PROP       1
   CC_POL            3
   CC_DIIS_SIZE      15
   CC_FULLRESPONSE   0
   SYM_IGNORE        TRUE
   no_reorient       true
$end   

$rem
   JOBTYPE           POLARIZABILITY
   METHOD            CCSD
   BASIS             aug-cc-pVDZ
   MEM_TOTAL         100000
   RESPONSE_POLAR    -1
   FDIFF_STEPSIZE    10
   SYM_IGNORE        TRUE
   no_reorient       true
$end

the results are,

 Polarizability Tensor ALPHA_IJ, I,J=X,Y,Z (a.u.):
 ------------------------------------------------------------------------------
         44.82016642   -0.00040250   -0.00022886
         -0.00040209   46.87450347    0.00056227
         -0.00022878    0.00055402   33.88883125
 ------------------------------------------------------------------------------
 Isotropic polarizability: 1/3 * (a_xx + a_yy + a_zz)
 Isotropic polarizability in a.u. = 41.86116705
 Anisotropic polarizability: SQRT(0.5*((a_xx-a_yy)^2 + (a_yy-a_zz)^2 + (a_zz-a_xx)^2))
 Anisotropic polarizability in a.u. = 12.09012143

Polarizability Tensor ALPHA_IJ, I,J=X,Y,Z (a.u.):
 ------------------------------------------------------------------------------
        -45.25278662    0.00049956    0.00022950
          0.00049956  -47.70229233   -0.00063772
          0.00022950   -0.00063772  -34.39252191
 ------------------------------------------------------------------------------
 Isotropic polarizability: 1/3 * (a_xx + a_yy + a_zz)
 Isotropic polarizability in a.u. = -42.44920029
 Anisotropic polarizability: SQRT(0.5*((a_xx-a_yy)^2 + (a_yy-a_zz)^2 + (a_zz-a_xx)^2))
 Anisotropic polarizability in a.u. = 12.26978843

and

 Polarizability tensor      [a.u.]
     43.9396279     -0.0007206     -0.0002346
     -0.0003291     45.8813873      0.0006981
     -0.0002276      0.0009505     33.3528175
  Principal components       [a.u.]
     33.3528175     43.9396276     45.8813876
  Isotropic polarizability   [a.u.]
     41.0579442
  Anisotropic polarizability [a.u.]
     11.6793846

I know the two analytical polarizability results should be different as the incomplete basis, but I want to ask if I can trust the numerical resut?