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