petite.h
//
// petite.h --- definition of the PetiteList class
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Edward Seidl <seidl@janed.com>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB. If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//
#ifndef _chemistry_qc_basis_petite_h
#define _chemistry_qc_basis_petite_h
#ifdef __GNUC__
#pragma interface
#endif
#include <scconfig.h>
#include <iostream>
#include <scconfig.h>
#include <util/misc/scint.h>
#include <util/ref/ref.h>
#include <math/scmat/blocked.h>
#include <math/scmat/offset.h>
#include <chemistry/molecule/molecule.h>
#include <chemistry/qc/basis/gaussbas.h>
#include <chemistry/qc/basis/integral.h>
// //////////////////////////////////////////////////////////////////////////
namespace sc {
inline sc_int_least64_t
ij_offset64(sc_int_least64_t i, sc_int_least64_t j)
{
return (i>j) ? (((i*(i+1)) >> 1) + j) : (((j*(j+1)) >> 1) + i);
}
inline sc_int_least64_t
i_offset64(sc_int_least64_t i)
{
return ((i*(i+1)) >> 1);
}
// //////////////////////////////////////////////////////////////////////////
// These are helper functions for PetiteList and GPetite4
int **compute_atom_map(const Ref<GaussianBasisSet> &);
void delete_atom_map(int **atom_map, const Ref<GaussianBasisSet> &);
int **compute_shell_map(int **atom_map, const Ref<GaussianBasisSet> &);
void delete_shell_map(int **shell_map, const Ref<GaussianBasisSet> &);
// //////////////////////////////////////////////////////////////////////////
struct contribution {
int bfn;
double coef;
contribution();
contribution(int b, double c);
~contribution();
};
struct SO {
int len;
int length;
contribution *cont;
SO();
SO(int);
~SO();
SO& operator=(const SO&);
void set_length(int);
void reset_length(int);
// is this equal to so to within a sign
int equiv(const SO& so);
};
struct SO_block {
int len;
SO *so;
SO_block();
SO_block(int);
~SO_block();
void set_length(int);
void reset_length(int);
int add(SO& s, int i);
void print(const char *title);
};
// //////////////////////////////////////////////////////////////////////////
// this should only be used from within a SymmGaussianBasisSet
class PetiteList : public RefCount {
private:
int natom_;
int nshell_;
int ng_;
int nirrep_;
int nblocks_;
int c1_;
Ref<GaussianBasisSet> gbs_;
Ref<Integral> ints_;
char *p1_; // p1[n] is 1 if shell n is in the group P1
int **atom_map_; // atom_map[n][g] is the atom that symop g maps atom n
// into
int **shell_map_; // shell_map[n][g] is the shell that symop g maps shell n
// into
char *lamij_; // see Dupuis & King, IJQC 11,613,(1977)
int *nbf_in_ir_;
void init();
public:
PetiteList(const Ref<GaussianBasisSet>&, const Ref<Integral>&);
~PetiteList();
Ref<GaussianBasisSet> basis() { return gbs_; }
Ref<Integral> integral() { return ints_; }
Ref<PetiteList> clone() { return new PetiteList(gbs_, ints_); }
int nirrep() const { return nirrep_; }
int order() const { return ng_; }
int atom_map(int n, int g) const { return (c1_) ? n : atom_map_[n][g]; }
int shell_map(int n, int g) const { return (c1_) ? n : shell_map_[n][g]; }
int lambda(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
int lambda(int i, int j) const
{ return (c1_) ? 1 : (int) lamij_[ij_offset(i,j)]; }
int in_p1(int n) const { return (c1_) ? 1 : (int) p1_[n]; }
int in_p2(int ij) const { return (c1_) ? 1 : (int) lamij_[ij]; }
int in_p2(int i, int j) const { return (c1_) ? 1 : (int) lamij_[ij_offset(i,j)]; }
int in_p4(int ij, int kl, int i, int j, int k, int l) const;
int in_p4(int i, int j, int k, int l) const;
int nfunction(int i) const
{ return (c1_) ? gbs_->nbasis() : nbf_in_ir_[i]; }
int nblocks() const { return nblocks_; }
void print(std::ostream& =ExEnv::out0(), int verbose=1);
// these return blocked dimensions
RefSCDimension AO_basisdim();
RefSCDimension SO_basisdim();
// return the basis function rotation matrix R(g)
RefSCMatrix r(int g);
// return information about the transformation from AOs to SOs
SO_block * aotoso_info();
RefSCMatrix aotoso();
RefSCMatrix sotoao();
// given a skeleton matrix, form the symmetrized matrix in the SO basis
void symmetrize(const RefSymmSCMatrix& skel, const RefSymmSCMatrix& sym);
// transform a matrix from AO->SO or SO->AO.
// this can take either a blocked or non-blocked AO basis matrix.
RefSymmSCMatrix to_SO_basis(const RefSymmSCMatrix&);
// this returns a non-blocked AO basis matrix.
RefSymmSCMatrix to_AO_basis(const RefSymmSCMatrix&);
// these two are just for eigenvectors
// returns non-blocked AO basis eigenvectors
RefSCMatrix evecs_to_AO_basis(const RefSCMatrix&);
// returns blocked SO basis eigenvectors
RefSCMatrix evecs_to_SO_basis(const RefSCMatrix&);
};
inline int
PetiteList::in_p4(int ij, int kl, int i, int j, int k, int l) const
{
if (c1_)
return 1;
sc_int_least64_t ijkl = i_offset64(ij)+kl;
int nijkl=1;
for (int g=1; g < ng_; g++) {
int gij = ij_offset(shell_map_[i][g],shell_map_[j][g]);
int gkl = ij_offset(shell_map_[k][g],shell_map_[l][g]);
sc_int_least64_t gijkl = ij_offset64(gij,gkl);
if (gijkl > ijkl)
return 0;
else if (gijkl == ijkl)
nijkl++;
}
return ng_/nijkl;
}
inline int
PetiteList::in_p4(int i, int j, int k, int l) const
{
if (c1_)
return 1;
int ij = ij_offset(i,j);
int kl = ij_offset(k,l);
sc_int_least64_t ijkl = ij_offset64(ij,kl);
int nijkl=1;
for (int g=1; g < ng_; g++) {
int gij = ij_offset(shell_map_[i][g],shell_map_[j][g]);
int gkl = ij_offset(shell_map_[k][g],shell_map_[l][g]);
sc_int_least64_t gijkl = ij_offset64(gij,gkl);
if (gijkl > ijkl)
return 0;
else if (gijkl == ijkl)
nijkl++;
}
return ng_/nijkl;
}
}
#endif
// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: