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
16import numpy
17from pyscf import gto
18
19mol = gto.Mole()
20mol.atom = [
21    ['N', (0.,0.,0.)],
22    ['H', (0.,1.,1.)],
23    ['H', (1.,0.,1.)],
24    ['H', (1.,1.,0.)], ]
25mol.basis = {
26    "N": [(0, 0, (15, 1)), ],
27    "H": [(0, 0, (1, 1, 0), (3, 3, 1), (5, 1, 0)),
28          (1, 0, (1, 1)), (2, 0, (.8, 1)), ]}
29mol.basis['N'].extend(gto.mole.expand_etbs(((0, 4, 1, 1.8),
30                                            (1, 3, 2, 1.8),
31                                            (2, 2, 1, 1.8),)))
32
33mol.verbose = 0
34mol.output = None
35mol.build()
36
37tao = mol.time_reversal_map()
38
39# . . . /
40# . . ./
41#     /. .
42#    / . .
43#   /  . .
44#
45def rotatesub(smat, i0, j0, tao):
46    di = abs(tao[i0]) - i0
47    dj = abs(tao[j0]) - j0
48    rmat = numpy.empty((dj,di),dtype=smat.dtype)
49    for i in range(di):
50        for j in range(dj):
51            rmat[dj-1-j,di-1-i] = smat[i,j]
52    if tao[j0] < 0: # |j,-j> => -|j,j>
53        rmat[1::2] = -rmat[1::2]
54    else:
55        rmat[::2] = -rmat[::2]
56    if tao[i0] < 0:
57        rmat[:,1::2] = -rmat[:,1::2]
58    else:
59        rmat[:,::2] = -rmat[:,::2]
60    return rmat
61
62def rotatesub1(smat, i0, j0, tao):
63    di = abs(tao[i0]) - i0
64    dj = abs(tao[j0]) - j0
65    rmat = numpy.empty((dj,di),dtype=smat.dtype)
66    for j in range(dj):
67        for i in range(di):
68            rmat[j,i] = smat[di-1-i,dj-1-j]
69    if tao[j0] < 0: # |j,-j> => -|j,j>
70        rmat[1::2] = -rmat[1::2]
71    else:
72        rmat[::2] = -rmat[::2]
73    if tao[i0] < 0:
74        rmat[:,1::2] = -rmat[:,1::2]
75    else:
76        rmat[:,::2] = -rmat[:,::2]
77    return rmat
78
79if __name__ == '__main__':
80    ao_loc = mol.ao_loc_2c()
81    s = mol.intor('cint1e_spnucsp')
82    for ish in range(mol.nbas):
83        for jsh in range(mol.nbas):
84            i0 = ao_loc[ish]
85            i1 = abs(tao[i0])
86            j0 = ao_loc[jsh]
87            j1 = abs(tao[j0])
88            rmat = rotatesub(s[i0:i1,j0:j1], i0, j0, tao)
89            dd = abs(s[j0:j1,i0:i1]-rmat).sum()
90            assert(numpy.allclose(rmat,
91                                  rotatesub1(s[i0:i1,j0:j1], i0, j0, tao)))
92            if dd > 1e-12:
93                print(ish, jsh, dd)
94                assert(0)
95    print('pass')
96
97    nao = mol.nao_2c()
98    j = 0
99    idx0 = []
100    for i in range(nao):
101        if abs(tao[i]) > j:
102            idx0.append(i)
103            j = abs(tao[i])
104    j = 0
105    idx1 = [0]
106    while idx1[-1] < nao:
107        j += 1
108        idx1.append(abs(tao[idx1[-1]]))
109
110    print(all(numpy.array(idx0[:j-1]) == numpy.array(idx1[:j-1])))
111