1#!/usr/bin/env python
2# Copyright 2014-2021 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 copy
19import numpy
20from pyscf import gto
21from pyscf import scf
22from pyscf import lib
23from pyscf.mp.dfump2_native import DFUMP2, SCSUMP2
24
25
26mol = gto.Mole()
27mol.verbose = 0
28mol.output = None
29mol.atom = '''
30O    0.000   0.000  -1.141
31O    0.000   0.000   1.141
32'''
33mol.unit = 'Bohr'
34mol.basis = 'def2-SVP'
35mol.spin = 2
36mol.build()
37
38mf = scf.UHF(mol)
39mf.conv_tol = 1.0e-12
40mf.kernel()
41
42
43def check_orth(obj, mol, mo_coeff, thresh=1.0e-12):
44    sao = mol.intor_symmetric('int1e_ovlp')
45    sno = numpy.linalg.multi_dot([mo_coeff.T, sao, mo_coeff])
46    I = numpy.eye(mo_coeff.shape[1])
47    obj.assertTrue(numpy.allclose(sno, I, atol=1.0e-12))
48
49
50class KnownValues(unittest.TestCase):
51
52    def setUp(self):
53        self.assertTrue(mf.converged)
54        self.mf = copy.copy(mf)
55        self.mf.mol = mf.mol.copy()
56        self.mf.mo_coeff = mf.mo_coeff.copy()
57        self.mf.mo_occ = mf.mo_occ.copy()
58        self.mf.mo_energy = mf.mo_energy.copy()
59
60    def test_energy(self):
61        with DFUMP2(self.mf) as pt:
62            pt.kernel()
63            self.assertAlmostEqual(pt.e_corr, -0.347887316046, delta=1.0e-8)
64            self.assertAlmostEqual(pt.e_tot, -149.838271530605, delta=1.0e-8)
65
66    def test_energy_fc(self):
67        with DFUMP2(self.mf, frozen=2) as pt:
68            pt.kernel()
69            self.assertAlmostEqual(pt.e_corr, -0.343816318675, delta=1.0e-8)
70            self.assertAlmostEqual(pt.e_tot, -149.834200533235, delta=1.0e-8)
71
72    def test_energy_fclist(self):
73        self.mf.mo_coeff[0, :, [1, 3]] = self.mf.mo_coeff[0, :, [3, 1]]
74        self.mf.mo_energy[0, [1, 3]] = self.mf.mo_energy[0, [3, 1]]
75        self.mf.mo_coeff[1, :, [0, 4]] = self.mf.mo_coeff[1, :, [4, 0]]
76        self.mf.mo_energy[1, [0, 4]] = self.mf.mo_energy[1, [4, 0]]
77        with DFUMP2(self.mf, frozen=[[0, 3], [1, 4]]) as pt:
78            pt.kernel()
79            self.assertAlmostEqual(pt.e_corr, -0.343816318675, delta=1.0e-8)
80            self.assertAlmostEqual(pt.e_tot, -149.834200533235, delta=1.0e-8)
81
82    def test_natorbs(self):
83        mol = self.mf.mol
84        with DFUMP2(self.mf) as pt:
85            natocc, natorb = pt.make_natorbs()
86            # number of electrons conserved
87            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
88            # orbitals orthogonal
89            check_orth(self, mol, natorb)
90            # selected values
91            self.assertAlmostEqual(natocc[0], 1.9999191951, delta=1.0e-7)
92            self.assertAlmostEqual(natocc[6], 1.9473283296, delta=1.0e-7)
93            self.assertAlmostEqual(natocc[7], 1.0168954406, delta=1.0e-7)
94            self.assertAlmostEqual(natocc[8], 1.0168954406, delta=1.0e-7)
95            self.assertAlmostEqual(natocc[9], 0.0478262909, delta=1.0e-7)
96            self.assertAlmostEqual(natocc[27], 0.0002326288, delta=1.0e-7)
97
98    def test_natorbs_fc(self):
99        mol = self.mf.mol
100        with DFUMP2(self.mf, frozen=2) as pt:
101            natocc, natorb = pt.make_natorbs()
102            # number of electrons conserved
103            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
104            # orbital orthogonal
105            check_orth(self, mol, natorb)
106            # selected values
107            self.assertAlmostEqual(natocc[0], 1.9999987581, delta=1.0e-7)
108            self.assertAlmostEqual(natocc[1], 1.9999987356, delta=1.0e-7)
109            self.assertAlmostEqual(natocc[2], 1.9882629065, delta=1.0e-7)
110            self.assertAlmostEqual(natocc[6], 1.9473585838, delta=1.0e-7)
111            self.assertAlmostEqual(natocc[7], 1.0168965649, delta=1.0e-7)
112            self.assertAlmostEqual(natocc[8], 1.0168965649, delta=1.0e-7)
113            self.assertAlmostEqual(natocc[9], 0.0477790944, delta=1.0e-7)
114            self.assertAlmostEqual(natocc[27], 0.0002307322, delta=1.0e-7)
115
116    def test_natorbs_fclist(self):
117        self.mf.mo_coeff[0, :, [1, 3]] = self.mf.mo_coeff[0, :, [3, 1]]
118        self.mf.mo_energy[0, [1, 3]] = self.mf.mo_energy[0, [3, 1]]
119        self.mf.mo_coeff[1, :, [0, 4]] = self.mf.mo_coeff[1, :, [4, 0]]
120        self.mf.mo_energy[1, [0, 4]] = self.mf.mo_energy[1, [4, 0]]
121        with DFUMP2(self.mf, frozen=[[0, 3], [1, 4]]) as pt:
122            # also check the density matrix
123            rdm1 = pt.make_rdm1()
124            self.assertAlmostEqual(rdm1[0, 0, 0], 1.0, delta=1.0e-12)
125            self.assertAlmostEqual(rdm1[0, 3, 3], 1.0, delta=1.0e-12)
126            self.assertAlmostEqual(rdm1[1, 1, 1], 1.0, delta=1.0e-12)
127            self.assertAlmostEqual(rdm1[1, 4, 4], 1.0, delta=1.0e-12)
128            # now calculate the natural orbitals
129            natocc, natorb = pt.make_natorbs()
130            # number of electrons conserved
131            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
132            # orbital orthogonal
133            check_orth(self, mol, natorb)
134            # selected values
135            self.assertAlmostEqual(natocc[0], 1.9999987581, delta=1.0e-7)
136            self.assertAlmostEqual(natocc[1], 1.9999987356, delta=1.0e-7)
137            self.assertAlmostEqual(natocc[2], 1.9882629065, delta=1.0e-7)
138            self.assertAlmostEqual(natocc[6], 1.9473585838, delta=1.0e-7)
139            self.assertAlmostEqual(natocc[7], 1.0168965649, delta=1.0e-7)
140            self.assertAlmostEqual(natocc[8], 1.0168965649, delta=1.0e-7)
141            self.assertAlmostEqual(natocc[9], 0.0477790944, delta=1.0e-7)
142            self.assertAlmostEqual(natocc[27], 0.0002307322, delta=1.0e-7)
143
144    def test_natorbs_relaxed(self):
145        mol = self.mf.mol
146        with DFUMP2(self.mf) as pt:
147            pt.cphf_tol = 1e-12
148            natocc, natorb = pt.make_natorbs(relaxed=True)
149            # number of electrons conserved
150            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
151            # orbitals orthogonal
152            check_orth(self, mol, natorb)
153            # selected values
154            self.assertAlmostEqual(natocc[0], 1.9999198031, delta=1.0e-7)
155            self.assertAlmostEqual(natocc[6], 1.9478407509, delta=1.0e-7)
156            self.assertAlmostEqual(natocc[7], 1.0169668947, delta=1.0e-7)
157            self.assertAlmostEqual(natocc[8], 1.0169668947, delta=1.0e-7)
158            self.assertAlmostEqual(natocc[9], 0.0453923546, delta=1.0e-7)
159            self.assertAlmostEqual(natocc[27], 0.0002225494, delta=1.0e-7)
160
161    def test_natorbs_relaxed_fc(self):
162        mol = self.mf.mol
163        with DFUMP2(self.mf, frozen=2) as pt:
164            pt.cphf_tol = 1e-12
165            natocc, natorb = pt.make_natorbs(relaxed=True)
166            # number of electrons conserved
167            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
168            # orbital orthogonal
169            check_orth(self, mol, natorb)
170            # selected values
171            self.assertAlmostEqual(natocc[0], 2.0000050774, delta=1.0e-7)
172            self.assertAlmostEqual(natocc[1], 2.0000042352, delta=1.0e-7)
173            self.assertAlmostEqual(natocc[2], 1.9889171379, delta=1.0e-7)
174            self.assertAlmostEqual(natocc[6], 1.9478689720, delta=1.0e-7)
175            self.assertAlmostEqual(natocc[7], 1.0169674773, delta=1.0e-7)
176            self.assertAlmostEqual(natocc[8], 1.0169674773, delta=1.0e-7)
177            self.assertAlmostEqual(natocc[9], 0.0453427169, delta=1.0e-7)
178            self.assertAlmostEqual(natocc[27], 0.0002207476, delta=1.0e-7)
179
180    def test_natorbs_relaxed_fclist(self):
181        self.mf.mo_coeff[0, :, [0, 5]] = self.mf.mo_coeff[0, :, [5, 0]]
182        self.mf.mo_energy[0, [0, 5]] = self.mf.mo_energy[0, [5, 0]]
183        self.mf.mo_coeff[1, :, [1, 4]] = self.mf.mo_coeff[1, :, [4, 1]]
184        self.mf.mo_energy[1, [1, 4]] = self.mf.mo_energy[1, [4, 1]]
185        mol = self.mf.mol
186        with DFUMP2(self.mf, frozen=[[1, 5], [0, 4]]) as pt:
187            pt.cphf_tol = 1e-12
188            natocc, natorb = pt.make_natorbs(relaxed=True)
189            # number of electrons conserved
190            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
191            # orbital orthogonal
192            check_orth(self, mol, natorb)
193            # selected values
194            self.assertAlmostEqual(natocc[0], 2.0000050774, delta=1.0e-7)
195            self.assertAlmostEqual(natocc[1], 2.0000042352, delta=1.0e-7)
196            self.assertAlmostEqual(natocc[2], 1.9889171379, delta=1.0e-7)
197            self.assertAlmostEqual(natocc[6], 1.9478689720, delta=1.0e-7)
198            self.assertAlmostEqual(natocc[7], 1.0169674773, delta=1.0e-7)
199            self.assertAlmostEqual(natocc[8], 1.0169674773, delta=1.0e-7)
200            self.assertAlmostEqual(natocc[9], 0.0453427169, delta=1.0e-7)
201            self.assertAlmostEqual(natocc[27], 0.0002207476, delta=1.0e-7)
202
203    def test_memory(self):
204        # Dummy class to set a very low memory limit.
205        class fakeDFUMP2(DFUMP2):
206            _mem_kb = 0
207            @property
208            def max_memory(self):
209                return lib.current_memory()[0] + 1e-3 * self._mem_kb
210            @max_memory.setter
211            def max_memory(self, val):
212                pass
213        with fakeDFUMP2(self.mf) as pt:
214            E, natocc_ur, natocc_re = None, None, None
215            pt.cphf_tol = 1e-12
216            # Try very low amounts of memory (in kB) until there is no failure.
217            # Assume it should certainly work before 1 MB is reached.
218            for m in range(8, 1000, 8):
219                pt._mem_kb = m
220                try:
221                    E = pt.kernel()
222                except MemoryError:
223                    pass
224                else:
225                    break
226            for m in range(4, 1000, 4):
227                pt._mem_kb = m
228                try:
229                    natocc_ur = pt.make_natorbs()[0]
230                except MemoryError:
231                    pass
232                else:
233                    break
234            for m in range(20, 1000, 20):
235                pt._mem_kb = m
236                try:
237                    natocc_re = pt.make_natorbs(relaxed=True)[0]
238                except MemoryError:
239                    pass
240                else:
241                    break
242        self.assertAlmostEqual(E, -0.347887316046, delta=1.0e-8)
243        self.assertAlmostEqual(natocc_ur[6], 1.9473283296, delta=1.0e-7)
244        self.assertAlmostEqual(natocc_ur[7], 1.0168954406, delta=1.0e-7)
245        self.assertAlmostEqual(natocc_ur[9], 0.0478262909, delta=1.0e-7)
246        self.assertAlmostEqual(natocc_re[6], 1.9478407509, delta=1.0e-7)
247        self.assertAlmostEqual(natocc_re[7], 1.0169668947, delta=1.0e-7)
248        self.assertAlmostEqual(natocc_re[9], 0.0453923546, delta=1.0e-7)
249
250    def test_scs_energy(self):
251        with SCSUMP2(self.mf) as pt:
252            pt.kernel()
253            self.assertAlmostEqual(pt.e_corr, -0.324631353397, delta=1.0e-8)
254            self.assertAlmostEqual(pt.e_tot, -149.815015567956, delta=1.0e-8)
255
256    def test_scs_natorbs(self):
257        mol = self.mf.mol
258        with SCSUMP2(self.mf) as pt:
259            natocc, natorb = pt.make_natorbs()
260            # number of electrons conserved
261            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
262            # orbitals orthogonal
263            check_orth(self, mol, natorb)
264            # selected values
265            self.assertAlmostEqual(natocc[6], 1.9512132898, delta=1.0e-7)
266            self.assertAlmostEqual(natocc[7], 1.0092563934, delta=1.0e-7)
267            self.assertAlmostEqual(natocc[9], 0.0451631109, delta=1.0e-7)
268
269    def test_scs_natorbs_relaxed(self):
270        mol = self.mf.mol
271        with SCSUMP2(self.mf) as pt:
272            pt.cphf_tol = 1e-12
273            natocc, natorb = pt.make_natorbs(relaxed=True)
274            # number of electrons conserved
275            self.assertAlmostEqual(numpy.sum(natocc), mol.nelectron, delta=1.0e-10)
276            # orbitals orthogonal
277            check_orth(self, mol, natorb)
278            # selected values
279            self.assertAlmostEqual(natocc[6], 1.9516484920, delta=1.0e-7)
280            self.assertAlmostEqual(natocc[7], 1.0093068572, delta=1.0e-7)
281            self.assertAlmostEqual(natocc[9], 0.0434261850, delta=1.0e-7)
282
283
284if __name__ == "__main__":
285    print("Full Tests for native DF-UMP2")
286    unittest.main()
287