nao.cc

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

//
// nao.cc
//
// Copyright (C) 1997 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[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 <util/misc/formio.h>
#include <chemistry/qc/wfn/wfn.h>
#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/basis/transform.h>

using namespace std;
using namespace sc;

#undef DEBUG

namespace sc {
static RefSCMatrix
operator *(const RefDiagSCMatrix &d, const RefSymmSCMatrix &s)
{
  RefSCMatrix ret(s.dim(), s.dim(), s.kit());
  int n = s.dim()->n();
  for (int i=0; i<n; i++) {
      for (int j=0; j<n; j++) {
          ret.set_element(i,j, d.get_element(i)*s.get_element(i,j));
        }
    }
  return ret;
}
}

static RefSymmSCMatrix
weight_matrix(const RefDiagSCMatrix &d, const RefSymmSCMatrix &s)
{
  RefSymmSCMatrix ret = s.clone();
  int n = s.dim()->n();
  for (int i=0; i<n; i++) {
      for (int j=0; j<=i; j++) {
          ret.set_element(i,j, s.get_element(i,j)
                               *d.get_element(i)*d.get_element(j));
        }
    }
  return ret;
}

static int
nnmb_atom(int z, int l)
{
  if (l==0) {
      if (z <= 2) return 1;
      else if (z <= 10) return 2;
      else if (z <= 18) return 3;
    }
  else if (l==1) {
      if (z <= 4) return 0;
      else if (z <= 12) return 1;
      else if (z <= 20) return 2;
    }
  else if (l==2) {
      if (z <= 20) return 0;
    }
  else if (l==3) {
      if (z <= 56) return 0;
    }
  else {
      return 0;
    }
  ExEnv::errn() << "NAO: z too big" << endl;
  abort();
  return 0;
}

static int
nnmb_all_atom(int z, int maxl)
{
  int ret = 0;
  for (int i=0; i<=maxl; i++) {
      ret += nnmb_atom(z,i) * (2*i+1);
    }
  return ret;
}

static RefSymmSCMatrix
mhalf(const RefSymmSCMatrix &S)
{
  RefSCDimension tdim = S.dim();
  Ref<SCMatrixKit> kit = S.kit();

  // find a symmetric orthogonalization transform
  RefSCMatrix trans(tdim,tdim,kit);
  RefDiagSCMatrix eigval(tdim,kit);

  S.diagonalize(eigval,trans);

  Ref<SCElementOp> squareroot = new SCElementSquareRoot;
  eigval.element_op(squareroot);

  Ref<SCElementOp> invert = new SCElementInvert(1.0e-12);
  eigval.element_op(invert);

  RefSymmSCMatrix OL(tdim,kit);
  OL.assign(0.0);
  // OL = trans * eigval * trans.t();
  OL.accumulate_transform(trans, eigval);
  return OL;
}

static void
delete_partition_info(int natom, int *maxam_on_atom,
                      int **nam_on_atom, int ***amoff_on_atom)
{
  int i, j;
  for (i=0; i<natom; i++) {
      for (j=0; j<=maxam_on_atom[i]; j++) {
          delete[] amoff_on_atom[i][j];
        }
      delete[] nam_on_atom[i];
      delete[] amoff_on_atom[i];
    }
  delete[] maxam_on_atom;
  delete[] nam_on_atom;
  delete[] amoff_on_atom;
}

#ifdef DEBUG
static double
ttrace(const RefSCMatrix &N,
       const RefSymmSCMatrix &P,
       const RefSymmSCMatrix &S)
{
  RefSCMatrix Nt = N.t();
  RefSymmSCMatrix Pt = P.clone();
  Pt.assign(0.0);
  Pt.accumulate_transform(Nt, P);
  RefSymmSCMatrix St = S.clone();
  St.assign(0.0);
  St.accumulate_transform(Nt, S);
  return (mhalf(St)*Pt*mhalf(St)).trace();
}

// for N giving an orthonormal basis
static double
ttrace(const RefSCMatrix &N,
       const RefSymmSCMatrix &P)
{
  RefSCMatrix Nt = N.t();
  RefSymmSCMatrix Pt = P.clone();
  Pt.assign(0.0);
  Pt.accumulate_transform(Nt, P);
  return Pt.trace();
}
#endif

static RefSCMatrix
assemble(const RefSCDimension dim,
         const RefSCMatrix &Nm, int *Nm_map,
         const RefSCMatrix &Nr1, int *Nr1_map,
         const RefSCMatrix &Nr2 = 0,  int *Nr2_map = 0)
{
  int nnmb = Nm.ncol();
  int nr1 = Nr1.ncol();
  int nr2 = (Nr2.null()?0:Nr2.ncol());
  int nb = dim.n();
  if (nb != nnmb + nr1 + nr2) {
      ExEnv::errn() << "assemble: dim mismatch" << endl;
      abort();
    }
  RefSCMatrix N(Nm.rowdim(), Nm.rowdim(), Nm.kit());
  // collect Nm, Nr1, and Nr2 back into N
  int i;
  for (i=0; i<nnmb; i++) {
      if (Nm_map[i] < 0 || Nm_map[i] >= nb) {
          ExEnv::errn() << "assemble: bad Nm_map" << endl;
          abort();
        }
      N.assign_column(Nm.get_column(i), Nm_map[i]);
    }
  for (i=0; i<nr1; i++) {
      if (Nr1_map[i] < 0 || Nr1_map[i] >= nb) {
          ExEnv::errn() << "assemble: bad Nr1_map" << endl;
          abort();
        }
      N.assign_column(Nr1.get_column(i), Nr1_map[i]);
    }
  for (i=0; i<nr2; i++) {
      if (Nr2_map[i] < 0 || Nr2_map[i] >= nb) {
          ExEnv::errn() << "assemble: bad Nr2_map" << endl;
          abort();
        }
      N.assign_column(Nr2.get_column(i), Nr2_map[i]);
    }
  return N;
}

// form symmetry average NAO for each atom
static void
form_nao(const RefSymmSCMatrix &P, const RefSymmSCMatrix &S,
         const RefSCMatrix &N, const RefDiagSCMatrix &W, int natom,
         int *maxam_on_atom, int **nam_on_atom, int ***amoff_on_atom,
         const Ref<SCMatrixKit>& kit)
{
  int i,j,k,l,m;

  N.assign(0.0);
  W.assign(0.0);

  for (i=0; i<natom; i++) {
      for (j=0; j<=maxam_on_atom[i]; j++) {
          int nfunc = 2*j + 1;
          double oonfunc = 1.0/nfunc;
          int nt = nam_on_atom[i][j];
          RefSCDimension tdim(new SCDimension(nt));
          RefSymmSCMatrix Pt(tdim, kit);
          RefSymmSCMatrix St(tdim, kit);
          Pt.assign(0.0);
          St.assign(0.0);
          for (k=0; k<nt; k++) {
              for (l=0; l<nt; l++) {
                  double Stmp = 0.0;
                  double Ptmp = 0.0;
                  for (m=0; m<nfunc; m++) {
                      int ii = amoff_on_atom[i][j][k] + m;
                      int jj = amoff_on_atom[i][j][l] + m;
                      Stmp += S.get_element(ii,jj);
                      Ptmp += P.get_element(ii,jj);
                    }
                  St.set_element(k,l,Stmp*oonfunc);
                  Pt.set_element(k,l,Ptmp*oonfunc);
                }
            }
          // find a symmetric orthogonalization transform
          RefSymmSCMatrix OL = mhalf(St);

          // transform Pt to the orthogonal basis
          RefSymmSCMatrix PtL(tdim,kit);
          PtL.assign(0.0);
          PtL.accumulate_transform(OL, Pt);

          // diagonalize PtL
          RefSCMatrix trans(tdim,tdim,kit);
          RefDiagSCMatrix eigval(tdim,kit);
          PtL.diagonalize(eigval, trans);

          // transform trans to the nonortho basis
          trans = OL * trans;

#         ifdef DEBUG
          eigval.print("eigval");
#         endif
          // fill in the elements of W
          for (k=0; k<nt; k++) {
              // the eigenvalues come out backwards: reverse them
              int krev = nt-k-1;
              double elem = eigval.get_element(krev);
              for (m=0; m<nfunc; m++) {
                  int ii = amoff_on_atom[i][j][k] + m;
#                 ifdef DEBUG
                  ExEnv::outn().form("W(%2d) = %12.8f\n", ii, elem);
#                 endif
                  W.set_element(ii, elem);
                }
            }

          // fill in the elements of N
          for (k=0; k<nt; k++) {
              for (l=0; l<nt; l++) {
                  // the eigenvalues come out backwards: reverse them
                  int lrev = nt-l-1;
                  double elem = trans.get_element(k,lrev);
                  for (m=0; m<nfunc; m++) {
                      int ii = amoff_on_atom[i][j][k] + m;
                      int jj = amoff_on_atom[i][j][l] + m;
                      N.set_element(ii,jj, elem);
                    }
                }
            }
        }
    }
}

// From "Natural Population Analysis", Alan E. Reed, Robert B. Weinstock,
// Frank Weinhold, JCP, 83 (1985), p 735.
RefSCMatrix
Wavefunction::nao(double *atom_charges)
{

  Ref<GaussianBasisSet> b = basis();
  Ref<PetiteList> pl = integral()->petite_list();

  // compute S, the ao basis overlap
  RefSymmSCMatrix blockedS = pl->to_AO_basis(overlap());
  RefSymmSCMatrix S
      = dynamic_cast<BlockedSymmSCMatrix*>(blockedS.pointer())->block(0);
  blockedS = 0;
# ifdef DEBUG
  S.print("S");
# endif

  // compute P, the ao basis density
  RefSymmSCMatrix P
      = dynamic_cast<BlockedSymmSCMatrix*>(ao_density().pointer())->block(0);

  // why?  good question.
  RefSymmSCMatrix Ptmp = P->clone();
  Ptmp.assign(0.0);
  Ptmp->accumulate_transform(S, P);
# ifdef DEBUG
  P.print("P");
  ExEnv::out0() << "nelec = " << (mhalf(S) * Ptmp * mhalf(S)).trace() << endl;
  ExEnv::out0() << "nelec(2) = " << (P * S).trace() << endl;
# endif
  P = Ptmp;
  Ptmp = 0;

  int i,j,k,l;
  int nb = b->nbasis();
  int nsh = b->nshell();
  int natom = molecule()->natom();

# ifdef DEBUG
  ExEnv::out0() << "nb = " << nb << endl;
  ExEnv::out0() << "nsh = " << nsh << endl;
  ExEnv::out0() << "natom = " << natom << endl;
# endif

  // Step 2a. Transform to solid harmonics.
  // -- for now program will abort if basis does not use only S.H and cart d.
  RefSCDimension aodim = P.dim();
  RefSCMatrix Tdfg(aodim, aodim, matrixkit());
  Tdfg->unit();
  for (i=0; i<nsh; i++) {
      const GaussianShell &shell = b->shell(i);
      int off = b->shell_to_function(i);
      for (j=0; j<shell.ncontraction(); j++) {
          if (shell.am(j) == 2 && ! shell.is_pure(j)) {
              for (k=0; k<6; k++) {
                  for (l=0; l<6; l++) {
                      Tdfg.set_element(off+k,off+l,0.0);
                    }
                }
              // this will put the s function first and the d second
              // first grab the s function
              SphericalTransformIter *sti;
              sti = integral()->new_spherical_transform_iter(2,0,0);
              for (sti->begin(); sti->ready(); sti->next()) {
                  Tdfg->set_element(off + sti->pureindex(),
                                    off + sti->cartindex(),
                                    sti->coef());
                }
              delete sti;
              // now for the pure d part of the cartesian d shell
              sti = integral()->new_spherical_transform_iter(2,0,2);
              for (sti->begin(); sti->ready(); sti->next()) {
                  Tdfg->set_element(off + sti->pureindex() + 1,
                                    off + sti->cartindex(),
                                    sti->coef());
                }
              delete sti;
            }
          else if (shell.am(j) > 2 && ! shell.is_pure(j)) {
              ExEnv::errn() << "NAOs can only be computed for puream if am > 2" << endl;
              abort();
            }
          off += shell.nfunction(j);
        }
    }

  // Tdfg should already be orthogonal, normalize them
//   RefSCMatrix Tdfgo = Tdfg*Tdfg.t();
//   RefDiagSCMatrix Tdfg_norm(Tdfg.rowdim(), matrixkit());
//   for (i=0; i<nb; i++) {
//       double o = Tdfgo.get_element(i,i);
//       Tdfg_norm.set_element(i,1.0/sqrt(o));
//     }
//   Tdfgo = 0;
//   Tdfg = Tdfg_norm * Tdfg;

# ifdef DEBUG
  Tdfg.print("Tdfg");
  (Tdfg.t() * Tdfg).print("Tdfg.t() * Tdfg");
  (Tdfg * Tdfg.t()).print("Tdfg * Tdfg.t()");
# endif

  RefSymmSCMatrix Pdfg(aodim, matrixkit());
  // Pdfp = Tdfp.t() * P * Tdfp
  Pdfg.assign(0.0); Pdfg.accumulate_transform(Tdfg, P);
  RefSymmSCMatrix Sdfg(aodim, matrixkit());
  // Sdfp = Tdfp.t() * S * Tdfp
  Sdfg.assign(0.0); Sdfg.accumulate_transform(Tdfg, S);
# ifdef DEBUG
  ExEnv::out0() << "nelec = " << (mhalf(Sdfg) * Pdfg * mhalf(Sdfg)).trace() << endl;
# endif

  // Step 2b. Partitioning and symmetry averaging of P and S
  // Partitioning:
  int *maxam_on_atom = new int[natom];
  int **nam_on_atom = new int*[natom];
  int ***amoff_on_atom = new int**[natom];
  int maxam = -1;
  for (i=0; i<natom; i++) {
      maxam_on_atom[i] = -1;
      for (j=0; j<b->nshell_on_center(i); j++) {
          GaussianShell &shell = b->shell(i,j);
          int maxam_on_shell = shell.max_angular_momentum();
          if (maxam_on_atom[i] < maxam_on_shell)
              maxam_on_atom[i] = maxam_on_shell;
        }
      if (maxam_on_atom[i] > maxam) maxam = maxam_on_atom[i];
      nam_on_atom[i] = new int[maxam_on_atom[i]+1];
      for (j=0; j<=maxam_on_atom[i]; j++) {
          nam_on_atom[i][j] = 0;
          for (k=0; k<b->nshell_on_center(i); k++) {
              GaussianShell &shell = b->shell(i,k);
              for (l=0; l<shell.ncontraction(); l++) {
                  if (shell.am(l) == j) nam_on_atom[i][j]++;
                  if (shell.am(l) == 2 && !shell.is_pure(l) && j == 0) {
                      // the s component of a cartesian d
                      nam_on_atom[i][0]++;
                    }
                }
            }
        }
      amoff_on_atom[i] = new int*[maxam_on_atom[i]+1];
      for (j=0; j<=maxam_on_atom[i]; j++) {
          amoff_on_atom[i][j] = new int[nam_on_atom[i][j]];
          int nam = 0;