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:

Generated at Mon Dec 3 23:23:40 2007 for MPQC 2.3.1 using the documentation package Doxygen 1.5.2.