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