# petite.h

    00001 //
    00002 // petite.h --- definition of the PetiteList class
    00003 //
    00004 // Copyright (C) 1996 Limit Point Systems, Inc.
    00005 //
    00006 // Author: Edward Seidl <seidl@janed.com>
    00007 // Maintainer: LPS
    00008 //
    00009 // This file is part of the SC Toolkit.
    00010 //
    00011 // The SC Toolkit is free software; you can redistribute it and/or modify
    00012 // it under the terms of the GNU Library General Public License as published by
    00013 // the Free Software Foundation; either version 2, or (at your option)
    00014 // any later version.
    00015 //
    00016 // The SC Toolkit is distributed in the hope that it will be useful,
    00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
    00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    00019 // GNU Library General Public License for more details.
    00020 //
    00021 // You should have received a copy of the GNU Library General Public License
    00022 // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
    00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
    00024 //
    00025 // The U.S. Government is granted a limited license as per AL 91-7.
    00026 //
    00027
    00028 #ifndef _chemistry_qc_basis_petite_h
    00029 #define _chemistry_qc_basis_petite_h
    00030
    00031 #ifdef __GNUC__
    00032 #pragma interface
    00033 #endif
    00034
    00035 #include <scconfig.h>
    00036 #include <iostream>
    00037 #include <scconfig.h>
    00038
    00039 #include <util/misc/scint.h>
    00040 #include <util/ref/ref.h>
    00041 #include <math/scmat/blocked.h>
    00042 #include <math/scmat/offset.h>
    00043 #include <chemistry/molecule/molecule.h>
    00044 #include <chemistry/qc/basis/gaussbas.h>
    00045 #include <chemistry/qc/basis/integral.h>
    00046
    00047 // //////////////////////////////////////////////////////////////////////////
    00048
    00049 namespace sc {
    00050
    00051 inline sc_int_least64_t
    00052 ij_offset64(sc_int_least64_t i, sc_int_least64_t j)
    00053 {
    00054   return (i>j) ? (((i*(i+1)) >> 1) + j) : (((j*(j+1)) >> 1) + i);
    00055 }
    00056
    00057 inline sc_int_least64_t
    00058 i_offset64(sc_int_least64_t i)
    00059 {
    00060   return ((i*(i+1)) >> 1);
    00061 }
    00062
    00063 // //////////////////////////////////////////////////////////////////////////
    00064 // These are helper functions for PetiteList and GPetite4
    00065
    00066 int **compute_atom_map(const Ref<GaussianBasisSet> &);
    00067 void delete_atom_map(int **atom_map, const Ref<GaussianBasisSet> &);
    00068
    00069 int **compute_shell_map(int **atom_map, const Ref<GaussianBasisSet> &);
    00070 void delete_shell_map(int **shell_map, const Ref<GaussianBasisSet> &);
    00071
    00072 // //////////////////////////////////////////////////////////////////////////
    00073
    00074 struct contribution {
    00075     int bfn;
    00076     double coef;
    00077
    00078     contribution();
    00079     contribution(int b, double c);
    00080     ~contribution();
    00081 };
    00082
    00083 struct SO {
    00084     int len;
    00085     int length;
    00086     contribution *cont;
    00087
    00088     SO();
    00089     SO(int);
    00090     ~SO();
    00091
    00092     SO& operator=(const SO&);
    00093
    00094     void set_length(int);
    00095     void reset_length(int);
    00096
    00097     // is this equal to so to within a sign
    00098     int equiv(const SO& so);
    00099 };
    00100
    00101 struct SO_block {
    00102     int len;
    00103     SO *so;
    00104
    00105     SO_block();
    00106     SO_block(int);
    00107     ~SO_block();
    00108
    00109     void set_length(int);
    00110     void reset_length(int);
    00111
    00112     int add(SO& s, int i);
    00113     void print(const char *title);
    00114 };
    00115
    00116 // //////////////////////////////////////////////////////////////////////////
    00117 // this should only be used from within a SymmGaussianBasisSet
    00118
    00119 class PetiteList : public RefCount {
    00120   private:
    00121     int natom_;
    00122     int nshell_;
    00123     int ng_;
    00124     int nirrep_;
    00125     int nblocks_;
    00126     int c1_;
    00127
    00128     Ref<GaussianBasisSet> gbs_;
    00129     Ref<Integral> ints_;
    00130
    00131     char *p1_;        // p1[n] is 1 if shell n is in the group P1
    00132     int **atom_map_;  // atom_map[n][g] is the atom that symop g maps atom n
    00133                       // into
    00134     int **shell_map_; // shell_map[n][g] is the shell that symop g maps shell n
    00135                       // into
    00136     char *lamij_;     // see Dupuis & King, IJQC 11,613,(1977)
    00137
    00138     int *nbf_in_ir_;
    00139
    00140     void init();
    00141
    00142   public:
    00143     PetiteList(const Ref<GaussianBasisSet>&, const Ref<Integral>&);
    00144     ~PetiteList();
    00145
    00146     Ref<GaussianBasisSet> basis() { return gbs_; }
    00147     Ref<Integral> integral() { return ints_; }
    00148     Ref<PetiteList> clone() { return new PetiteList(gbs_, ints_); }
    00149
    00150     int nirrep() const { return nirrep_; }
    00151     int order() const { return ng_; }
    00152     int atom_map(int n, int g) const { return (c1_) ? n : atom_map_[n][g]; }
    00153     int shell_map(int n, int g) const { return (c1_) ? n : shell_map_[n][g]; }
    00154     int lambda(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
    00155     int lambda(int i, int j) const
    00156                           { return (c1_) ? 1 : (int) lamij_[ij_offset(i,j)]; }
    00157
    00158     int in_p1(int n) const { return (c1_) ? 1 : (int) p1_[n]; }
    00159     int in_p2(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
    00161     int in_p2(int i, int j) const { return (c1_) ? 1 : (int) lamij_[ij_offset(i,j)]; }
    00162     int in_p4(int ij, int kl, int i, int j, int k, int l) const;
    00164     int in_p4(int i, int j, int k, int l) const;
    00165
    00166     int nfunction(int i) const
    00167                             { return (c1_) ? gbs_->nbasis() : nbf_in_ir_[i]; }
    00168
    00169     int nblocks() const { return nblocks_; }
    00170
    00171     void print(std::ostream& =ExEnv::out0(), int verbose=1);
    00172
    00173     // these return blocked dimensions
    00174     RefSCDimension AO_basisdim();
    00175     RefSCDimension SO_basisdim();
    00176
    00177     // return the basis function rotation matrix R(g)
    00178     RefSCMatrix r(int g);
    00179
    00180     // return information about the transformation from AOs to SOs
    00181     SO_block * aotoso_info();
    00182     RefSCMatrix aotoso();
    00183     RefSCMatrix sotoao();
    00184
    00185     // given a skeleton matrix, form the symmetrized matrix in the SO basis
    00186     void symmetrize(const RefSymmSCMatrix& skel, const RefSymmSCMatrix& sym);
    00187
    00188     // transform a matrix from AO->SO or SO->AO.
    00189     // this can take either a blocked or non-blocked AO basis matrix.
    00190     RefSymmSCMatrix to_SO_basis(const RefSymmSCMatrix&);
    00191
    00192     // this returns a non-blocked AO basis matrix.
    00193     RefSymmSCMatrix to_AO_basis(const RefSymmSCMatrix&);
    00194
    00195     // these two are just for eigenvectors
    00196     // returns non-blocked AO basis eigenvectors
    00197     RefSCMatrix evecs_to_AO_basis(const RefSCMatrix&);
    00198     // returns blocked SO basis eigenvectors
    00199     RefSCMatrix evecs_to_SO_basis(const RefSCMatrix&);
    00200 };
    00201
    00202 inline int
    00203 PetiteList::in_p4(int ij, int kl, int i, int j, int k, int l) const
    00204 {
    00205   if (c1_)
    00206     return 1;
    00207
    00208   sc_int_least64_t ijkl = i_offset64(ij)+kl;
    00209   int nijkl=1;
    00210
    00211   for (int g=1; g < ng_; g++) {
    00212     int gij = ij_offset(shell_map_[i][g],shell_map_[j][g]);
    00213     int gkl = ij_offset(shell_map_[k][g],shell_map_[l][g]);
    00214     sc_int_least64_t gijkl = ij_offset64(gij,gkl);
    00215
    00216     if (gijkl > ijkl)
    00217       return 0;
    00218     else if (gijkl == ijkl)
    00219       nijkl++;
    00220   }
    00221
    00222   return ng_/nijkl;
    00223 }
    00224
    00225 inline int
    00226 PetiteList::in_p4(int i, int j, int k, int l) const
    00227 {
    00228   if (c1_)
    00229     return 1;
    00230
    00231   int ij = ij_offset(i,j);
    00232   int kl = ij_offset(k,l);
    00233   sc_int_least64_t ijkl = ij_offset64(ij,kl);
    00234   int nijkl=1;
    00235
    00236   for (int g=1; g < ng_; g++) {
    00237     int gij = ij_offset(shell_map_[i][g],shell_map_[j][g]);
    00238     int gkl = ij_offset(shell_map_[k][g],shell_map_[l][g]);
    00239     sc_int_least64_t gijkl = ij_offset64(gij,gkl);
    00240
    00241     if (gijkl > ijkl)
    00242       return 0;
    00243     else if (gijkl == ijkl)
    00244       nijkl++;
    00245   }
    00246
    00247   return ng_/nijkl;
    00248 }
    00249
    00250 }
    00251
    00252
    00253
    00254 #endif
    00255
    00256 // Local Variables:
    00257 // mode: c++
    00258 // c-file-style: "ETS"
    00259 // End: