aotoso.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

//
// aotoso.cc --- more symmetry stuff
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Edward Seidl <[e-mail]>
// 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.
//

#include <float.h>

#include <util/misc/formio.h>

#include <chemistry/qc/basis/basis.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/basis/shellrot.h>
#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/basis/f77sym.h>

using namespace std;
using namespace sc;

extern "C" {
  void
  F77_DGESVD(const char * JOBU, const char *JOBVT,
             int *M, int *N, double *A, int *LDA,
             double *S, double *U, int *LDU, double *VT, int *LDVT,
             double *WORK, int *LWORK, int *INFO );
}

////////////////////////////////////////////////////////////////////////////

contribution::contribution()
{
}

contribution::~contribution()
{
}

contribution::contribution(int b, double c) : bfn(b), coef(c)
{
}

////////////////////////////////////////////////////////////////////////////

SO::SO() : len(0), length(0), cont(0)
{
}

SO::SO(int l) : len(0), length(0), cont(0)
{
  set_length(l);
}

SO::~SO()
{
  set_length(0);
}

SO&
SO::operator=(const SO& so)
{
  set_length(so.length);
  length = so.length;
  for (int i=0; i < length; i++)
    cont[i] = so.cont[i];
  return *this;
}

void
SO::set_length(int l)
{
  len=l;
  length=l;
  if (cont) {
    delete[] cont;
    cont=0;
  }

  if (l)
    cont = new contribution[l];
}

void
SO::reset_length(int l)
{
  length=l;

  if (l <= len)
    return;
  
  l=l+10;
  
  contribution *newcont = new contribution[l];
  
  if (cont) {
    for (int i=0; i < len; i++)
      newcont[i] = cont[i];
    
    delete[] cont;
  }

  cont=newcont;
  len=l;
}

int
SO::equiv(const SO& so)
{
  int i;
      
  if (so.length != length)
    return 0;

  double c=0;
  for (i=0; i < length; i++) {
    if (cont[i].bfn != so.cont[i].bfn)
      return 0;
    c += cont[i].coef*so.cont[i].coef;
  }
      
  // if the overlap == 1.0, they're equal (SO's should have been
  // normalized by now)
  if (fabs(fabs(c)-1.0) < 1.0e-3)
    return 1;

  return 0;
}

////////////////////////////////////////////////////////////////////////////

SO_block::SO_block() : len(0), so(0)
{
}

SO_block::SO_block(int l) : len(0), so(0)
{
  set_length(l);
}

SO_block::~SO_block()
{
  set_length(0);
}

void
SO_block::set_length(int l)
{
  len=l;
  if (so) {
    delete[] so;
    so=0;
  }

  if (l)
    so = new SO[l];
}

void
SO_block::reset_length(int l)
{
  if (len == l) return;

  SO *newso = new SO[l];
  
  if (so) {
    for (int i=0; i < len; i++)
      newso[i] = so[i];
    
    delete[] so;
  }

  so=newso;
  len=l;
}

int
SO_block::add(SO& s, int i)
{
  // first check to see if s is already here
  for (int j=0; j < ((i < len) ? i : len); j++)
    if (so[j].equiv(s))
      return 0;
      
  if (i >= len)
    reset_length(i+1);
  so[i] = s;

  return 1;
}

void
SO_block::print(const char *title)
{
  int i,j;
  ExEnv::out0() << indent << "SO block " << title << endl;
  for (i=0; i < len; i++) {
    ExEnv::out0() << indent << "SO " << i+1 << endl << indent;
    for (j=0; j < so[i].length; j++)
      ExEnv::out0() << scprintf(" %10d",so[i].cont[j].bfn);
    ExEnv::out0() << endl << indent;
    for (j=0; j < so[i].length; j++)
      ExEnv::out0() << scprintf(" %10.7f",so[i].cont[j].coef);
    ExEnv::out0() << endl;
  }
}

////////////////////////////////////////////////////////////////////////////

struct lin_comb {
    int ns;
    int f0;
    int mapf0;
    double **c;

    lin_comb(int ins, int if0, int imf0) : ns(ins), f0(if0), mapf0(imf0) {
      int i;
      
      c = new double*[ns];
      for (i=0; i < ns; i++) {
        c[i] = new double[ns];
        memset(c[i],0,sizeof(double)*ns);
      }
    }

    ~lin_comb() {
      if (c) {
        for (int i=0; i < ns; i++)
          if (c[i])
            delete[] c[i];
        delete[] c;
        c=0;
      }
    }

    void print() const {
      int i;
      ExEnv::out0() << indent;
      for (i=0; i < ns; i++)
        ExEnv::out0() << scprintf(" %10d",mapf0+i);
      ExEnv::out0() << endl;
      
      for (i=0; i < ns; i++) {
        ExEnv::out0() << indent << scprintf("%2d",f0+i);
        for (int j=0; j < ns; j++)
          ExEnv::out0() << scprintf(" %10.7f",c[i][j]);
        ExEnv::out0() << endl;
      }
    }
};

////////////////////////////////////////////////////////////////////////////

SO_block *
PetiteList::aotoso_info()
{
  int iuniq, i, j, d, ii, jj, g, s, c, ir;

  GaussianBasisSet& gbs = *gbs_.pointer();
  Molecule& mol = *gbs.molecule().pointer();

  // create the character table for the point group
  CharacterTable ct = mol.point_group()->char_table();
  SymmetryOperation so;

  if (c1_) {
    SO_block *SOs = new SO_block[1];
    SOs[0].set_length(gbs.nbasis());
    for (i=0; i < gbs.nbasis(); i++) {
      SOs[0].so[i].set_length(1);
      SOs[0].so[i].cont[0].bfn=i;
      SOs[0].so[i].cont[0].coef=1.0;
    }
    return SOs;
  }

  // ncomp is the number of symmetry blocks we have. for point groups with
  // complex E representations, this will be cut in two eventually
  int ncomp=0;
  for (i=0; i < nirrep_; i++)
    ncomp += ct.gamma(i).degeneracy();
  
  // saoelem is the current SO in a block
  int *saoelem = new int[ncomp];
  memset(saoelem,0,sizeof(int)*ncomp);

  int *whichir = new int[ncomp];
  int *whichcmp = new int[ncomp];
  for (i=ii=0; i < nirrep_; i++) {
    for (int j=0; j < ct.gamma(i).degeneracy(); j++,ii++) {
      whichir[ii] = i;
      whichcmp[ii] = j;
    }
  }
  
  // SOs is an array of SO_blocks which holds the redundant SO's
  SO_block *SOs = new SO_block[ncomp];

  for (i=0; i < ncomp; i++) {
    ir = whichir[i];
    int len = (ct.gamma(ir).complex()) ? nbf_in_ir_[ir]/2 : nbf_in_ir_[ir];
    SOs[i].set_length(len);
  }
  
  // loop over all unique shells
  for (iuniq=0; iuniq < gbs.molecule()->nunique(); iuniq++) {
    int nequiv = gbs.molecule()->nequivalent(iuniq);
    i = gbs.molecule()->unique(iuniq);
    for (s=0; s < gbs.nshell_on_center(i); s++) {
      int shell_i = gbs.shell_on_center(i,s);
      
      // test to see if there are any high am cartesian functions in this
      // shell.  for now don't allow symmetry with cartesian functions...I
      // just can't seem to get them working.
      for (c=0; c < gbs(i,s).ncontraction(); c++) {
        if (gbs(i,s).am(c) > 1 && gbs(i,s).is_cartesian(c)) {
          if (ng_ != nirrep_) {
            ExEnv::err0() << indent
                         << "PetiteList::aotoso: cannot yet handle"
                         << " symmetry for angular momentum >= 2\n";
            abort();
          }
        }
      }

      // the functions do not mix between contractions
      // so the contraction loop can be done outside the symmetry
      // operation loop
      int bfn_offset_in_shell = 0;
      for (c=0; c < gbs(i,s).ncontraction(); c++) {
        int nfuncuniq = gbs(i,s).nfunction(c);
        int nfuncall = nfuncuniq * nequiv;

        // allocate an array to store linear combinations of orbitals
        // this is destroyed by the SVD routine
        double **linorb = new double*[nfuncuniq];
        linorb[0] = new double[nfuncuniq*nfuncall];
        for (j=1; j<nfuncuniq; j++) {
          linorb[j] = &linorb[j-1][nfuncall];
        }

        // a copy of linorb
        double **linorbcop = new double*[nfuncuniq];
        linorbcop[0] = new double[nfuncuniq*nfuncall];
        for (j=1; j<nfuncuniq; j++) {
          linorbcop[j] = &linorbcop[j-1][nfuncall];
        }

        // allocate an array for the SVD routine
        double **u = new double*[nfuncuniq];
        u[0] = new double[nfuncuniq*nfuncuniq];
        for (j=1; j<nfuncuniq; j++) {
          u[j] = &u[j-1][nfuncuniq];
        }

        // loop through each irrep to form the linear combination
        // of orbitals of that symmetry
        int irnum = 0;
        for (ir=0; ir<ct.nirrep(); ir++) {
          int cmplx = (ct.complex() && ct.gamma(ir).complex());
          for (int comp=0; comp < ct.gamma(ir).degeneracy(); comp++) {
            memset(linorb[0], 0, nfuncuniq*nfuncall*sizeof(double));
            for (d=0; d < ct.gamma(ir).degeneracy(); d++) {
              // if this is a point group with a complex E representation,
              // then only do the first set of projections for E
              if (d && cmplx) break;

              // operate on each function in this contraction with each
              // symmetry operation to form symmetrized linear combinations
              // of orbitals

              for (g=0; g<ng_; g++) {
                // the character
                double ccdg = ct.gamma(ir).p(comp,d,g);

                so = ct.symm_operation(g);
                int equivatom = atom_map_[i][g];
                int atomoffset
                  = gbs.molecule()->atom_to_unique_offset(equivatom);
        
                ShellRotation rr
                  = ints_->shell_rotation(gbs(i,s).am(c),
                                          so,gbs(i,s).is_pure(c));

                for (ii=0; ii < rr.dim(); ii++) {
                  for (jj=0; jj < rr.dim(); jj++) {
                    linorb[ii][atomoffset*nfuncuniq+jj] += ccdg * rr(ii,jj);
                  }
                }
              }
            }