CCSD density matrix

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.

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)
set_global_option('DERTYPE','FIRST')
set_global_option('WFN','CCSD')
ccdensity(ccwfn)

Da_so = ccwfn.Da()
Da_so.print_out()

The 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.

Thank you for your reply! That’s just perfect!
I think I just have one last question:
I tried with this code and found that it only worked with ‘Frozen core’ false, which makes sense because I think that it is still under the control of ‘gradient’. Because

So I am wondering if I can call this function with ‘Frozen core’ on? If so, how can I do it?

Unfortunately, there is no way of getting relaxed CCSD densities/properties with frozen core as of now since we haven’t coded frozen core contributions to orbital relaxations. Frozen core CC gradients is definitely on our TODO list and we would try to implement it soon in PSI4.

Thank you very much for all the reply and explanation! They really helped a lot!

A post was split to a new topic: Can I get a CCSD(T) relaxed density matrix to hack a property calc?