#This is the script to analyse our calculations. It uses the lazy module of the pyhande library.
#It uses the pyhande/complex_analysis branch, git hash 6a5fdc7a4d0b80cb363e90afacb908cacfcdc12e
import pyhande
import collections
calcs = [["k111/ccsd.out"],["k211/ccsd.out"],["k221/ccsd.out","k221/ccsd_2.out"], ["k222/ccsd.out", "k222/ccsd_2.out", "k222/ccsd_3.out", "k222/ccsd_4.out", "k222/ccsd_5.out"], ["k111/ccsdt/ccsdt.out"],["k211/ccsdt/ccsdt.out", "k211/ccsdt/ccsdt_2.out", "k211/ccsdt/ccsdt_3.out", "k211/ccsdt/ccsdt_4.out", "k211/ccsdt/ccsdt_5.out"], ["k111/ccsdtq/ccsdtq.out"], ["k111/fciqmc/fciqmc.5759735.sdb.out"]]
results = collections.OrderedDict({})
fout = open("analysis.txt","w")
for files in calcs:
    info = pyhande.lazy.std_analysis(files, verbosity = -1)
    eref = info[0].metadata['reference']['H00']
    nkpoint = info[0].metadata['system']['nel']/8
    print >> fout, files
    print >> fout, "Starting iteration: ", info[0].metadata['pyhande']['reblock_start']
    print >> fout, "Reference energy: ", eref
    print >> fout, info[0].opt_block
    
    ecorr = info[0].opt_block['mean']['Proj. Energy']
    err = info[0].opt_block['standard error']['Proj. Energy']
    nex = info[0].metadata['reference']['ex_level']
    if not nkpoint in results:
        hf = open(files[0].split('/')[0] + '/HF.out','r')
        for line in hf:
            if "converged" in line:
                ehf = float(line.split()[4])
        hf.close()
        cc = open(files[0].split('/')[0] + '/CC_pyscf.out','r')
        for line in cc:
            if "E(CCSD)" in line:
                ecc = float(line.split()[5])
        cc.close()
        results[nkpoint] = [round(eref/nkpoint,8), '', '', '', '', '', '', '', '', round(ehf,8), round(ecc,8)]
    
    #Set fciqmc result to nex=5.
    nex = min(nex,5)
    results[nkpoint][2*nex-3] = round(ecorr/nkpoint,8)
    results[nkpoint][2*nex-2] = round(err/nkpoint,8)

dout = open("energies.dat","w")
dout.write("#kpoint nmr, Ref. E, Ecorr SD, std err, Ecorr SDT, std err, Ecorr SDTQ, std err, Ecorr FCI, std err, HF E, Ecorr SD (PySCF)\n")
for i, l in results.items():
    print >> dout, i, ', ', ', '.join(map(str,l))
