1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: petite.cc
4 // Copyright (C) 2009 Toru Shiozaki
5 //
6 // Author: Toru Shiozaki <shiozaki@northwestern.edu>
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // This program is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program 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 General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program. If not, see <http://www.gnu.org/licenses/>.
23 //
24
25
26 #include <src/molecule/petite.h>
27 #include <array>
28 #include <memory>
29 #include <iostream>
30 #include <cassert>
31 #include <cmath>
32 #include <algorithm>
33
34 using namespace std;
35 using namespace bagel;
36
matmul33(const vector<double> & a,const array<double,3> & b)37 static inline array<double,3> matmul33(const vector<double>& a, const array<double,3>& b) {
38 assert(a.size() == 9 && b.size() == 3);
39 // note that the array is writen in C style (not in fortran style)!!!
40 array<double,3> out;
41 out[0] = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
42 out[1] = a[3] * b[0] + a[4] * b[1] + a[5] * b[2];
43 out[2] = a[6] * b[0] + a[7] * b[1] + a[8] * b[2];
44 return out;
45 };
46
Petite(const vector<shared_ptr<const Atom>> & atoms,const string sym)47 Petite::Petite(const vector<shared_ptr<const Atom>>& atoms, const string sym) : sym_(sym) {
48 const string c1("c1");
49 const string cs("cs");
50 const string ci("ci");
51 const string c2("c2");
52 const string d2("d2");
53 const string c2h("c2h");
54 const string c2v("c2v");
55 const string d2h("d2h");
56
57 natom_ = atoms.size();
58
59 vector<shared_ptr<const Shell>> vbb;
60 vector<int> offset;
61 int cnt = 0;
62 for (auto& a : atoms) {
63 vector<shared_ptr<const Shell>> tmp = a->shells();
64 vbb.insert(vbb.end(), tmp.begin(), tmp.end());
65 offset.push_back(cnt);
66 cnt += tmp.size();
67 }
68 nshell_ = vbb.size();
69
70 if (sym == c1) {
71 nirrep_ = 1;
72 } else {
73 if (sym == c2v){
74 SymC2v datc2v;
75 symop_ = datc2v.symop();
76 nirrep_ = datc2v.nirrep();
77 } else if (sym == d2h) {
78 SymD2h datd2h;
79 symop_ = datd2h.symop();
80 nirrep_ = datd2h.nirrep();
81 } else if (sym == cs) {
82 SymCs datcs;
83 symop_ = datcs.symop();
84 nirrep_ = datcs.nirrep();
85 } else if (sym == ci) {
86 SymCi datci;
87 symop_ = datci.symop();
88 nirrep_ = datci.nirrep();
89 } else if (sym == c2) {
90 SymC2 datc2;
91 symop_ = datc2.symop();
92 nirrep_ = datc2.nirrep();
93 } else if (sym == d2) {
94 SymD2 datd2;
95 symop_ = datd2.symop();
96 nirrep_ = datd2.nirrep();
97 } else if (sym == c2h) {
98 SymC2h datc2h;
99 symop_ = datc2h.symop();
100 nirrep_ = datc2h.nirrep();
101 } else {
102 assert(false);
103 }
104 nsymop_ = symop_.size();
105
106 // making map for atoms
107 for (int iatom = 0; iatom != natom_; ++iatom) {
108 const array<double,3> position = atoms[iatom]->position();
109 vector<int> tmp(nsymop_);
110
111 for (int iop = 0; iop != nsymop_; ++iop) {
112 const array<double,3> target = matmul33(symop_[iop], position);
113 bool found = false;
114 for (int jatom = 0; jatom != natom_; ++jatom) {
115 const array<double,3> current = atoms[jatom]->position();
116 if (current == target) {
117 found = true;
118 tmp[iop] = jatom;
119 break;
120 }
121 }
122 if (!found)
123 throw logic_error("Petite constructor error");
124 }
125 sym_atommap_.push_back(tmp);
126
127 // making map for shells
128 for (int i = 0; i != atoms[iatom]->nshell(); ++i) {
129 vector<int> stmp(nsymop_);
130 for (int iop = 0; iop != nsymop_; ++iop) {
131 stmp[iop] = offset[sym_atommap_[iatom][iop]] + i;
132 }
133 sym_shellmap_.push_back(stmp);
134 }
135 } // end of atom loop
136
137 // now we determine p1 and p2
138 p1_.resize(nshell_);
139 lambda_.resize(nshell_ * nshell_);
140
141 fill(p1_.begin(), p1_.end(), 0);
142
143 for (int i = 0; i != nshell_; ++i) {
144 bool skipp1 = false;
145 for (int iop = 0; iop != nsymop_; ++iop)
146 if (sym_shellmap_[i][iop] < i) skipp1 = true;
147 if (skipp1) continue;
148
149 p1_[i] = 1;
150
151 for (int j = i; j != nshell_; ++j) {
152 const int ij = i * nshell_ + j;
153 int nij = 0;
154
155 bool skipp2 = false;
156 for (int iop = 0; iop != nsymop_; ++iop) {
157 const int ci = sym_shellmap_[i][iop];
158 const int cj = sym_shellmap_[j][iop];
159 const int cij = ci * nshell_ + cj;
160 if (ij == cij) ++nij;
161 else if (cij < ij) skipp2 = true;
162 }
163 if (skipp2) continue;
164 lambda_[ij] = (nsymop_ / nij > 0) ? true : false;
165 }
166 } // end of shell loop
167
168 } // end of if_c1
169
170 }
171
172
~Petite()173 Petite::~Petite() {
174
175 }
176
177