csgrad.cc

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

                       biggest_ints_3a.indices(i)[2],
                       biggest_ints_3a.indices(i)[3],
                       biggest_ints_3a.val(i)
                       )
           << endl;
      }
    ExEnv::outn() << "biggest 3/4 transformed ints (in 4.)" << endl;
    for (i=0; i<biggest_ints_3.ncontrib(); i++) {
      ExEnv::outn() << scprintf("%3d %3d %3d %3d %16.12f",
                       biggest_ints_3.indices(i)[0],
                       biggest_ints_3.indices(i)[1],
                       biggest_ints_3.indices(i)[2],
                       biggest_ints_3.indices(i)[3],
                       biggest_ints_3.val(i)
                       )
           << endl;
      }
#endif

    if (restart_orbital_memgrp_ == 0) {
      if (biggest_coefs.ncontrib()) {
        ExEnv::out0() << endl << indent
             << "Largest first order coefficients (unique):"
             << endl;
        }
      for (i=0; i<biggest_coefs.ncontrib(); i++) {
        int i0 = biggest_coefs.indices(i)[0];
        int i1 = biggest_coefs.indices(i)[1];
        int i2 = biggest_coefs.indices(i)[2] + nocc;
        int i3 = biggest_coefs.indices(i)[3] + nocc;
        int spincase = biggest_coefs.indices(i)[4];
        ExEnv::out0() << indent
             << scprintf("  %2d %12.8f %2d %3s %2d %3s -> %2d %3s %2d %3s (%s)",
                         i+1, biggest_coefs.val(i),
                         symorb_num_[i0]+1,
                         ct.gamma(symorb_irrep_[i0]).symbol(),
                         symorb_num_[i1]+1,
                         ct.gamma(symorb_irrep_[i1]).symbol(),
                         symorb_num_[i2]+1,
                         ct.gamma(symorb_irrep_[i2]).symbol(),
                         symorb_num_[i3]+1,
                         ct.gamma(symorb_irrep_[i3]).symbol(),
                         (spincase==1111?"++++":"+-+-")
               )
             << endl;
        }
      }

    // Print out various energies etc.

    if (debug_) {
      ExEnv::out0() << indent << "Number of shell quartets for which AO integrals\n"
           << indent << "(or integral derivatives) would have been computed\n"
           << indent << "without bounds checking: "
           << npass*nshell*nshell*(nshell+1)*(nshell+1)/2
           << endl;

      ExEnv::out0() << indent << "Number of shell quartets for which AO integrals\n"
           << indent << "were computed: " << aoint_computed
           << endl;

      if (dograd) {
        ExEnv::out0() << indent
             << "Number of shell quartets for which AO integral derivatives\n"
             << indent << "were computed: " << aointder_computed
             << endl;
        }
      }

    ExEnv::out0()<<endl<<indent
        <<scprintf("RHF energy [au]:                   %17.12lf\n", escf);
    ExEnv::out0()<<indent
        <<scprintf("MP2 correlation energy [au]:       %17.12lf\n", ecorr_mp2);
    ExEnv::out0()<<indent
        <<scprintf("MP2 energy [au]:                   %17.12lf\n", emp2);
    ExEnv::out0().flush();
    }
  if (method_ && strcmp(method_,"mp")) {
    ExEnv::out0() << indent
         << "MBPT2: bad method for closed shell case: " << method_
         << ", using mp" << endl;
    }
  set_energy(emp2);
  set_actual_value_accuracy(reference_->actual_value_accuracy()
                            *ref_to_mp2_acc);

  RefSCDimension nocc_act_dim(new SCDimension(nocc_act,1));
  nocc_act_dim->blocks()->set_subdim(0, new SCDimension(nocc_act));
  RefSCDimension nvir_act_dim(new SCDimension(nvir_act,1));
  nvir_act_dim->blocks()->set_subdim(0, new SCDimension(nvir_act));
  RefSCDimension nocc_dim(new SCDimension(nocc,1));
  nocc_dim->blocks()->set_subdim(0, new SCDimension(nocc));
  RefSCDimension nvir_dim(new SCDimension(nvir,1));
  nvir_dim->blocks()->set_subdim(0, new SCDimension(nvir));
  RefSCDimension nbasis_dim = ao_dimension()->blocks()->subdim(0);
  RefSCDimension noso_dim(new SCDimension(noso,1));

  if (dograd || do_d1_) {
    msg_->sum(Laj,nvir*nocc);

    RefSCMatrix T1_mat(nocc_act_dim, nvir_act_dim, kit);
    // the elements of T1_mat are the single-substitution amplitudes
    BiggestContribs biggest_t1(2,10);
    // compute the S2 norm of Lee et al. (s2_diag)
    double s2_diag = 0.0;
    for (j=nfzc; j<nocc; j++) {
      laj_ptr = &Laj[j*nvir];
      for (a=0; a<nvir_act; a++) {
        tmpval = 0.5**laj_ptr++/(evals[a+nocc]-evals[j]);
        T1_mat.set_element(j-nfzc,a,tmpval);
        biggest_t1.insert(tmpval,j,a);
        s2_diag += tmpval*tmpval;
        }
      }
    s2_diag = sqrt(s2_diag/(2*nocc_act));
    // compute the T1 matrix 1-norm
    double t1onenorm = 0.0;
    Ref<SCElementSumAbs> sumabs = new SCElementSumAbs;
    Ref<SCElementOp> genop = sumabs.pointer();
    for (a=0; a < nvir_act; a++) {
      sumabs->init();
      T1_mat.get_column(a).element_op(genop);
      if (t1onenorm < sumabs->result()) t1onenorm = sumabs->result();
      }
    // compute the T1 matrix inf-norm
    double t1infnorm = 0.0;
    for (j=0; j < nocc_act; j++) {
      sumabs->init();
      T1_mat.get_row(j).element_op(genop);
      if (t1infnorm < sumabs->result()) t1infnorm = sumabs->result();
      }
    // compute the T1 matrix 2-norm ( = D1(MP2) )
    RefSymmSCMatrix D1_mat(nocc_act_dim,kit);
    D1_mat.assign(0.0);
    D1_mat.accumulate_symmetric_product(T1_mat);
    T1_mat = 0;
    double d1_diag = sqrt(D1_mat.eigvals().get_element(nocc_act-1));
    D1_mat = 0;
    // print the norms
    ExEnv::out0()<<endl;
    ExEnv::out0()
         <<indent<<scprintf("D1(MP2)                = %12.8f", d1_diag)
         <<endl
         <<indent<<scprintf("S2 matrix 1-norm       = %12.8f", t1onenorm)
         <<endl
         <<indent<<scprintf("S2 matrix inf-norm     = %12.8f", t1infnorm)
         <<endl
         <<indent<<scprintf("S2 diagnostic          = %12.8f", s2_diag)
         <<endl;
    if (biggest_t1.ncontrib()) {
      ExEnv::out0() << endl
           << indent << "Largest S2 values (unique determinants):" << endl;
      }
    for (i=0; i<biggest_t1.ncontrib(); i++) {
      int i0 = biggest_t1.indices(i)[0];
      int i1 = biggest_t1.indices(i)[1] + nocc;
      ExEnv::out0()
           << indent << scprintf("  %2d %12.8f %2d %3s -> %2d %3s",
                                 i+1, biggest_t1.val(i),
                                 symorb_num_[i0]+1,
                                 ct.gamma(symorb_irrep_[i0]).symbol(),
                                 symorb_num_[i1]+1,
                                 ct.gamma(symorb_irrep_[i1]).symbol())
           << endl;
      }

    } // if (dograd || do_d1_)

  for (i=0; i<thr_->nthread(); i++) {
      delete e12thread[i];
    }
  delete[] e12thread;

  // quit here if only the energy is needed
  if (!dograd) {
    delete[] tbints_; tbints_ = 0;
    if (do_d1_) delete[] Laj;
    delete[] scf_vector;
    delete[] scf_vector_dat;
    delete[] evals;
    tim_exit("mp2-mem");
    return;
    }

  // Accumulate intermediate gradients on node 0
  sum_gradients(msg_, ginter, natom, 3);

  // Add intermediate gradients to the gradient on node 0
  if (me==0)
    accum_gradients(gradient, ginter, natom, 3);

  // Print out contribution to the gradient from non-sep. 2PDM
  if (debug_) {
    print_natom_3(ginter,
              "Contribution to MP2 gradient from non-separable 2PDM [au]:");
    }

  ////////////////////////////////////////////////////////
  // Add contributions from all nodes to various matrices
  ////////////////////////////////////////////////////////
  tmpint = (nvir > nocc ? nvir:nocc);
  double *tmpmat = new double[tmpint*tmpint];
  msg_->sum(Pkj,nocc*(nocc+1)/2,tmpmat); // Pkj is now complete
  msg_->sum(Pab,nvir*(nvir+1)/2,tmpmat); // Pab is now complete
  msg_->sum(Wab,nvir*nvir,tmpmat);
  msg_->sum(Wkj,nocc*nocc,tmpmat);
  msg_->sum(Waj,nvir*nocc,tmpmat);
  if (do_d2_) {
    msg_->sum(d2occ_mat,nocc_act*(nocc_act+1)/2,tmpmat); 
    msg_->sum(d2vir_mat,nvir_act*(nvir_act+1)/2,tmpmat); 
    }
  delete[] tmpmat;

  // Compute D2 diagnostic (d2_diag) from matrices d2_occ_mat and d2_vir_mat
  if (do_d2_) {
    RefSymmSCMatrix D2occ_mat(nocc_act_dim, kit);
    RefSymmSCMatrix D2vir_mat(nvir_act_dim,kit);
    D2occ_mat->assign(d2occ_mat);
    D2vir_mat->assign(d2vir_mat);
    d2o = sqrt(D2occ_mat.eigvals().get_element(nocc_act-1));
    d2v = sqrt(D2vir_mat.eigvals().get_element(nvir_act-1));
    d2_diag = (d2o > d2v ? d2o:d2v);
    ExEnv::out0() << endl
         << indent <<  scprintf("D2(MP1) = %12.8f", d2_diag) << endl << endl;
    delete[] d2occ_mat;
    delete[] d2vir_mat;
  }

  // Finish computation of Wab
  tim_enter("Pab and Wab");
  pab_ptr = Pab;
  for (a=0; a<nvir_act; a++) {  // active-active part of Wab
    wba_ptr = &Wab[a];
    wab_ptr = &Wab[a*nvir];
    for (b=0; b<=a; b++) {
      if (a==b) {
        *wab_ptr -= evals[nocc+a]**pab_ptr;
        }
      else {
        *wab_ptr -= evals[nocc+a]**pab_ptr;
        *wba_ptr -= evals[nocc+b]**pab_ptr;
        } 
      pab_ptr++;
      wab_ptr++;
      wba_ptr += nvir;
      } // exit b loop
    }   // exit a loop
  for (a=0; a<nfzv; a++) {  // active-frozen part of Wab
    wba_ptr = &Wab[nvir_act+a];
    wab_ptr = &Wab[(nvir_act+a)*nvir];
    pab_ptr = &Pab[(nvir_act+a)*(nvir_act+a+1)/2];
    for (b=0; b<nvir_act; b++) {
      *wab_ptr -= evals[nocc+b]**pab_ptr;
      *wba_ptr -= evals[nocc+b]**pab_ptr;
      pab_ptr++;
      wab_ptr++;
      wba_ptr += nvir;
      } // exit b loop
    }   // exit a loop
  // Wab is now complete
  tim_exit("Pab and Wab");
  RefSCMatrix Wab_matrix(nvir_dim, nvir_dim, kit);
  Wab_matrix->assign(Wab); // Put elements of Wab into Wab_matrix
  free(Wab);

  // Update Wkj with contribution from Pkj
  tim_enter("Pkj and Wkj");
  pkj_ptr = Pkj;
  for (k=0; k<nocc; k++) {
    wjk_ptr = &Wkj[k];
    wkj_ptr = &Wkj[k*nocc];
    for (j=0; j<=k; j++) {
      if (k<nfzc && j<nfzc) {   // don't want both j and k frozen
        wkj_ptr++;
        wjk_ptr += nocc;
        pkj_ptr++;
        continue;
        }
      if (j==k) {
        *wkj_ptr++ -= evals[k]**pkj_ptr++;
        }
      else if (j<nfzc) {
        *wkj_ptr++ -= evals[k]**pkj_ptr;
        *wjk_ptr   -= evals[k]**pkj_ptr;
        pkj_ptr++;
        }
      else {
        *wkj_ptr++ -= evals[k]**pkj_ptr;
        *wjk_ptr   -= evals[j]**pkj_ptr;
        pkj_ptr++;
        }
      wjk_ptr += nocc;
      }  // exit j loop
    }    // exit k loop
  tim_exit("Pkj and Wkj");

  /////////////////////////////////
  // Finish the computation of Laj
  /////////////////////////////////

  tim_enter("Laj");

  RefSCMatrix Cv(nbasis_dim, nvir_dim, kit); // virtual block of scf_vector
  RefSCMatrix Co(nbasis_dim, nocc_dim, kit); // occupied block of scf_vector
  for (p=0; p<nbasis; p++) {
    c_pq = scf_vector[p];
    for (q=0; q<noso; q++) {
      if (q<nocc) Co->set_element(p, q, *c_pq++);
      else Cv->set_element(p, q-nocc, *c_pq++);
      }
    }

  // Compute the density-like Dmat_matrix
  // (Cv*Pab_matrix*Cv.t() + Co*Pkj_matrix*Co.t())
  RefSymmSCMatrix Pab_matrix(nvir_dim,kit);
  RefSymmSCMatrix Pkj_matrix(nocc_dim,kit);
  RefSymmSCMatrix Dmat_matrix(nbasis_dim,kit);
  Pab_matrix->assign(Pab); // fill in elements of Pab_matrix from Pab
  free(Pab);
  Pkj_matrix->assign(Pkj); // fill in elements of Pkj_matrix from Pkj
  free(Pkj);
  Dmat_matrix.assign(0.0);
  Dmat_matrix.accumulate_transform(Cv,Pab_matrix);
  Dmat_matrix.accumulate_transform(Co,Pkj_matrix);
  // We now have the density-like matrix Dmat_matrix

  // Compute the G matrix

  RefSymmSCMatrix Gmat(nbasis_dim,kit);
  init_cs_gmat();
  tim_enter("make_gmat for Laj");
  make_cs_gmat_new(Gmat, Dmat_matrix); 
  if (debug_ > 1) {
    Dmat_matrix.print("Dmat");
    Gmat.print("Gmat");
    }
  tim_exit("make_gmat for Laj");

  // Finish computation of Laj
  RefSCMatrix Laj_matrix(nocc_dim,nvir_dim,kit); // elements are ordered as j*nvir+a
  Laj_matrix->assign(Laj);
  if (debug_ > 1) Laj_matrix->print("Laj (first bit)");
  Laj_matrix = Laj_matrix - 2*Co.t()*Gmat*Cv;
  if (debug_ > 1) Laj_matrix->print("Laj (all of it)");
  Laj_matrix->convert(Laj);  // Put new Laj_matrix elements into Laj

  tim_exit("Laj");

  //////////////////////////////////////
  // Computation of Laj is now complete
  //////////////////////////////////////

  ////////////////////////////
  // Solve the CPHF equations
  ////////////////////////////
  RefSCMatrix Paj_matrix(nvir_dim, nocc_dim, kit);
  tim_enter("cphf");
  cs_cphf(scf_vector, Laj, evals, Paj_matrix);
  tim_exit("cphf");

  free(Laj);

  // Finish computation of Waj
  for (a=0; a<nvir; a++) {
    waj_ptr = &Waj[a];
    for (j=0; j<nocc; j++) {
      *waj_ptr -= evals[j]*Paj_matrix->get_element(a,j);
      waj_ptr += nvir;
      }
    }
  // Waj is now complete
  RefSCMatrix Waj_matrix(nocc_dim, nvir_dim, kit);
  Waj_matrix->assign(Waj); // Put elements of Waj into Waj_matrix
  // NB. Waj_matrix elements are ordered as j*nvir+a
  free(Waj);


  // Finish computation of Wkj
  tim_enter("Pkj and Wkj");
  // Compute Dmat_matrix =
  // Co*Pkj_matrix*Co.t() + Co*Paj_matrix.t()*Cv.t()
  // + Cv*Paj_matrix*Co.t() + Cv*Pab_matrix*Cv.t();
  Dmat_matrix.assign(0.0);
  Dmat_matrix.accumulate_symmetric_sum(Cv*Paj_matrix*Co.t());
  Dmat_matrix.accumulate_transform(Co,Pkj_matrix);
  Dmat_matrix.accumulate_transform(Cv,Pab_matrix);
  tim_enter("make_gmat for Wkj");
  make_cs_gmat_new(Gmat, Dmat_matrix);
  tim_exit("make_gmat for Wkj");
  done_cs_gmat();
  for (i=0; i<thr_->nthread(); i++) tbints_[i] = 0;
  delete[] tbints_