detopt_data/plot_bo.py at master · cqsl/detopt_data · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
from plot_common import *
prefix = 'plots/'
matplotlib.use('Agg')

plt.figure(figsize=(8.8/2.54,6.6/2.54))

plt.subplots_adjust(
    left=0.125,
    right=0.99,
    bottom=0.15,
    top=0.98
)


molecules = ['c2.xyz', 'bn.xyz', 'beo.xyz', 'lif.xyz', 'bh3.xyz', 'n2-2-118.xyz']
files_bo = []
for m in molecules:
    files_bo = files_bo + list(glob.glob(f'output/normal/{m.replace('.xyz', '')}_*.h5'))

files_bo = list(sorted(files_bo))
data_bo = load(files_bo)


markers = ['v', 'o', 's', '^', 'x', '+']
basis = 'ccpvdz'

mit = iter(markers)
cit = iter(colors[1:])

ndets = [1,2,4,8,16,32,64,128,256,512,768]

refs = [e_fci_gao, e_fci_psi4, e_dmrg_chan, e_dmrg_block2]
refs_is_fci = [True, True, False, False]
spin=0
df_rows = ndets + ['FCI', 'DMRG']
df_cols = [molecule_dict.get(m, m).replace('$\\mathregular{', '').replace('$', '').replace('}', '').replace('^{\\star', '') for m in molecules]
df_data = np.full((len(df_rows), len(df_cols)), None)

for i, molecule in enumerate(molecules):
    col = next(cit)
    mar = next(mit)
    E = []
    ndet = []
    for k, v in data_bo.items():
        if k[0] != molecule or k[1] != basis or k[3] not in ndets:
            continue
        E.append(v['e'].real)
        ndet.append(k[3])
    assert len(E) > 0
    E = np.array(E)
    ndet = np.array(ndet)
    perm = np.argsort(ndet)
    E = E[perm]
    ndet = ndet[perm]

    for ref, is_fci in zip(refs, refs_is_fci):
        eref = ref.get(basis, {}).get(molecule, None)
        if eref is not None:

            if is_fci:
                df_data[-2,i] = eref
            else:
                df_data[-1,i] = eref

            break

    err = E - eref
    molname = molecule_dict.get(molecule, molecule)
    nee = sum(data_bo[molecule, basis,spin, ndet[-1], None, None]['nelec'])
    noo = data_bo[molecule, basis, spin, ndet[-1], None, None]['nao']
    lll = f"{molname} ({noo}, {nee})"
    l, = plt.plot(1/np.array(ndet), err, label=lll,
                  marker=mar, markerfacecolor=matplotlib.colors.to_rgb(col)+(0.75,),
                  linewidth=0.,
                 )

    for ndeti, ei in zip(ndet, E):
        l, = np.where(ndeti == np.array(ndets))
        df_data[l,i] = ei.real

plt.ylabel(r'$\mathregular{E-E_{ref}\quad[a.u.]}$')


# move N2 to first pos
handles, labels = plt.gca().get_legend_handles_labels()
handles = [handles[-1], *handles[:-1]]
labels = [labels[-1], *labels[:-1]]
plt.legend(handles, labels, title=r'MOLECULE ($\mathregular{m, n}$)', loc=4, frameon=False)

plt.yscale('log')
plt.xscale('log', base=10)
plt.xlabel(r'$\mathregular{N_D^{-1}}$')
plt.xlim(7e-4, 2)
plt.ylim(2e-5, 6e-1)


plt.savefig(prefix+'fig2.pdf')


if True:
    import pandas as pd
    from pathlib import Path
    df = pd.DataFrame(df_data, df_rows, df_cols)
    df.index.name = r"N_D \ molecule"
    path = Path(prefix+"source_data.xlsx")
    with pd.ExcelWriter(path, mode="a" if path.exists() else "w", engine="openpyxl", if_sheet_exists="replace" if path.exists() else None) as writer:
        df.to_excel(writer, sheet_name="fig2", index=True)