nao.cc

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

          for (k=0; k<b->nshell_on_center(i); k++) {
              GaussianShell &shell = b->shell(i,k);
              int function_offset
                  = b->shell_to_function(b->shell_on_center(i,k));
              int conoffset = 0;
              for (l=0; l<shell.ncontraction(); l++) {
                  if (shell.am(l) == j) {
                      amoff_on_atom[i][j][nam]
                          = function_offset + conoffset;
                      if (j == 2 && !shell.is_pure(l)) {
                          // the pure d part of a cartesian d shell is offset
                          amoff_on_atom[i][j][nam]++;
                        }
                      nam++;
                    }
                  if (shell.am(l) == 2 && (!shell.is_pure(l)) && j == 0) {
                      // the s component of a cartesian d
                      amoff_on_atom[i][j][nam]
                          = function_offset + conoffset;
                      nam++;
                    }
                  conoffset += shell.nfunction(l);
                }
            }
        }
    }

# ifdef DEBUG
  ExEnv::out0() << indent << "Basis set partitioning:" << endl;
  ExEnv::out0() << incindent;
  for (i=0; i<natom; i++) {
      ExEnv::out0() << indent <<  "atom " << i
           << " maxam = " << maxam_on_atom[i] << endl;
      ExEnv::out0() << incindent;
      for (j=0; j<=maxam_on_atom[i]; j++) {
          ExEnv::out0() << indent <<  "am = " << j
               << " n = " << nam_on_atom[i][j] << endl;
          ExEnv::out0() << incindent;
          ExEnv::out0() << indent << "offsets =";
          for (k=0; k<nam_on_atom[i][j]; k++) {
              ExEnv::out0() << " " << amoff_on_atom[i][j][k];
            }
          ExEnv::out0() << endl;
          ExEnv::out0() << decindent;
        }
      ExEnv::out0() << decindent;
    }
  ExEnv::out0() << decindent;
# endif

  // Symmetry averaging and Step 2c: Formation of pre-NAO's
  RefSCMatrix N(aodim, aodim, matrixkit());
  RefDiagSCMatrix W(aodim, matrixkit());
  form_nao(Pdfg, Sdfg, N, W, natom, maxam_on_atom, nam_on_atom, amoff_on_atom,
           matrixkit());
# ifdef DEBUG
  N.print("N");
  W.print("W");
  ExEnv::out0() << "nelec = " << ttrace(N, Pdfg, Sdfg) << endl;
# endif

  // Step 3a: selection of NMB orbitals

  // count the size of nmb
  int nnmb = 0;
  for (i=0; i<natom; i++) {
      nnmb += nnmb_all_atom(molecule()->Z(i),
                            maxam_on_atom[i]);
    }
  int nnrb = nb - nnmb;

# ifdef DEBUG
  ExEnv::out0() << "nnmb = " << nnmb << endl;
  ExEnv::out0() << "nnrb = " << nnrb << endl;
# endif

  RefSCDimension nmbdim = new SCDimension(nnmb);
  RefSCDimension nrbdim = new SCDimension(nnrb);

  // split N into the nmb and nrb parts
  RefSCMatrix Nm(aodim, nmbdim, matrixkit());
  RefSCMatrix Nr(aodim, nrbdim, matrixkit());
  RefDiagSCMatrix Wm(nmbdim, matrixkit());
  RefDiagSCMatrix Wr(nrbdim, matrixkit());
  int *Nm_map = new int[nnmb];
  int *Nr_map = new int[nnrb];

  int im = 0;
  int ir = 0;
  for (i=0; i<natom; i++) {
      int z = molecule()->Z(i);
      for (j=0; j<=maxam_on_atom[i]; j++) {
          int nnmb_zj = nnmb_atom(z,j);
          int nt = nam_on_atom[i][j];
          for (k=0; k<nt; k++) {
              int iN = amoff_on_atom[i][j][k];
              if (k<nnmb_zj) {
                  for (l=0; l<(2*j+1); l++) {
                      Nm_map[im] = iN;
                      Wm.set_element(im, W.get_element(iN));
                      Nm.assign_column(N.get_column(iN++),im++);
                    }
                }
              else {
                  for (l=0; l<(2*j+1); l++) {
                      Nr_map[ir] = iN;
                      Wr.set_element(ir, W.get_element(iN));
                      Nr.assign_column(N.get_column(iN++),ir++);
                    }
                }
            }
        }
    }
# ifdef DEBUG
  ExEnv::out0() << "Nmmap:"; for (i=0;i<nnmb;i++) ExEnv::out0()<<" "<<Nm_map[i]; ExEnv::out0()<<endl;
  ExEnv::out0() << "Nrmap:"; for (i=0;i<nnrb;i++) ExEnv::out0()<<" "<<Nr_map[i]; ExEnv::out0()<<endl;
  Wm.print("Wm");
  Wr.print("Wr");
  Nm.print("Nm");
  Nr.print("Nr");
  (Nm.t() * Sdfg * Nr).print("3a Smr");
  ExEnv::out0() << "nelec = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr,Nr_map), Pdfg, Sdfg) << endl;
# endif

  // Step 3b: Schmidt interatomic orthogonalization of NRB to NMB orbs

  // orthogonalize the NMB orbs (temporarily, to project them out of NRB)
  int ii=0;
  for (i=0; i<nnmb; i++,ii++) {
      N.assign_column(Nm.get_column(i),ii);
    }
  for (i=0; i<nnrb; i++,ii++) {
      N.assign_column(Nr.get_column(i),ii);
    }
  N->schmidt_orthog(Sdfg.pointer(),nnmb);

  RefSCMatrix Nmo = Nm.clone();
  for (i=0; i<nnmb; i++) {
      Nmo.assign_column(N.get_column(i),i);
    }
  RefSCMatrix OSmr = Nmo.t() * Sdfg * Nr;
  OSmr.scale(-1.0);
  Nr.accumulate(Nmo * OSmr);
# ifdef DEBUG
  OSmr.print("OSmr");
  Nmo.print("Nmo = Nm after temporay orthog");
  Nr.print("Nr after orthogonalization to NMB");
  (Nm.t() * Sdfg * Nr).print("3b Smr");
# endif
  Nmo = 0;

  // Step 3c: Restoration of natural character of the NRB
  // Partitioning:
  int *r_maxam_on_atom = new int[natom];
  int **r_nam_on_atom = new int*[natom];
  int ***r_amoff_on_atom = new int**[natom];
  int r_offset = 0;
  for (i=0; i<natom; i++) {
      int z = molecule()->Z(i);
      r_maxam_on_atom[i] = maxam_on_atom[i];
      r_nam_on_atom[i] = new int[r_maxam_on_atom[i]+1];
      for (j=0; j<=r_maxam_on_atom[i]; j++) {
          r_nam_on_atom[i][j] = nam_on_atom[i][j] - nnmb_atom(z,j);
          if (r_nam_on_atom[i][j] < 0) {
              ExEnv::errn() << "NAO: < 0 rydberg orbitals of a given type" << endl;
              abort();
            }
        }
      r_amoff_on_atom[i] = new int*[r_maxam_on_atom[i]+1];
      for (j=0; j<=r_maxam_on_atom[i]; j++) {
          r_amoff_on_atom[i][j] = new int[r_nam_on_atom[i][j]];
          for (k=0; k<r_nam_on_atom[i][j]; k++) {
              r_amoff_on_atom[i][j][k] = r_offset;
              r_offset += 2*j + 1;
            }
        }
    }
  RefSymmSCMatrix Pr(nrbdim, matrixkit());
  // Pr = Nr.t() * Tdfg.t() * P * Tdfg * Nr;
  Pr.assign(0.0); Pr.accumulate_transform(Nr.t(), Pdfg);
  RefSymmSCMatrix Sr(nrbdim, matrixkit());
  // Sr = Nr.t() * Tdfg.t() * S * Tdfg * Nr;
  Sr.assign(0.0); Sr.accumulate_transform(Nr.t(), Sdfg);

  // Symmetry averaging and restoration of natural character of NRB
  RefSCMatrix Nrr(nrbdim, nrbdim, matrixkit());
  form_nao(Pr, Sr, Nrr, Wr,
           natom, r_maxam_on_atom, r_nam_on_atom, r_amoff_on_atom,
           matrixkit());
  Nr = Nr * Nrr;
  // these are out-of-date
  Pr = 0; Sr = 0;
# ifdef DEBUG
  Wr.print("Wr after restoring natural character");
  Nr.print("Nr after restoring natural character");
  (Nm.t() * Sdfg * Nr).print("3c Smr");
# endif

  // Step 4a: Weighted interatomic orthogonalization
  // nmb part of OW
  RefSymmSCMatrix Sm(nmbdim, matrixkit());
  Sm.assign(0.0); Sm.accumulate_transform(Nm.t(), Sdfg);
  RefSymmSCMatrix SWm = weight_matrix(Wm, Sm);
  RefSCMatrix OWm = Wm * mhalf(SWm);
# ifdef DEBUG
  Sm.print("Sm before 4a");
  OWm.print("OWm");
  (OWm.t() * Sm * OWm).print("Sm after 4a");
  ExEnv::out0() << "nelec = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr,Nr_map), Pdfg, Sdfg) << endl;
# endif

  // put OWm into Nm
  Nm = Nm * OWm;

# ifdef DEBUG
  Nm.print("Nm after interatomic orthog");
  (Nm.t() * Sdfg * Nr).print("4a Smr before r orthog");
  ExEnv::out0() << "nelec = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr,Nr_map), Pdfg, Sdfg)
       << endl;
# endif

  // nrb part of OW
  // based on Wr, r is split into r1 and r2
  double tw = 1.0e-4; // the tolerance used for the split
  int nr1 = 0;
  int nr2 = 0;
  for (i=0; i<nnrb; i++) {
      if (fabs(Wr.get_element(i)) >= tw) nr1++;
      else nr2++;
    }
  RefSCDimension r1dim(new SCDimension(nr1));
  RefSCDimension r2dim(new SCDimension(nr2));
  RefSCMatrix Nr1(aodim, r1dim, matrixkit());
  RefSCMatrix Nr2(aodim, r2dim, matrixkit());
  RefDiagSCMatrix Wr1(r1dim, matrixkit());
  int *Nr1_map = new int[nr1];
  int *Nr2_map = new int[nr2];
  int ir1 = 0;
  int ir2 = 0;
  for (i=0; i<nnrb; i++) {
      if (fabs(Wr.get_element(i)) >= tw) {
          Nr1_map[ir1] = Nr_map[i];
          Wr1.set_element(ir1, Wr.get_element(i));
          Nr1.assign_column(Nr.get_column(i),ir1++);
        }
      else {
          Nr2_map[ir2] = Nr_map[i];
          Nr2.assign_column(Nr.get_column(i),ir2++);
        }
    }
# ifdef DEBUG
  ExEnv::out0() << "Nr1map:"; for (i=0;i<nr1;i++) ExEnv::out0()<<" "<<Nr1_map[i]; ExEnv::out0()<<endl;
  ExEnv::out0() << "Nr2map:"; for (i=0;i<nr2;i++) ExEnv::out0()<<" "<<Nr2_map[i]; ExEnv::out0()<<endl;
  Nr1.print("Nr1");
  Nr2.print("Nr2");
  (Nm.t() * Sdfg * Nr1).print("4a Smr1 before r orthog");
  (Nm.t() * Sdfg * Nr2).print("4a Smr2 before r orthog");
  ExEnv::out0() << "nelec = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
       << endl;
# endif

  // Schmidt orthogonalization of r2 to r1
  // Collect Nr together again (but in the order: r1, r2)
  ii=0;
  for (i=0; i<nr1; i++,ii++) {
      Nr.assign_column(Nr1.get_column(i),ii);
    }
  for (i=0; i<nr2; i++,ii++) {
      Nr.assign_column(Nr2.get_column(i),ii);
    }
  Nr->schmidt_orthog(Sdfg.pointer(),nr1);
  RefSCMatrix Nr1o = Nr1.copy();
  for (i=0; i<nr1; i++) {
      Nr1o.assign_column(Nr.get_column(i), i);
    }
  RefSCMatrix Or1r2 = Nr1o.t() * Sdfg * Nr2;
  Or1r2.scale(-1.0);
# ifdef DEBUG
  (Nm.t() * Sdfg * Nr2).print("4a Smr2 before orthog of r2 to r1");
  (Nr1.t() * Sdfg * Nr2).print("4a Sr1r2 before orthog of r2 to r1");
# endif
  Nr2.accumulate(Nr1o * Or1r2);
# ifdef DEBUG
  Nr2.print("Nr2 after orthogonalization to r1");
  (Nm.t() * Sdfg * Nr2).print("4a Smr2 after orthog of r2 to r1");
  (Nr1.t() * Sdfg * Nr2).print("4a Sr1r2 after orthog of r2 to r1");
  ExEnv::out0() << "nelec = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
       << endl;
# endif

  // weighted symmetric orthog of r1
  RefSymmSCMatrix Sr1(r1dim, matrixkit());
  Sr1.assign(0.0); Sr1.accumulate_transform(Nr1.t(), Sdfg);
  RefSymmSCMatrix SWr1 = weight_matrix(Wr1, Sr1);
  RefSCMatrix OWr1 = Wr1 * mhalf(SWr1);
# ifdef DEBUG
  OWr1.print("OWr1");
  (Nr1.t() * Sdfg * Nr1).print("Nr1.t() * Sdfg * Nr1");
# endif
  // Put OWr1 into Nr1
  Nr1 = Nr1 * OWr1;
# ifdef DEBUG
  Nr1.print("Nr1 after weighted symmetric orthogonalization");
  ExEnv::out0() << "nelec = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
       << endl;
# endif

  // symmetric orthog of r1
  RefSymmSCMatrix Sr2(r2dim, matrixkit());
  Sr2.assign(0.0); Sr2.accumulate_transform(Nr2.t(), Sdfg);
  RefSymmSCMatrix OWr2 = mhalf(Sr2);
# ifdef DEBUG
  OWr2.print("OWr2");
# endif

  // Put OWr2 into Nr2
  Nr2 = Nr2 * OWr2;
# ifdef DEBUG
  Nr2.print("Nr2 after weighted symmetric orthogonalization");
  ExEnv::out0() << "nelec = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg, Sdfg)
       << endl;
  ExEnv::out0() << "nelec(o) = "
       << ttrace(assemble(aodim,Nm,Nm_map,Nr1,Nr1_map,Nr2,Nr2_map), Pdfg)
       << endl;
# endif

  // Step 4b. restoration of the natural character of the naos

  // collect Nm, Nr1, and Nr2 back into N
  N = assemble(aodim, Nm,Nm_map, Nr1,Nr1_map, Nr2,Nr2_map);
# ifdef DEBUG
  N.print("N after 4a");
# endif

  // compute the density and overlap in the current basis
  // N currently has the entire transform, starting from the dfg basis
  P.assign(0.0);
  P.accumulate_transform(N.t(), Pdfg);
  S.assign(0.0);
  S.accumulate_transform(N.t(), Sdfg);
# ifdef DEBUG
  P.print("P after 4a");
  S.print("S after 4a");
  (Nm.t() * Sdfg * Nm).print("4a Sm");
  (Nr1.t() * Sdfg * Nr1).print("4a Sr1");
  (Nr2.t() * Sdfg * Nr2).print("4a Sr2");
  (Nm.t() * Sdfg * Nr1).print("4a Smr1");
  (Nm.t() * Sdfg * Nr2).print("4a Smr2");
  (Nr1.t() * Sdfg * Nr2).print("4a Sr1r2");
# endif

  RefSCMatrix Nred(aodim, aodim, matrixkit());
  form_nao(P, S, Nred, W, natom, maxam_on_atom, nam_on_atom, amoff_on_atom,
           matrixkit());
  N = N * Nred;

  RefSymmSCMatrix Pfinal(aodim, matrixkit());
  Pfinal.assign(0.0);
  Pfinal.accumulate_transform(N.t(), Pdfg);
# ifdef DEBUG
  Nred.print("Nred");
  N.print("N after 4b");
  ExEnv::out0() << "nelec = " << ttrace(N, Pdfg, Sdfg) << endl;
  ExEnv::out0() << "nelec(o) = " << ttrace(N, Pdfg) << endl;
  Pfinal.print("final P");
  (N.t() * Sdfg * N).print("final S");
  ExEnv::out0().form("nelec = trace(final P) = %14.8f", (N.t() * Pdfg * N).trace());

  (mhalf(Sdfg) * Pdfg * mhalf(Sdfg)).print("P in symm orth basis");
# endif

# ifdef DEBUG
  ExEnv::out0() << "nb   = " << nb << endl;
  ExEnv::out0() << "nnmb = " << nnmb << endl;
  ExEnv::out0() << "nnrb = " << nnrb << endl;
  ExEnv::out0() << "nr1  = " << nr1 << endl;
  ExEnv::out0() << "nr2  = " << nr2 << endl;
# endif

  ExEnv::out0() << indent << "Natural Population Analysis:" << endl;
  ExEnv::out0() << incindent;
  ExEnv::out0() << indent << " n   atom    charge ";
  for (i=0; i<=maxam; i++) {
      const char *am = "SPDFGH?";
      int index;
      if (i>6) index = 6;
      else index = i;
      ExEnv::out0() << "    ne(" << am[index] << ") ";
    }
  ExEnv::out0() << endl;
  for (i=0; i<natom; i++) {
      double e = 0.0;
      for (j=0; j<=maxam_on_atom[i]; j++) {
          for (k=0; k<nam_on_atom[i][j]; k++) {
              for (l=0; l<(2*j+1); l++) {
                  e += Pfinal.get_element(amoff_on_atom[i][j][k] + l,
                                          amoff_on_atom[i][j][k] + l);
                }
            }
        }
      ExEnv::out0() << indent
           << scprintf("%3d   %2s   % 8.6f",i + 1,
                       AtomInfo::symbol(molecule()->Z(i)),
                       double(molecule()->Z(i)) - e);
      if (atom_charges) {
          atom_charges[i] = molecule()->Z(i) - e;
        }
      for (j=0; j<=maxam_on_atom[i]; j++) {
          e = 0.0;
          for (k=0; k<nam_on_atom[i][j]; k++) {
              for (l=0; l<(2*j+1); l++) {
                  e += Pfinal.get_element(amoff_on_atom[i][j][k] + l,
                                          amoff_on_atom[i][j][k] + l);
                }
            }
          ExEnv::out0() << scprintf(" % 8.6f",e);
        }
      ExEnv::out0() << endl;
    }
  ExEnv::out0() << endl;
  ExEnv::out0() << decindent;

  delete[] Nm_map;
  delete[] Nr_map;
  delete[] Nr1_map;
  delete[] Nr2_map;
  delete_partition_info(natom,maxam_on_atom,nam_on_atom,amoff_on_atom);
  delete_partition_info(natom,r_maxam_on_atom,r_nam_on_atom,r_amoff_on_atom);

  return N;
}

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

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: