CCSD density matrix

Hi,

I managed to get the one particle density matrix from a CCSD calculation, but I’m not sure what basis it is using. Is it expressed in the basis set or the fock molecular orbitals basis?

the input script is just doing a scf run, saving the wavefunction, and feeding it to a CCSD gradient run. I print out ccsdwfn.Da().toarray with the dense option enabled.

I apprecite the help.

Actually, it might just be the AO-basis SCF density. Unfortunately, CCSD density doesn’t get passed on to the ccsdwfn and hence reference wfn’s density is printed when you do ccsdwfn.Da().toarray. If you set the print flag of the ccdensity module to be 5 in your input file, you can get the onepdm printed in the output. You can add the line below in your input file.
set ccdensity{ print 5 }

1 Like

I see.
I had that problem before, where the Da of the SCF wfn was the same as the Da of the CCSD wfn. After searching around for solutions, I found that calling gradient() instead of energy() forces the code to update the density matrix.
Infact, I tested it with a simple H2 molecule and a fake basis set with just two s orbitals, and the SCF and CCSD density matrixes were different. Now I’m wondering:

  1. are they the same matrix in a different basis (fock-AO, ccsd-MO)?
  2. are they different?
  3. are they both in AO?

@felix00 The wfn returned from CC calculations will hold the ccsd density if it was calculated. Since the density is not required to compute the energy it will not be constructed in a simple energy('ccsd') calculation. As you have seen requesting a gradient will require the density to be constructed so the wfn that is returned will be updated and wfn.Da() will now return the ccsd density.

In all cases the density returned by wfn.Da() is in the SO basis, you can obtain the MO basis density by transforming using wfn.Ca(). Or using the wfn.D_subset_helper() function.

So to answer your three questions they are different matrices the first would be the SCF density in the SO basis, the second is the CCSD density also in the SO basis.

Thanks a lot for the info.
I am not so familiar with group theory, so just to be clear, if the molecule has C1 symmetry (in my limited understanding that means “a mess with no symmetries”), the SO is equivalent to the AO? or MO?
In any case, is there a quick way to get it in AO?

thanks a lot!

Yes if the molecule has C1 symmetry, or if the symmetry was set to C1 manually SO==AO.

@amjames we should change the name of the density appropriately. If I print out the density returned from a ccsd gradient calculation wfn.Da().print_out(), it still says SCF density (Symmetry 0) ## in the output file, which is misleading.

Just out of curiosity… how would one go about printing the two-particle DM after a CCSD calculation?

Hi, I really appreciate your reply here. They are quite helpful to me! And I still got a problem with getting access to the ccdensity matrix during the calculation. Could you please give me some specific instruction about that? What I did basically looks like this:
molecule mol {
0 1
C .0000030150 .0000036400 .0000139900
H 1.0962824850 .0002186300 .0088333700
H -.3657450550 1.0334718700 .0088329700
H -.3725226950 -.5271066900 .8861702200
H -.3580328250 -.5066056500 -.9039205000

symmetry c1
no_reorient
no_com
}

set {
basis d-aug-cc-pVDZ
guess sad
scf_type direct
freeze_core true
}
property(‘ccsd’, properties=[‘dipole’,‘quadrupole’,‘polarizability’])

mints = psi4.core.MintsHelper(ccsd.wfn.basisset())
print_out(mints.Da())

But it told me that “ccsd” was not defined… What should I do? Thank you so much for your help…

Hi,

seems you are accessing ccsd.wfn.basisset() but the base object ccsd is not defined as a python variable. I doubt it is defined as global either, but the devs will know better.
In my case I hacked it this way:

# — HARTREE-FOCK ENERGY AND WF —#
E, wf = energy(‘scf’, return_wfn=True)
print “HF energy:”, E

# — CCSD starting from the HFWF —#
grad, ccwf = gradient(‘ccsd’, return_wfn=True, ref_wfn=wf)
oeprop(ccwf, ‘DIPOLE’, title=‘MYDIPOLE’)

now it is possible to access the D matrix of the ccwf object (CCSD wave function). This will return the symmetric matrix:
ccwf.Da().to_array(copy=False, dense=True)

When I started playing with this stuff, using gradient was the only (easy) way to force it update the DM. I think the property function could not return the wavefunction object. Unfortunately gradients with frozen_core were not supported so I had to do the full CCSD, and I kinda wanted the forces anyway.
Now I can live without forces and I’d love if the devs could suggest a script to get DM in AO basis from a CCSD frozen_core run.

1 Like

Thank you so much for your reply!

Wait… Is there a way to update the DM when we just calculate the ccsd energy with frozen core?

I am trying to figure that out in another thread I started since I got some errors when trying to run larger systems.
Perhaps with:
ccwf = properties(‘ccsd’, properties=[‘dipole’], return_wfn=True)
the DM will be updated. Need devs to confirm or test.

@felix00 calling properties(‘ccsd’, properties=[‘dipole’], return_wfn=True) will indeed
update the DM. This way you don’t need to calculate gradients. However, if you require
the relaxed 1PDM, then you do need to use the gradients approach instead. Also, you can’t
use frozen core for calculating relaxed densities as it has not been implemented in PSI4 .

Thank you for your reply. But don’t I need to relax the density matrix when I am calculating the “polarizability”?

I did some test calculations with H2O under 6-31G and got the following result(just the first five items of the first rows of DM’s were showed ):

SCF(FC=0) 1.04247606429881 -0.09150908508311 -0.00 -0.00 -0.02016733000915
gradient(FC=0) 1.04220942491249 -0.09128336492263 -0.00 -0.00 -0.01913359194153
SCF(FC=1) 1.04247606429881 -0.09150908508311 -0.00 -0.00 -0.02016733000915
dipole(FC=1) 1.04225500411870 -0.09151171578025 -0.00 0.00 -0.01910792905284
polarizability(FC=1) 1.04247604635203 -0.09150899173888 -0.00 -0.00 -0.02016731398249

Well from the result above I do see some differences between them. However, here are my questions:

  1. Since they look so similar to each other(they are the same until the fourth digital number), is the DM really updated in gradient/dipole calculations?
  2. The DM form polarizability calculation looks the same as that form the SCF calculation w/o CF. Why?
  3. I noticed that you used the keyword “properties” as well as in some Github discussion. However, I got an error when I used it in my input file… Does it mean that I am not using the latest version Psi4?

For clearance, here are the two input files I used:


memory 90000 mb

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H 0.75718909 0.00000000 -0.58651475
H -0.75718909 0.00000000 -0.58651475

symmetry c1
no_reorient
no_com
}

set {
basis 6-31G
guess sad
scf_type direct
freeze_core false
}

e, wfn=energy(‘SCF’, return_wfn=True)
mints = psi4.core.MintsHelper(wfn.basisset())

import numpy
numpy.set_printoptions(threshold=numpy.nan)
print_out("\nOverlap matrix")
mints.ao_overlap().print_out()

Da_so = wfn.Da()
Da_so.print_out()

grad, ccwfn=gradient(‘ccsd’,return_wfn=True,ref_wfn=wfn)
ccDa_so = ccwfn.Da()
ccDa_so.print_out()


memory 90000 mb

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H 0.75718909 0.00000000 -0.58651475
H -0.75718909 0.00000000 -0.58651475

symmetry c1
no_reorient
no_com
}

set {
basis 6-31G
guess sad
scf_type direct
freeze_core True
}

e, wfn=energy(‘SCF’, return_wfn=True)
mints = psi4.core.MintsHelper(wfn.basisset())

import numpy
numpy.set_printoptions(threshold=numpy.nan)
print_out("\nOverlap matrix")
mints.ao_overlap().print_out()

Da_so = wfn.Da()
Da_so.print_out()

dp, ccwfn=property(‘ccsd’, properties=[‘dipole’], return_wfn=True)
ccDa_so = ccwfn.Da()
ccDa_so.print_out()

pl, ccwfn1=property(‘ccsd’, properties=[‘polarizability’], return_wfn=True)
ccDa_so = ccwfn1.Da()
ccDa_so.print_out()


I am sorry that I have so many questions… I am so confused about the DM now. Any kind help would be strongly appreciated.

You are right! You don’t want to include orbital relaxation or relax the density while calculating ccsd polarizabilties as it results in spurious poles appearing in the linear response function.

The DMs should indeed be updated in both ccsd gradient and dipole calculations. Here, we are adding the correlation density part to the SCF density (MO basis) and then back-transforming them to the ao basis. The total density looks very similar to the SCF density since the correlated density elements would be very small in magnitude.

Yes, the latest version of psi4 uses properties instead of property. You can do a psi4 --version to know which version you are using.

You are right. Actually calculating ccsd polarizabilities does not require ground state coupled cluster densities. So, the density in this case would be just the SCF density. Let me know if you have any more questions.

Now I am a bit confused about DM too : )
what is the practical difference between relaxed and not?
Let’s say I want a description of the electron density… which one should I use?
In my previous calculations I computed Mulliken charges using the density that came our after gradients (I guess relaxed). Even though Mulliken/Bader/Löwdin charges are not quantum observables and the results should be taken with a bit of salt, which DM should one use?

In order to calculate properties based on derivative of the energy wrt a perturbation, like dipoles with non-variational methods like MP2, CCSD one should ideally use the relaxed density approach as it will be identical to the finite difference calculations of energies in the presence of the perturbation. However, I am not sure if its worth the effort to calculate dipoles using relaxed density in CCSD though as the differences might be small for practical purposes. Look at this paper by Dr. Bartlett http://aip.scitation.org/doi/abs/10.1063/1.453596

For properties based on an expectation value like Mulliken charges, I don’t see the reasoning behind using a relaxed density. (Relaxed CCSD densities have been shown to have unphysical occupation numbers) So, using unrelaxed densities in my view is fine for calculating Mulliken charges.

As you already know, calculating orbital relaxation in ccsd requires gradient(...) calculations which are very expensive compared to a properties(...) calculation. So by default we don’t relax the DM when calculating any CCSD one electron property using properties(...) as its not worth the effort.

Thank you for your reply! I am now a little bit clear about the stuff. Now I am wondering if there is a single function that I can call to relax given wavefunction and density matrix. Because what I want to do is:

pl, ccsd_wfn = property(‘ccsd’, properties=[‘polarizability’],return_wfn=True)
(The unrelaxed ccsd wfn and DM are calculated)

ccsd_wfn_relaxed = some_function(ccsd_wfn)
(Relax the ccsd wfn and DM)

oeprop(ccwf, ‘DIPOLE’,‘QUADRUPOLE’, title=‘MYMULTIPOLE’)
(The relaxed ccsd wfn and DM are used to calculate dipole and quadrupole)


Otherwise, I need to replace ‘some_function’ above with ‘gradient’, which means I have to calculate ccsd again. This will make my calculation at least as twice expensive as the above calculation…

Thanks a lot for the explanation.