1#!/usr/bin/env python
2import numpy
3from pyscf import gto
4from pyscf import scf
5from pyscf import mcscf
6
7
8def run(b, mo0=None, dm0=None):
9    mol = gto.Mole()
10    mol.build(
11        verbose = 5,
12        output = 'o2rhf-%3.2f.out' % b,
13        atom = [
14            ['O', (0, 0,  b/2)],
15            ['O', (0, 0, -b/2)],],
16        basis = 'cc-pvdz',
17        spin = 2,
18        symmetry = 1,
19    )
20
21    mf = scf.RHF(mol)
22    mf.scf(dm0)
23
24    mc = mcscf.CASSCF(mf, 12, 8)
25    if mo0 is not None:
26        #from pyscf import lo
27        #mo0 = lo.orth.vec_lowdin(mo0, mf.get_ovlp())
28        mo0 = mcscf.project_init_guess(mc, mo0)
29    else:
30        mo0 = mcscf.sort_mo(mc, mf.mo_coeff, [5,6,7,8,9,11,12,13,14,15,16,17])
31        mc.max_orb_stepsize = .02
32    mc.kernel(mo0)
33    mc.analyze()
34    return mf, mc
35
36
37def urun(b, mo0=None, dm0=None):
38    mol = gto.Mole()
39    mol.build(
40        verbose = 5,
41        output = 'o2uhf-%3.2f.out' % b,
42        atom = [
43            ['O', (0, 0,  b/2)],
44            ['O', (0, 0, -b/2)],],
45        basis = 'cc-pvdz',
46        spin = 2,
47    )
48
49    mf = scf.UHF(mol)
50    mf.scf(dm0)
51
52    mc = mcscf.CASSCF(mf, 12, 8)
53    if mo0 is not None:
54        #from pyscf import lo
55        #mo0 =(lo.orth.vec_lowdin(mo0[0], mf.get_ovlp()),
56        #      lo.orth.vec_lowdin(mo0[1], mf.get_ovlp()))
57        mo0 = mcscf.project_init_guess(mc, mo0)
58    mc.kernel(mo0)
59    mc.analyze()
60    return mf, mc
61
62x = numpy.hstack((numpy.arange(0.9, 2.01, 0.1),
63                  numpy.arange(2.1, 4.01, 0.1)))
64
65dm0 = mo0 = None
66eumc = []
67euhf = []
68s = []
69for b in reversed(x):
70    mf, mc = urun(b, mo0, dm0)
71    mo0 = mc.mo_coeff
72    dm0 = mf.make_rdm1()
73    s.append(mc.spin_square()[1])
74    euhf.append(mf.hf_energy)
75    eumc.append(mc.e_tot)
76
77euhf.reverse()
78eumc.reverse()
79s.reverse()
80#print s
81
82dm0 = mo0 = None
83ermc = []
84erhf = []
85for b in x:
86    mf, mc = run(b, mo0, dm0)
87    mo0 = mc.mo_coeff
88    dm0 = mf.make_rdm1()
89    erhf.append(mf.hf_energy)
90    ermc.append(mc.e_tot)
91
92with open('o2-scan.txt', 'w') as fout:
93    fout.write('  ROHF 0.9->4.0   RCAS(12,8)    UHF 4.0->0.9  UCAS(12,8) \n')
94    for i, xi in enumerate(x):
95        fout.write('%2.1f  %12.8f  %12.8f  %12.8f  %12.8f\n'
96                   % (xi, erhf[i], ermc[i], euhf[i], eumc[i]))
97
98import matplotlib.pyplot as plt
99plt.plot(x, erhf, label='ROHF,0.9->4.0')
100plt.plot(x, euhf, label='UHF, 4.0->0.9')
101plt.plot(x, ermc, label='RCAS(6,6),0.9->4.0')
102plt.plot(x, eumc, label='UCAS(6,6),4.0->0.9')
103plt.legend()
104plt.show()
105