1 //
2 // petite.cc --- implementation 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 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <stdio.h>
33
34 #include <util/misc/formio.h>
35 #include <chemistry/molecule/localdef.h>
36
37 #include <chemistry/qc/basis/basis.h>
38 #include <chemistry/qc/basis/integral.h>
39 #include <chemistry/qc/basis/petite.h>
40 #include <chemistry/qc/basis/shellrot.h>
41
42 using namespace std;
43
44 ////////////////////////////////////////////////////////////////////////////
45
46 PetiteList::PetiteList(const Ref<GaussianBasisSet> &gbs,
47 const Ref<Integral>& ints) :
48 gbs_(gbs),
49 ints_(ints)
50 {
51 init();
52 }
53
54 PetiteList::~PetiteList()
55 {
56 if (p1_)
57 delete[] p1_;
58
59 if (lamij_)
60 delete[] lamij_;
61
62 if (nbf_in_ir_)
63 delete[] nbf_in_ir_;
64
65 if (atom_map_) {
66 for (int i=0; i < natom_; i++)
67 delete[] atom_map_[i];
68 delete[] atom_map_;
69 }
70
71 if (shell_map_) {
72 1 0 0 0 0 0 0 0 0 0 for (int i=0; i < nshell_; i++)
73 1 0 0 0 0 0 0 0 0 0 delete[] shell_map_[i];
74 delete[] shell_map_;
75 }
76
77 natom_=0;
78 nshell_=0;
79 ng_=0;
80 nblocks_=0;
81 nirrep_=0;
82 p1_=0;
83 atom_map_=0;
84 shell_map_=0;
85 lamij_=0;
86 nbf_in_ir_=0;
87 }
88
89 void
90 PetiteList::init()
91 {
92 int i;
93
94 // grab references to the Molecule and BasisSet for convenience
95 GaussianBasisSet& gbs = *gbs_.pointer();
96 Molecule& mol = *gbs.molecule().pointer();
97
98 // create the character table for the point group
99 CharacterTable ct = mol.point_group()->char_table();
100
101 // initialize private members
102 c1_=0;
103 ng_ = ct.order();
104 natom_ = mol.natom();
105 nshell_ = gbs.nshell();
106 nirrep_ = ct.nirrep();
107
108 // if point group is C1, then zero everything
109 if (ng_==1) {
110 c1_=1;
111 nblocks_=1;
112
113 p1_=0;
114 atom_map_=0;
115 shell_map_=0;
116 lamij_=0;
117 nbf_in_ir_=0;
118 return;
119 }
120
121 // allocate storage for arrays
122 p1_ = new char[nshell_];
123 lamij_ = new char[i_offset(nshell_)];
124
125 atom_map_ = new int*[natom_];
126 for (i=0; i < natom_; i++)
127 0 0 0 1 0 0 0 0 0 0 atom_map_[i] = new int[ng_];
128
129 shell_map_ = new int*[nshell_];
130 for (i=0; i < nshell_; i++)
131 1 0 0 0 0 0 0 0 0 0 shell_map_[i] = new int[ng_];
132
133 // set up atom and shell mappings
134 double np[3];
135 0 0 0 1 0 0 0 0 0 0 SymmetryOperation so;
136
137 // loop over all centers
138 for (i=0; i < natom_; i++) {
139 SCVector3 ac(mol.r(i));
140
141 // then for each symop in the pointgroup, transform the coordinates of
142 // center "i" and see which atom it maps into
143 for (int g=0; g < ng_; g++) {
144 so = ct.symm_operation(g);
145
146 for (int ii=0; ii < 3; ii++) {
147 np[ii] = 0;
148 for (int jj=0; jj < 3; jj++)
149 1 0 0 0 0 0 0 0 0 0 np[ii] += so(ii,jj) * ac[jj];
150 }
151
152 0 0 0 0 1 0 0 0 0 0 atom_map_[i][g] = mol.atom_at_position(np, 0.05);
153 if (atom_map_[i][g] < 0) {
154 ExEnv::out() << "ERROR: Symmetry operation " << g << " did not map atom "
155 << i+1 << " to another atom:" << endl;
156 ExEnv::out() << indent << "Molecule:" << endl;
157 ExEnv::out() << incindent;
158 mol.print();
159 ExEnv::out() << decindent;
160 ExEnv::out() << indent << "attempted to find atom at" << endl;
161 ExEnv::out() << incindent;
162 ExEnv::out() << indent << np[0] << " " << np[1] << " " << np[2] << endl;
163 abort();
164 }
165 }
166
167 // hopefully, shells on equivalent centers will be numbered in the same
168 // order
169 for (int s=0; s < gbs.nshell_on_center(i); s++) {
170 int shellnum = gbs.shell_on_center(i,s);
171 1 0 0 0 0 0 0 0 0 0 for (int g=0; g < ng_; g++) {
172 1 0 3 0 5 0 0 0 0 0 shell_map_[shellnum][g] = gbs.shell_on_center(atom_map_[i][g],s);
173 }
174 1 0 1 0 0 0 0 0 0 0 }
175 }
176
177 memset(p1_,0,nshell_);
178 memset(lamij_,0,i_offset(nshell_));
179
180 // now we do p1_ and lamij_
181 1 0 0 0 0 0 0 0 0 0 for (i=0; i < nshell_; i++) {
182 int g;
183
184 // we want the highest numbered shell in a group of equivalent shells
185 for (g=0; g < ng_; g++)
186 if (shell_map_[i][g] > i)
187 break;
188
189 if (g < ng_)
190 continue;
191
192 // i is in the group P1
193 0 0 1 0 0 0 0 0 0 0 p1_[i] = 1;
194
195 1 0 0 0 0 0 0 0 0 0 for (int j=0; j <= i; j++) {
196 1 0 0 0 0 0 0 0 0 0 int ij = i_offset(i)+j;
197 int nij = 0;
198
199 // test to see if IJ is in the group P2, if it is, then set lambda(ij)
200 // equal to the number of equivalent shell pairs. This number is
201 // just the order of the group divided by the number of times ij is
202 // mapped into itself
203 int gg;
204 0 0 2 0 0 0 0 0 0 0 for (gg=0; gg < ng_; gg++) {
205 1 0 3 2 2 1 0 0 0 0 int gi = shell_map_[i][gg];
206 1 0 0 1 0 1 0 0 0 0 int gj = shell_map_[j][gg];
207 7 0 3 2 2 0 0 0 0 0 int gij = ij_offset(gi,gj);
208 0 0 1 0 1 0 0 0 0 0 if (gij > ij)
209 break;
210 1 0 1 0 1 0 0 0 0 0 else if (gij == ij)
211 2 0 3 2 0 1 0 0 0 0 nij++;
212 1 0 0 0 0 0 0 0 0 0 }
213
214 if (gg < ng_)
215 continue;
216
217 7 0 5 2 5 2 0 0 0 0 lamij_[ij] = (char) (ng_/nij);
218 0 0 1 0 0 0 0 0 0 0 }
219 1 0 0 0 1 0 0 0 0 0 }
220
221 // form reducible representation of the basis functions
222 double *red_rep = new double[ng_];
223 memset(red_rep,0,sizeof(double)*ng_);
224
225 for (i=0; i < natom_; i++) {
226 for (int g=0; g < ng_; g++) {
227 so = ct.symm_operation(g);
228 1 0 0 0 0 0 0 0 0 0 int j= atom_map_[i][g];
229
230 if (i!=j)
231 continue;
232
233 0 0 0 0 1 0 0 0 0 0 for (int s=0; s < gbs.nshell_on_center(i); s++) {
234 for (int c=0; c < gbs(i,s).ncontraction(); c++) {
235 int am=gbs(i,s).am(c);
236
237 if (am==0)
238 red_rep[g] += 1.0;
239 else {
240 ShellRotation r(am,so,ints_,gbs(i,s).is_pure(c));
241 red_rep[g] += r.trace();
242 }
243 }
244 2 0 0 1 1 1 0 0 0 0 }
245 }
246 }
247
248 // and then use projection operators to figure out how many SO's of each
249 // symmetry type there will be
250 nblocks_ = 0;
251 nbf_in_ir_ = new int[nirrep_];
252 for (i=0; i < nirrep_; i++) {
253 double t=0;
254 for (int g=0; g < ng_; g++)
255 t += ct.gamma(i).character(g)*red_rep[g];
256
257 0 0 2 0 0 0 0 0 0 0 nbf_in_ir_[i] = ((int) (t+0.5))/ng_;
258 if (ct.gamma(i).complex()) {
259 nblocks_++;
260 nbf_in_ir_[i] *= 2;
261 } else {
262 nblocks_ += ct.gamma(i).degeneracy();
263 }
264 }
265
266 delete[] red_rep;
267 }
268
269 RefSCDimension
270 PetiteList::AO_basisdim()
271 {
272 if (c1_)
273 return SO_basisdim();
274
275 RefSCDimension dim = new SCDimension(gbs_->nbasis(),1);
276 dim->blocks()->set_subdim(0, gbs_->basisdim());
277 return dim;
278 }
279
280 RefSCDimension
281 PetiteList::SO_basisdim()
282 {
283 int i, j, ii;
284
285 // grab a reference to the basis set
286 GaussianBasisSet& gbs = *gbs_.pointer();
287
288 // create the character table for the point group
289 CharacterTable ct = gbs.molecule()->point_group()->char_table();
290
291 // ncomp is the number of symmetry blocks we have
292 int ncomp=nblocks();
293
294 // saoelem is the current SO in a block
295 int *nao = new int [ncomp];
296 memset(nao,0,sizeof(int)*ncomp);
297
298 if (c1_)
299 nao[0] = gbs.nbasis();
300 else {
301 for (i=ii=0; i < nirrep_; i++) {
302 int je = ct.gamma(i).complex() ? 1 : ct.gamma(i).degeneracy();
303 for (j=0; j < je; j++,ii++)
304 nao[ii] = nbf_in_ir_[i];
305 0 0 0 0 1 0 0 0 0 0 }
306 }
307
308 RefSCDimension ret = new SCDimension(gbs.nbasis(),ncomp,nao);
309 delete[] nao;
310
311 for (i=ii=0; i < nirrep_; i++) {
312 Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
313 int np=grp->n();
314 2 0 0 1 0 0 0 0 0 0 int *subblksize = new int[np];
315 0 0 1 0 0 0 0 0 0 0 int nbas=(c1_) ? gbs.nbasis() : nbf_in_ir_[i];
316 for (j=0; j < np; j++) {
317 1 0 0 0 0 0 0 0 0 0 if (j < nbas%np)
318 subblksize[j] = nbas/np + 1;
319 else
320 subblksize[j] = nbas/np;
321 }
322
323 int je = ct.gamma(i).complex() ? 1 : ct.gamma(i).degeneracy();
324 for (j=0; j < je; j++,ii++) {
325 char lab[24];
326 sprintf(lab,"irrep %s comp %d", ct.gamma(i).symbol(), j);
327 ret->blocks()->set_subdim(ii, new SCDimension(nbas, np, subblksize));
328 }
329
330 delete[] subblksize;
331 }
332
333 return ret;
334 }
335
336 void
337 PetiteList::print(ostream& o, int verbose)
338 {
339 int i;
340
341 o << node0 << indent << "PetiteList:" << endl << incindent;
342
343 if (c1_) {
344 o << node0 << indent << "is c1\n" << decindent;
345 return;
346 }
347
348 if (verbose) {
349 o << node0
350 << indent << "natom_ = " << natom_ << endl
351 << indent << "nshell_ = " << nshell_ << endl
352 << indent << "ng_ = " << ng_ << endl
353 << indent << "nirrep_ = " << nirrep_ << endl << endl
354 << indent << "atom_map_ =" << endl << incindent;
355
356 for (i=0; i < natom_; i++) {
357 o << node0 << indent;
358 for (int g=0; g < ng_; g++)
359 o << node0 << scprintf("%5d ",atom_map_[i][g]);
360 o << node0 << endl;
361 }
362
363 o << node0 << endl << decindent
364 << indent << "shell_map_ =" << endl << incindent;
365 for (i=0; i < nshell_; i++) {
366 o << node0 << indent;
367 for (int g=0; g < ng_; g++)
368 o << node0 << scprintf("%5d ",shell_map_[i][g]);
369 o << node0 << endl;
370 }
371
372 o << node0 << endl << decindent
373 << indent << "p1_ =" << endl << incindent;
374 for (i=0; i < nshell_; i++)
375 o << node0 << indent << scprintf("%5d\n",p1_[i]);
376
377 o << node0 << decindent
378 << indent << "lamij_ =" << endl << incindent;
379 for (i=0; i < nshell_; i++) {
380 o << node0 << indent;
381 for (int j=0; j <= i; j++)
382 o << node0 << scprintf("%5d ",lamij_[i_offset(i)+j]);
383 o << node0 << endl;
384 }
385 o << node0 << endl << decindent;
386 }
387
388 CharacterTable ct = gbs_->molecule()->point_group()->char_table();
389 for (i=0; i < nirrep_; i++)
390 o << node0 << indent
391 << scprintf("%5d functions of %s symmetry\n",
392 nbf_in_ir_[i], ct.gamma(i).symbol());
393 }
394
395 // forms the basis function rotation matrix for the g'th symmetry operation
396 // in the point group
397 RefSCMatrix
398 PetiteList::r(int g)
399 {
400 // grab a reference to the basis set
401 GaussianBasisSet& gbs = *gbs_.pointer();
402
403 SymmetryOperation so =
404 gbs.molecule()->point_group()->char_table().symm_operation(g);
405
406 RefSCMatrix ret = gbs.matrixkit()->matrix(gbs.basisdim(), gbs.basisdim());
407 ret.assign(0.0);
408
409 // this should be replaced with an element op at some point
410 if (c1_) {
411 for (int i=0; i < gbs.nbasis(); i++)
412 ret.set_element(i,i,1.0);
413 return ret;
414
415 } else {
416 for (int i=0; i < natom_; i++) {
417 int j = atom_map_[i][g];
418
419 for (int s=0; s < gbs.nshell_on_center(i); s++) {
420 int func_i = gbs.shell_to_function(gbs.shell_on_center(i,s));
421 int func_j = gbs.shell_to_function(gbs.shell_on_center(j,s));
422
423 for (int c=0; c < gbs(i,s).ncontraction(); c++) {
424 int am=gbs(i,s).am(c);
425
426 if (am==0) {
427 ret.set_element(func_j,func_i,1.0);
428 } else {
429 ShellRotation rr(am,so,ints_,gbs(i,s).is_pure(c));
430 for (int ii=0; ii < rr.dim(); ii++)
431 for (int jj=0; jj < rr.dim(); jj++)
432 ret.set_element(func_j+jj,func_i+ii,rr(ii,jj));
433 }
434
435 func_i += gbs(i,s).nfunction(c);
436 func_j += gbs(i,s).nfunction(c);
437 }
438 }
439 }
440 }
441 return ret;
442 }
443
444 /////////////////////////////////////////////////////////////////////////////
445
446 // Local Variables:
447 // mode: c++
448 // c-file-style: "ETS"
449 // End:
450