I recently used this piece of code to calculate the unrelaxed dipole moment (with either MP2 or CCSD) by finite field method. However, the result seems fairly off in magnitude. Is there anything wrong I did?
Check if scf_wfn.Fa() changes after Perturb_matrix.add(dipole[0]). My intuition is that you’ve added to a copy of the Fock matrix but haven’t actually changed the Fock matrix.
Figuring out if this actually does what you want will depend on how, exactly, the code in question determines its integrals… What Psi module are you trying to compute properties from? We have density fitted and non-density fitted versions of these methods.
Hi, thanks for replying. So I’m sure that the Fock matrix is changed because I used the same code to calculate the dipole polarizability by finite field method before and the resulting polarizability looks fine. The weird thing is that when I used the code to calculate the dipole moment, the results are not reasonable at all.
I can reproduce your unrelaxed and relaxed energies on master Psi. The CC code does pull the Fock integrals from the wavefunction’s Fock method. The CC code doesn’t seem to assume the occupied/virtual block of the Fock matrix is zero. In other words, there went all of my initial ideas.
Your unrelaxed dipole test computes the electronic dipole, and your finite field test (using perturb_with) computes the total dipole, electronic + nuclear. perturb_with knows to include the nuclear dipole terms with nuclear repulsion. Just changing the Fock matrix will not, so of course the two will have a large disagreement. There may be something else amiss. Let us know if accounting for that fixes your problem, or if there is still something else suspicious.
I will warn you that I am not familiar with the coupled cluster code, so you may get to the limits of my knowledge very quickly.
Thanks for the reply! I’ll try your suggestions, which is literally sum the calculated electronic dipole moment (by finite field method) and the nuclear dipole (from the output) together, right?
I’ll let you know if it woks once I test this suggestion! Thanks!
Right. You can also get the nuclear dipole from mol.nuclear_dipole(), so you can add the two together in your Psi input file.
And just to be clear, the two are already added together in the relaxed tests. It’s the unrelaxed tests where you modify the Fock matrix that need this added.
In terms of not working, the resulting numerical (electronic) dipole moment by using the “perturb_matrix” is (X: 1.039316426685904, Y: -0.785471768812530, Z: 0.000000000383653) in a.u. which is different from the analytical electronic dipole moment calculated by linear response CCSD calculation, which is (X: 1.7029, Y: -0.9185, Z: 0.0000) in a.u…
For comparison, the numerical (total) dipole moment by using “with_perturb” is (X: -0.500138727232233, Y: 0.538022021828692, Z: 0.000000005635737) in a.u., which is consistent with the analytical result (X: -0.5494, Y: 0.5656, Z: 0.0000) in a.u. (the little difference is believed to be because the former is the relaxed calculation and the latter is unrelaxed).
So it seems that the numerical result calculated by using “matrix_perturb” is neither the total dipole moment nor the electronic dipole moment…
After thinking it over, this discrepancy makes perfect sense.
Because CCSD (and MP2) use Hartree-Fock orbitals, we can easily separate the Lagrangian into a Hartree-Fock energy functional and a correlation energy functional. To get the orbital unrelaxed dipole moment, you want to compute the derivative of the Lagrangian energy with respect to an external electric field, except removing the orbital relaxation term in the correlation energy functional. You tried to accomplish this by eliminating all dependence of the Hartree-Fock orbitals and energies on the electric field. Unfortunately, this eliminates the dependence of the Hartree-Fock energy functional on your electric field, so you have completely eliminated the Hartree-Fock electronic dipole contribution. Unlike orbital relaxation, that is a very large contribution and goes far beyond the orbital unrelaxed approximation. As evidence, if you look at the breakdown of the final energy at the end of your unrelaxed output files, you see the final HF energy matches the final HF energy you computed with no perturbation.
You are not numerically differentiating the quantity you think you are. I can help with using Psi4, but devising the equations to perform this numerical differentiation for you is past what I have time to do.