1#!/usr/bin/env python 2# 3# Author: Qiming Sun <osirpt.sun@gmail.com> 4# 5 6import tempfile 7import numpy 8import h5py 9from pyscf import gto, scf, ao2mo 10 11''' 12Integral transformation for four different orbitals 13''' 14 15mol = gto.Mole() 16mol.build( 17 atom = [ 18 ["C", (-0.65830719, 0.61123287, -0.00800148)], 19 ["C", ( 0.73685281, 0.61123287, -0.00800148)], 20 ["C", ( 1.43439081, 1.81898387, -0.00800148)], 21 ["C", ( 0.73673681, 3.02749287, -0.00920048)], 22 ["C", (-0.65808819, 3.02741487, -0.00967948)], 23 ["C", (-1.35568919, 1.81920887, -0.00868348)], 24 ["H", (-1.20806619, -0.34108413, -0.00755148)], 25 ["H", ( 1.28636081, -0.34128013, -0.00668648)], 26 ["H", ( 2.53407081, 1.81906387, -0.00736748)], 27 ["H", ( 1.28693681, 3.97963587, -0.00925948)], 28 ["H", (-1.20821019, 3.97969587, -0.01063248)], 29 ["H", (-2.45529319, 1.81939187, -0.00886348)],], 30 31 basis = 'ccpvtz' 32) 33 34mf = scf.RHF(mol) 35mf.conv_tol = 1e-8 36e = mf.kernel() 37print('E = %.15g, ref -230.776765415' % e) 38 39# 40# Given four MOs, compute the MO-integrals and saved in dataset "mp2_bz" 41# 42eritmp = tempfile.NamedTemporaryFile() 43nocc = mol.nelectron // 2 44nvir = len(mf.mo_energy) - nocc 45co = mf.mo_coeff[:,:nocc] 46cv = mf.mo_coeff[:,nocc:] 47orbs = (co, cv, co, cv) 48# Depending on your hardware and BLAS library, it needs about 1 min on I5 3GHz 49# CPU with MKL library to transform the integrals 50ao2mo.general(mol, orbs, eritmp.name, dataname='mp2_bz')#, verbose=5) 51 52eia = mf.mo_energy[:nocc,None] - mf.mo_energy[None,nocc:] 53f = h5py.File(eritmp.name, 'r') 54eri = f['mp2_bz'] 55print('Note the shape of the transformed integrals (ij|kl) is %s.' % str(eri.shape)) 56print("It's a 2D array: the first index for compressed ij, the second index for compressed kl") 57 58emp2 = 0 59for i in range(nocc): 60 dajb = eia[i].reshape(-1,1) + eia.reshape(1,-1) 61 gi = numpy.array(eri[i*nvir:(i+1)*nvir]) 62 t2 = gi.flatten() / dajb.flatten() 63 gi = gi.reshape(nvir,nocc,nvir) 64 theta = gi*2 - gi.transpose(2,1,0) 65 emp2 += numpy.dot(t2, theta.flatten()) 66 67print('E_MP2 = %.15g, ref = -1.0435476768' % emp2) 68f.close() 69