1#!/usr/bin/env python 2# Copyright 2014-2018 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15# 16 17import unittest 18import numpy 19from pyscf import lib 20from pyscf.pbc import gto 21from pyscf.pbc import scf 22from pyscf.pbc.cc import kccsd_uhf 23from pyscf.pbc.cc import kccsd 24from pyscf.pbc.cc import eom_kccsd_ghf 25from pyscf.pbc.cc import eom_kccsd_uhf 26from pyscf.pbc.lib import kpts_helper 27 28cell = gto.Cell() 29cell.atom=''' 30He 0.000000000000 0.000000000000 0.000000000000 31He 1.685068664391 1.685068664391 1.685068664391 32''' 33cell.basis = [[0, (1., 1.)], [0, (.5, 1.)]] 34cell.a = ''' 350.000000000, 3.370137329, 3.370137329 363.370137329, 0.000000000, 3.370137329 373.370137329, 3.370137329, 0.000000000''' 38cell.unit = 'B' 39cell.mesh = [13,13,13] 40cell.build() 41 42 43def rand_t1_t2(cell, kpts, nocc, nvir): 44 nkpts = len(kpts) 45 nocca, noccb = nocc 46 nvira, nvirb = nvir 47 t1a = (numpy.random.random((nkpts,nocca,nvira)) + 48 numpy.random.random((nkpts,nocca,nvira))*1j - .5-.5j) 49 t1b = (numpy.random.random((nkpts,noccb,nvirb)) + 50 numpy.random.random((nkpts,noccb,nvirb))*1j - .5-.5j) 51 t2aa = (numpy.random.random((nkpts,nkpts,nkpts,nocca,nocca,nvira,nvira)) + 52 numpy.random.random((nkpts,nkpts,nkpts,nocca,nocca,nvira,nvira))*1j - .5-.5j) 53 kconserv = kpts_helper.get_kconserv(cell, kpts) 54 t2aa = t2aa - t2aa.transpose(1,0,2,4,3,5,6) 55 tmp = t2aa.copy() 56 for ki, kj, kk in kpts_helper.loop_kkk(nkpts): 57 kl = kconserv[ki, kk, kj] 58 t2aa[ki,kj,kk] = t2aa[ki,kj,kk] - tmp[ki,kj,kl].transpose(0,1,3,2) 59 t2ab = (numpy.random.random((nkpts,nkpts,nkpts,nocca,noccb,nvira,nvirb)) + 60 numpy.random.random((nkpts,nkpts,nkpts,nocca,noccb,nvira,nvirb))*1j - .5-.5j) 61 t2bb = (numpy.random.random((nkpts,nkpts,nkpts,noccb,noccb,nvirb,nvirb)) + 62 numpy.random.random((nkpts,nkpts,nkpts,noccb,noccb,nvirb,nvirb))*1j - .5-.5j) 63 t2bb = t2bb - t2bb.transpose(1,0,2,4,3,5,6) 64 tmp = t2bb.copy() 65 for ki, kj, kk in kpts_helper.loop_kkk(nkpts): 66 kl = kconserv[ki, kk, kj] 67 t2bb[ki,kj,kk] = t2bb[ki,kj,kk] - tmp[ki,kj,kl].transpose(0,1,3,2) 68 69 t1 = (t1a, t1b) 70 t2 = (t2aa, t2ab, t2bb) 71 return t1, t2 72 73class KnownValues(unittest.TestCase): 74 def test_amplitudes_to_vector(self): 75 numpy.random.seed(1) 76 kpts = cell.make_kpts([1,2,3]) 77 nmo = (8, 8) 78 nocc = (3, 2) 79 nvir = (nmo[0]-nocc[0], nmo[1]-nocc[1]) 80 t1, t2 = rand_t1_t2(cell, kpts, nocc, nvir) 81 vec = kccsd_uhf.amplitudes_to_vector(t1, t2) 82 r1, r2 = kccsd_uhf.vector_to_amplitudes(vec, nmo, nocc, len(kpts)) 83 self.assertAlmostEqual(abs(t1[0]-r1[0]).max(), 0, 12) 84 self.assertAlmostEqual(abs(t1[1]-r1[1]).max(), 0, 12) 85 self.assertAlmostEqual(abs(t2[0]-r2[0]).max(), 0, 12) 86 self.assertAlmostEqual(abs(t2[1]-r2[1]).max(), 0, 12) 87 self.assertAlmostEqual(abs(t2[2]-r2[2]).max(), 0, 12) 88 89 vec1 = kccsd_uhf.amplitudes_to_vector(r1, r2) 90 self.assertAlmostEqual(abs(vec-vec1).max(), 0, 12) 91 92 def test_spatial2spin(self): 93 numpy.random.seed(1) 94 kpts = cell.make_kpts([1,2,3]) 95 nmo = (8, 8) 96 nocc = (3, 2) 97 nvir = (nmo[0]-nocc[0], nmo[1]-nocc[1]) 98 t1, t2 = rand_t1_t2(cell, kpts, nocc, nvir) 99 100 orbspin = numpy.zeros((len(kpts),nmo[0]+nmo[1]), dtype=int) 101 orbspin[:,1::2] = 1 102 kconserv = kpts_helper.get_kconserv(cell, kpts) 103 104 r1 = kccsd.spatial2spin(t1, orbspin, kconserv) 105 r1 = kccsd.spin2spatial(r1, orbspin, kconserv) 106 self.assertAlmostEqual(abs(r1[0]-t1[0]).max(), 0, 12) 107 self.assertAlmostEqual(abs(r1[1]-t1[1]).max(), 0, 12) 108 109 r2 = kccsd.spatial2spin(t2, orbspin, kconserv) 110 r2 = kccsd.spin2spatial(r2, orbspin, kconserv) 111 self.assertAlmostEqual(abs(r2[0]-t2[0]).max(), 0, 12) 112 self.assertAlmostEqual(abs(r2[1]-t2[1]).max(), 0, 12) 113 self.assertAlmostEqual(abs(r2[2]-t2[2]).max(), 0, 12) 114 115 def test_eris(self): 116 numpy.random.seed(1) 117 kpts = cell.make_kpts([1,1,3]) 118 kmf = scf.KUHF(cell, kpts=kpts, exxdiv=None) 119 nmo = cell.nao_nr() 120 kmf.mo_occ = numpy.zeros((2,3,nmo)) 121 kmf.mo_occ[0,:,:3] = 1 122 kmf.mo_occ[1,:,:1] = 1 123 kmf.mo_energy = (numpy.arange(nmo) + 124 numpy.random.random((2,3,nmo)) * .3) 125 kmf.mo_energy[kmf.mo_occ == 0] += 2 126 127 kmf.mo_coeff = (numpy.random.random((2,3,nmo,nmo)) + 128 numpy.random.random((2,3,nmo,nmo))*1j - .5-.5j) 129 130 mycc = kccsd_uhf.UCCSD(kmf) 131 eris = mycc.ao2mo() 132 133 self.assertAlmostEqual(lib.finger(eris.fock[0]), 0.53719738596848132 +0.83031462049142502j, 11) 134 self.assertAlmostEqual(lib.finger(eris.fock[1]), 0.043623927844025398+0.20815796288019522j, 11) 135 136 self.assertAlmostEqual(lib.finger(eris.oooo), 0.10126616996580763 -0.89787173918481156j , 11) 137 self.assertAlmostEqual(lib.finger(eris.ooov), -0.035814310241888067 +0.20393025075274804j , 11) 138 self.assertAlmostEqual(lib.finger(eris.oovv), -0.032378345849800663 +0.015060519910448879j , 11) 139 self.assertAlmostEqual(lib.finger(eris.ovov), 0.017919215232962762 -0.37180556037878837j , 11) 140 self.assertAlmostEqual(lib.finger(eris.voov), -0.33038865500581482 +0.18384096784449852j , 11) 141 self.assertAlmostEqual(lib.finger(eris.vovv), 0.078104278754223946 +0.0004014143354997223j, 11) 142 self.assertAlmostEqual(lib.finger(eris.vvvv), -0.0199910973368542 -0.0019864189992825137j, 11) 143 144 self.assertAlmostEqual(lib.finger(eris.OOOO), 0.022687859086726745 +0.0076542690105189095j, 11) 145 self.assertAlmostEqual(lib.finger(eris.OOOV), -0.024119030579269278 -0.15249100640417029j , 11) 146 self.assertAlmostEqual(lib.finger(eris.OOVV), 0.085942751462484687 -0.27088394382044989j , 11) 147 self.assertAlmostEqual(lib.finger(eris.OVOV), 0.35291244981540776 +0.080119865729794376j , 11) 148 self.assertAlmostEqual(lib.finger(eris.VOOV), 0.0045484393536995267 +0.0094123990059577414j, 11) 149 self.assertAlmostEqual(lib.finger(eris.VOVV), -0.28341581692294759 +0.0022174023470048921j, 11) 150 self.assertAlmostEqual(lib.finger(eris.VVVV), 0.96007536729340814 -0.019410945571596398j , 11) 151 152 self.assertAlmostEqual(lib.finger(eris.ooOO), -0.32831508979976765 -0.32180378432620471j , 11) 153 self.assertAlmostEqual(lib.finger(eris.ooOV), 0.33617152217704632 -0.34130555912360216j , 11) 154 self.assertAlmostEqual(lib.finger(eris.ooVV), -0.00011230004797088339-1.2850251519380604j , 11) 155 self.assertAlmostEqual(lib.finger(eris.ovOV), 0.1365118156144336 +0.16999367231786541j , 11) 156 self.assertAlmostEqual(lib.finger(eris.voOV), 0.19736044623396901 -0.047060848969879901j , 11) 157 self.assertAlmostEqual(lib.finger(eris.voVV), 0.44513499758124858 +0.06343379901453805j , 11) 158 self.assertAlmostEqual(lib.finger(eris.vvVV), -0.070971875998391304 -0.31253893124900545j , 11) 159 160 #self.assertAlmostEqual(lib.finger(eris.OOoo), 0.031140414688898856 -0.23913617484062258j , 11) 161 self.assertAlmostEqual(lib.finger(eris.OOov), 0.20355552926191381 +0.18712171841650935j , 11) 162 self.assertAlmostEqual(lib.finger(eris.OOvv), 0.070789122903945706 -0.013360818695166678j , 11) 163 #self.assertAlmostEqual(lib.finger(eris.OVov), 0.38230103404493437 -0.019845885264560274j , 11) 164 #self.assertAlmostEqual(lib.finger(eris.VOov), 0.081760186267865437 -0.052409714443657308j , 11) 165 self.assertAlmostEqual(lib.finger(eris.VOvv), -0.036061642075282056 +0.019284185856121634j , 11) 166 #self.assertAlmostEqual(lib.finger(eris.VVvv), 0.13458896578260207 -0.11322854172459119j , 11) 167 168 def test_spatial2spin_ip(self): 169 numpy.random.seed(1) 170 kpts = cell.make_kpts([1,2,3]) 171 nkpts = len(kpts) 172 nmo = (8, 8) 173 nocc = (3, 2) 174 nvir = (nmo[0]-nocc[0], nmo[1]-nocc[1]) 175 176 orbspin = numpy.zeros((len(kpts),nmo[0]+nmo[1]), dtype=int) 177 orbspin[:,1::2] = 1 178 kconserv = kpts_helper.get_kconserv(cell, kpts) 179 180 kshift = 0 181 spin_r1_ip = (numpy.random.rand(nocc[0]+nocc[1])*1j + 182 numpy.random.rand(nocc[0]+nocc[1]) - 0.5 - 0.5*1j) 183 spin_r2_ip = (numpy.random.rand(nkpts**2 * (nocc[0]+nocc[1])**2 * (nvir[0]+nvir[1])) + 184 numpy.random.rand(nkpts**2 * (nocc[0]+nocc[1])**2 * (nvir[0]+nvir[1]))*1j - 185 0.5 - 0.5*1j) 186 spin_r2_ip = spin_r2_ip.reshape(nkpts, nkpts, (nocc[0]+nocc[1]), 187 (nocc[0]+nocc[1]), (nvir[0]+nvir[1])) 188 spin_r2_ip = eom_kccsd_ghf.enforce_2p_spin_ip_doublet(spin_r2_ip, kconserv, kshift, orbspin) 189 190 [r1a, r1b], [r2aaa, r2baa, r2abb, r2bbb] = \ 191 eom_kccsd_ghf.spin2spatial_ip_doublet(spin_r1_ip, spin_r2_ip, kconserv, kshift, orbspin) 192 193 r1, r2 = eom_kccsd_ghf.spatial2spin_ip_doublet([r1a, r1b], [r2aaa, r2baa, r2abb, r2bbb], 194 kconserv, kshift, orbspin=orbspin) 195 196 self.assertAlmostEqual(abs(r1-spin_r1_ip).max(), 0, 12) 197 self.assertAlmostEqual(abs(r2-spin_r2_ip).max(), 0, 12) 198 199 def test_spatial2spin_ea(self): 200 numpy.random.seed(1) 201 kpts = cell.make_kpts([1,2,3]) 202 nkpts = len(kpts) 203 nmo = (8, 8) 204 nocc = (3, 2) 205 nvir = (nmo[0]-nocc[0], nmo[1]-nocc[1]) 206 207 orbspin = numpy.zeros((len(kpts),nmo[0]+nmo[1]), dtype=int) 208 orbspin[:,1::2] = 1 209 kconserv = kpts_helper.get_kconserv(cell, kpts) 210 211 kshift = 0 212 spin_r1_ea = (numpy.random.rand(nvir[0]+nvir[1])*1j + 213 numpy.random.rand(nvir[0]+nvir[1]) - 0.5 - 0.5*1j) 214 spin_r2_ea = (numpy.random.rand(nkpts**2 * (nocc[0]+nocc[1])* (nvir[0]+nvir[1])**2) + 215 numpy.random.rand(nkpts**2 * (nocc[0]+nocc[1])* (nvir[0]+nvir[1])**2)*1j - 216 0.5 - 0.5*1j) 217 spin_r2_ea = spin_r2_ea.reshape(nkpts, nkpts, (nocc[0]+nocc[1]), 218 (nvir[0]+nvir[1]), (nvir[0]+nvir[1])) 219 spin_r2_ea = eom_kccsd_ghf.enforce_2p_spin_ea_doublet(spin_r2_ea, kconserv, kshift, orbspin) 220 221 [r1a, r1b], [r2aaa, r2baa, r2abb, r2bbb] = \ 222 eom_kccsd_ghf.spin2spatial_ea_doublet(spin_r1_ea, spin_r2_ea, kconserv, kshift, orbspin) 223 224 r1, r2 = eom_kccsd_ghf.spatial2spin_ea_doublet([r1a, r1b], [r2aaa, r2baa, r2abb, r2bbb], 225 kconserv, kshift, orbspin=orbspin) 226 227 self.assertAlmostEqual(abs(r1-spin_r1_ea).max(), 0, 12) 228 self.assertAlmostEqual(abs(r2-spin_r2_ea).max(), 0, 12) 229 230if __name__ == '__main__': 231 print("KUCCSD tests") 232 unittest.main() 233