jjgoings / McMurchie-Davidson

do a simple closed shell Hartree-Fock using McMurchie-Davidson to compute integrals
BSD 3-Clause "New" or "Revised" License
78 stars 17 forks source link

Do you have any plans to put in the capability of CCSD(T) some time? #4

Open jl2922 opened 5 years ago

jjgoings commented 5 years ago

No, I had not planned on implementing CCSD(T) at the moment. It certainly could be done by adding a new routine in postscf.py. I doubt it would be very fast, but the infrastructure is definitely there.

I'll look into it, but if you'd like to code it yourself and submit a pull request, I'm happy to take contributions!

jl2922 commented 5 years ago

I am not familiar with the detailed algorithm of coupled cluster, but I can try...

jl2922 commented 5 years ago

I followed your blog post on CCSD (http://joshuagoings.com/2013/07/17/coupled-cluster-with-singles-and-doubles-ccsd-in-python/), but the energy I got is slightly different from Molpro and PySCF by several microHartree, even though the convergence criterion is 1 nanoHartree. Did you make any approximations there? Molpro and PySCF agree all the way to the convergence tolerance for HeH+.

jjgoings commented 5 years ago

Check that the SCF and correlation energy match. I found that the nuclear repulsion energy was wrong (go figure) so I fixed that on the site. That explains most of the discrepancy. It's possible the "hard-coded" integrals do not agree 100% for various numerical precision issues.

jl2922 commented 5 years ago

I noticed that part and I thought maybe you were using a different convention. But, besides the nuclear repulsion, the CCSD correction itself is different, as mentioned in the previous message. I have merged you CCSD code to this repo, so the integrals are not hard coded.

jjgoings commented 5 years ago

I will need to see the code and some tests then. The case of HeH+ works for the values hard-coded, but unfortunately I cannot test it on anything else since it is not generalized, and I do not know how you interfaced the codes. How different is the SCF energy? Is the SCF different by a micro Hartree? A micro Hartree is small, and there could be plenty of reasons why they do not match exactly (e.g. integral screening, SCF algorithm, etc.). That's an error of 0.0006 kcal/mol. I'd be more concerned if you are off by milli Hartree.

jl2922 commented 5 years ago

The code is here: https://github.com/jl2922/McMurchie-Davidson/blob/ccsd/mmd/postscf.py

The SCF energy is exactly the same as Molpro. Only the CCSD correction is different. In this case, what might be the reason?

jjgoings commented 5 years ago

Okay, I ran the code. I agree something is not right, but I need some time to figure out where the error lies. It may be in the transition to spin orbital basis or it may be in the routines themselves.

Update: the initial guess T2 amplitudes give the correct MP2 correction, so the error must be in the CCSD iterations

erikkjellgren commented 5 years ago

Hi,

I think I some time ago noted a small mistake in the coding. Completely forgot to report, sorry. I'll see if I can track it again if I get time in the weekend. In the mean time, if a reference is needed to track down the code, you can try paste the python code from https://github.com/Melisius/SlowQuant/blob/master/slowquant/coupledcluster/CythonCC.pyx which should work.

jjgoings commented 5 years ago

I managed to track down the bugs in that blog post. Thanks @jl2922 for pointing them out. The error was in a couple permutations for constructing the T2. One messed up sign, and two cases of forgetting to permute indices. The code on the blog works for HeH+, as well as STO-3G water in the cases.

Anyway, I'll get to thinking about the best way of implementing CCSD(T). The Psi4NumPy have nice implementation that could basically be ported into the code as-is (https://github.com/psi4/psi4numpy)

In fact, their code it probably more functional than mine since they have access to Psi4 for all the integral evaluations, etc. Might be worth checking out if you haven't already.

I'll leave this issue open for now, and when I have time I'll consider making a cleaner/clearer implementation of CCSD(T).

jl2922 commented 5 years ago

Thank you! These are quite helpful.