You can write a little hack to calculate CCSD dipoles and quadrupoles using relaxed densities without having to calculate gradients:
dp, ccwfn = properties('ccsd', properties=['dipole'],return_wfn=True)
Da_so = ccwfn.Da()
ccdensity(ccwfn) call should give you the relaxed dipole which should be identical to the one obtained from a gradient calculation. Also, the density matrix should match as well.