//// nao.cc//// Copyright (C) 1997 Limit Point Systems, Inc.//// Author: Curtis Janssen <cljanss@limitpt.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.//#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 DEBUGnamespace sc {static RefSCMatrixoperator *(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 RefSymmSCMatrixweight_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 intnnmb_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 intnnmb_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 RefSymmSCMatrixmhalf(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 voiddelete_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 DEBUGstatic doublettrace(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 basisstatic doublettrace(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();}#endifstatic RefSCMatrixassemble(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 atomstatic voidform_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.RefSCMatrixWavefunction::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;