1 //
2 // petite.h --- definition of the PetiteList class
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Edward Seidl <seidl@janed.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifndef _chemistry_qc_basis_petite_h
29 #define _chemistry_qc_basis_petite_h
30
31 #ifdef __GNUC__
32 #pragma interface
33 #endif
34
35 #include <iostream>
36
37 #include <util/misc/scint.h>
38 #include <util/ref/ref.h>
39 #include <util/container/array.h>
40 #include <math/scmat/blocked.h>
41 #include <math/scmat/offset.h>
42 #include <chemistry/molecule/molecule.h>
43 #include <chemistry/qc/basis/gaussbas.h>
44 #include <chemistry/qc/basis/integral.h>
45
46 // //////////////////////////////////////////////////////////////////////////
47
48 inline sc_int_least64_t
49 ij_offset64(sc_int_least64_t i, sc_int_least64_t j)
50 2739 0 1947 693 1344 327 10 14 0 0 {
51 return (i>j) ? (((i*(i+1)) >> 1) + j) : (((j*(j+1)) >> 1) + i);
52 }
53
54 inline sc_int_least64_t
55 i_offset64(sc_int_least64_t i)
56 592 24 379 125 237 56 6 5 0 0 {
57 return ((i*(i+1)) >> 1);
58 }
59
60 // //////////////////////////////////////////////////////////////////////////
61
62 struct contribution {
63 int bfn;
64 double coef;
65
66 contribution();
67 contribution(int b, double c);
68 ~contribution();
69 };
70
71 struct SO {
72 int len;
73 int length;
74 contribution *cont;
75
76 SO();
77 SO(int);
78 ~SO();
79
80 SO& operator=(const SO&);
81
82 void set_length(int);
83 void reset_length(int);
84
85 // is this equal to so to within a sign
86 int equiv(const SO& so);
87 };
88
89 struct SO_block {
90 int len;
91 SO *so;
92
93 SO_block();
94 SO_block(int);
95 ~SO_block();
96
97 void set_length(int);
98 void reset_length(int);
99
100 int add(SO& s, int i);
101 void print(const char *title);
102 };
103
104 // //////////////////////////////////////////////////////////////////////////
105 // this should only be used from within a SymmGaussianBasisSet
106
107 class PetiteList : public RefCount {
108 private:
109 int natom_;
110 int nshell_;
111 int ng_;
112 int nirrep_;
113 int nblocks_;
114 int c1_;
115
116 Ref<GaussianBasisSet> gbs_;
117 Ref<Integral> ints_;
118
119 char *p1_; // p1[n] is 1 if shell n is in the group P1
120 int **atom_map_; // atom_map[n][g] is the atom that symop g maps atom n
121 // into
122 int **shell_map_; // shell_map[n][g] is the shell that symop g maps shell n
123 // into
124 char *lamij_; // see Dupuis & King, IJQC 11,613,(1977)
125
126 int *nbf_in_ir_;
127
128 void init();
129
130 public:
131 PetiteList(const Ref<GaussianBasisSet>&, const Ref<Integral>&);
132 ~PetiteList();
133
134 int nirrep() const { return nirrep_; }
135 int order() const { return ng_; }
136 int atom_map(int n, int g) const { return (c1_) ? n : atom_map_[n][g]; }
137 int shell_map(int n, int g) const { return (c1_) ? n : shell_map_[n][g]; }
138 int lambda(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
139 int lambda(int i, int j) const
140 2 0 4 0 2 0 0 0 0 0 { return (c1_) ? 1 : (int) lamij_[ij_offset(i,j)]; }
141
142 4 0 1 0 2 0 0 0 0 0 int in_p1(int n) const { return (c1_) ? 1 : (int) p1_[n]; }
143 14 0 5 3 1 0 0 0 0 0 int in_p2(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
144 int in_p4(int ij, int kl, int i, int j, int k, int l) const;
145
146 int nfunction(int i) const
147 { return (c1_) ? gbs_->nbasis() : nbf_in_ir_[i]; }
148
149 int nblocks() const { return nblocks_; }
150
151 void print(std::ostream& =ExEnv::out(), int verbose=1);
152
153 // these return blocked dimensions
154 RefSCDimension AO_basisdim();
155 RefSCDimension SO_basisdim();
156
157 // return the basis function rotation matrix R(g)
158 RefSCMatrix r(int g);
159
160 // return information about the transformation from AOs to SOs
161 SO_block * aotoso_info();
162 RefSCMatrix aotoso();
163 RefSCMatrix sotoao();
164
165 // given a skeleton matrix, form the symmetrized matrix in the SO basis
166 void symmetrize(const RefSymmSCMatrix& skel, const RefSymmSCMatrix& sym);
167
168 // transform a matrix from AO->SO or SO->AO.
169 // this can take either a blocked or non-blocked AO basis matrix.
170 RefSymmSCMatrix to_SO_basis(const RefSymmSCMatrix&);
171
172 // this returns a non-blocked AO basis matrix.
173 RefSymmSCMatrix to_AO_basis(const RefSymmSCMatrix&);
174
175 // these two are just for eigenvectors
176 // returns non-blocked AO basis eigenvectors
177 RefSCMatrix evecs_to_AO_basis(const RefSCMatrix&);
178 // returns blocked SO basis eigenvectors
179 RefSCMatrix evecs_to_SO_basis(const RefSCMatrix&);
180 };
181
182 inline int
183 PetiteList::in_p4(int ij, int kl, int i, int j, int k, int l) const
184 {
185 298 10 110 34 66 17 1 3 0 0 if (c1_)
186 return 1;
187
188 sc_int_least64_t ijkl = i_offset64(ij)+kl;
189 147 0 70 16 39 13 3 0 0 0 int nijkl=0;
190
191 672 2 277 55 184 13 0 5 0 0 for (int g=0; g < ng_; g++) {
192 1330 0 766 357 457 76 0 12 0 0 int gij = ij_offset(shell_map_[i][g],shell_map_[j][g]);
193 1088 0 752 288 534 102 3 20 0 0 int gkl = ij_offset(shell_map_[k][g],shell_map_[l][g]);
194 sc_int_least64_t gijkl = ij_offset64(gij,gkl);
195
196 281 0 247 80 186 46 0 1 0 0 if (gijkl > ijkl)
197 10 0 11 4 4 2 0 0 0 0 return 0;
198 211 0 143 36 78 21 0 0 0 0 else if (gijkl == ijkl)
199 172 0 119 42 68 17 0 1 0 0 nijkl++;
200 915 0 831 226 510 88 2 1 0 0 }
201
202 return ng_/nijkl;
203 }
204
205
206
207 #endif
208
209 // Local Variables:
210 // mode: c++
211 // c-file-style: "ETS"
212 // End:
213