mp2r12_energy.cc

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

  delete[] er12_aa_vec;

  //
  // Alpha-beta pairs
  //

  RefSCMatrix Bab_ij = Bab.clone();
  // In MP2-R12/A the B matrix is the same for all pairs
  if (stdapprox_ == LinearR12::StdApprox_A) {
    Bab_ij.assign(Bab);
    if (debug_ > 1)
      Bab_ij.print("Inverse alpha-beta MP2-R12/A B matrix");
    Bab_ij->gen_invert_this();
  }

  ij=0;
  for(int i=0; i<nocc_act; i++)
    for(int j=0; j<nocc_act; j++, ij++) {

      if (ij%ntasks != me)
        continue;

      RefSCVector Vab_ij = Vab.get_column(ij);

      // In MP2-R12/A' matrices B are pair-specific:
      // Form B(ij)kl,ow = Bkl,ow + 1/2(ek + el + eo + ew - 2ei - 2ej)Xkl,ow
      if (stdapprox_ == LinearR12::StdApprox_Ap) {
        Bab_ij.assign(Bab);
        int kl=0;
        for(int k=0; k<nocc_act; k++)
          for(int l=0; l<nocc_act; l++, kl++) {
            int ow=0;
            for(int o=0; o<nocc_act; o++)
              for(int w=0; w<nocc_act; w++, ow++) {
                double fx = 0.5 * (evals[k] + evals[l] + evals[o] + evals[w] - 2.0*evals[i] - 2.0*evals[j]) *
                Xab.get_element(kl,ow);
                Bab_ij.accumulate_element(kl,ow,fx);
              }
          }
        if (debug_ > 1)
	    Bab_ij.print("Alpha-beta MP2-R12/A' B matrix");

        Bab_ij->gen_invert_this();

        if (debug_ > 1)
          Bab_ij.print("Inverse alpha-beta MP2-R12/A' B matrix");
      }

      double eab_ij = -1.0*Vab_ij.dot(Bab_ij * Vab_ij);
      er12_ab_vec[ij] = eab_ij;
    }
  Bab_ij=0;
  msg->sum(er12_ab_vec,nab,0,-1);
  er12_ab_->assign(er12_ab_vec);
  emp2r12_ab_->assign(emp2_ab);
  emp2r12_ab_->accumulate(er12_ab_);
  delete[] er12_ab_vec;

  delete[] evals;

  evaluated_ = true;
  
  return;
}

void MP2R12Energy::print(std::ostream& so) const
{
} 

void MP2R12Energy::print_pair_energies(bool spinadapted, std::ostream& so)
{
  compute();
  
  char* SA_str;
  switch (stdapprox_) {
    case LinearR12::StdApprox_A:
      SA_str = strdup("A");
      break;

    case LinearR12::StdApprox_Ap:
      SA_str = strdup("A'");
      break;

    case LinearR12::StdApprox_B:
      SA_str = strdup("B");
      break;

    default:
      throw std::runtime_error("MP2R12Energy::print_pair_energies -- stdapprox_ is not valid");
  }

  Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
  int nocc_act = r12info->nocc_act();
  double escf = r12info->ref()->energy();

  double emp2tot_aa = 0.0;
  double emp2tot_ab = 0.0;
  double er12tot_aa = 0.0;
  double er12tot_ab = 0.0;
  double emp2tot_0 = 0.0;
  double emp2tot_1 = 0.0;
  double er12tot_0 = 0.0;
  double er12tot_1 = 0.0;

  RefSCVector emp2_aa = r12eval_->emp2_aa();
  RefSCVector emp2_ab = r12eval_->emp2_ab();

  /*---------------------------------------
    Spin-adapt pair energies, if necessary
   ---------------------------------------*/
  if (!spinadapted) {

    so << endl << indent << "Alpha-alpha MBPT2-R12/" << SA_str << " pair energies:" << endl;
    so << indent << scprintf("    i       j        mp2(ij)        r12(ij)      mp2-r12(ij)") << endl;
    so << indent << scprintf("  -----   -----   ------------   ------------   ------------") << endl;
    for(int i=0,ij=0;i<nocc_act;i++)
      for(int j=0;j<i;j++,ij++) {
        so << indent << scprintf("  %3d     %3d     %12.9lf   %12.9lf   %12.9lf",i+1,j+1,
                                 emp2_aa->get_element(ij),
                                 er12_aa_->get_element(ij),
                                 emp2r12_aa_->get_element(ij)) << endl;
      }
      
    so << endl << indent << "Alpha-beta MBPT2-R12/" << SA_str << " pair energies:" << endl;
    so << indent << scprintf("    i       j        mp2(ij)        r12(ij)      mp2-r12(ij)") << endl;
    so << indent << scprintf("  -----   -----   ------------   ------------   ------------") << endl;
    for(int i=0,ij=0;i<nocc_act;i++)
      for(int j=0;j<nocc_act;j++,ij++) {
        so << indent << scprintf("  %3d     %3d     %12.9lf   %12.9lf   %12.9lf",i+1,j+1,
                                 emp2_ab->get_element(ij),
                                 er12_ab_->get_element(ij),
                                 emp2r12_ab_->get_element(ij)) << endl;
      }

  }
  else {

    Ref<SCMatrixKit> localkit = er12_aa_.kit();
    RefSCVector emp2r12_0 = localkit->vector(r12eval_->dim_s());
    RefSCVector emp2r12_1 = localkit->vector(r12eval_->dim_t());
    RefSCVector emp2_0 = localkit->vector(r12eval_->dim_s());
    RefSCVector emp2_1 = localkit->vector(r12eval_->dim_t());
    RefSCVector er12_0 = localkit->vector(r12eval_->dim_s());
    RefSCVector er12_1 = localkit->vector(r12eval_->dim_t());

    // Triplet pairs are easy
    emp2r12_1->assign(emp2r12_aa_);
    emp2r12_1->scale(1.5);
    emp2_1->assign(emp2_aa);
    emp2_1->scale(1.5);
    er12_1->assign(er12_aa_);
    er12_1->scale(1.5);

    // Singlet pairs are a bit trickier
    int ij_s = 0;
    for(int i=0; i<nocc_act; i++)
      for(int j=0; j<=i; j++, ij_s++) {
        int ij_ab = i*nocc_act + j;
        int ij_aa = i*(i-1)/2 + j;
        double eab, eaa, e_s;

        eab = emp2r12_ab_->get_element(ij_ab);
        if (i != j)
          eaa = emp2r12_aa_->get_element(ij_aa);
        else
          eaa = 0.0;
        e_s = (i != j ? 2.0 : 1.0) * eab - 0.5 * eaa;
        emp2r12_0->set_element(ij_s,e_s);

        eab = emp2_ab->get_element(ij_ab);
        if (i != j)
          eaa = emp2_aa->get_element(ij_aa);
        else
          eaa = 0.0;
        e_s = (i != j ? 2.0 : 1.0) * eab - 0.5 * eaa;
        emp2_0->set_element(ij_s,e_s);

        eab = er12_ab_->get_element(ij_ab);
        if (i != j)
          eaa = er12_aa_->get_element(ij_aa);
        else
          eaa = 0.0;
        e_s = (i != j ? 2.0 : 1.0) * eab - 0.5 * eaa;
        er12_0->set_element(ij_s,e_s);
      }

    // compute total singlet and triplet energies
    RefSCVector unit_0 = localkit->vector(r12eval_->dim_s());
    RefSCVector unit_1 = localkit->vector(r12eval_->dim_t());
    unit_0->assign(1.0);
    unit_1->assign(1.0);
    emp2tot_0 = emp2_0.dot(unit_0);
    emp2tot_1 = emp2_1.dot(unit_1);
    er12tot_0 = er12_0.dot(unit_0);
    er12tot_1 = er12_1.dot(unit_1);

    so << endl << indent << "Singlet MBPT2-R12/" << SA_str << " pair energies:" << endl;
    so << indent << scprintf("    i       j        mp2(ij)        r12(ij)      mp2-r12(ij)") << endl;
    so << indent << scprintf("  -----   -----   ------------   ------------   ------------") << endl;
    for(int i=0,ij=0;i<nocc_act;i++)
      for(int j=0;j<=i;j++,ij++) {
        so << indent << scprintf("  %3d     %3d     %12.9lf   %12.9lf   %12.9lf",i+1,j+1,
                                 emp2_0->get_element(ij),
                                 er12_0->get_element(ij),
                                 emp2r12_0->get_element(ij)) << endl;
      }
      
    so << endl << indent << "Triplet MBPT2-R12/" << SA_str << " pair energies:" << endl;
    so << indent << scprintf("    i       j        mp2(ij)        r12(ij)      mp2-r12(ij)") << endl;
    so << indent << scprintf("  -----   -----   ------------   ------------   ------------") << endl;
    for(int i=0,ij=0;i<nocc_act;i++)
      for(int j=0;j<i;j++,ij++) {
        so << indent << scprintf("  %3d     %3d     %12.9lf   %12.9lf   %12.9lf",i+1,j+1,
                                 emp2_1->get_element(ij),
                                 er12_1->get_element(ij),
                                 emp2r12_1->get_element(ij)) << endl;
      }

  }


  double mp2_corr_energy_ = emp2tot_aa_() + emp2tot_ab_();
  double r12_corr_energy_ = er12tot_aa_() + er12tot_ab_();
  
  ///////////////////////////////////////////////////////////////
  // The computation of the MP2 energy is now complete on each
  // node;
  ///////////////////////////////////////////////////////////////

  if (spinadapted) {
    so <<endl<<indent
      <<scprintf("Singlet MP2 correlation energy [au]:          %17.12lf\n", emp2tot_0);
    so <<indent
      <<scprintf("Triplet MP2 correlation energy [au]:          %17.12lf\n", emp2tot_1);
    so <<indent
      <<scprintf("Singlet (MP2)-R12/%2s correlation energy [au]: %17.12lf\n", SA_str, er12tot_0);
    so <<indent
      <<scprintf("Triplet (MP2)-R12/%2s correlation energy [au]: %17.12lf\n", SA_str, er12tot_1);
    so <<indent
      <<scprintf("Singlet MP2-R12/%2s correlation energy [au]:   %17.12lf\n", SA_str,
                 emp2tot_0 + er12tot_0);
    so <<indent
      <<scprintf("Triplet MP2-R12/%2s correlation energy [au]:   %17.12lf\n", SA_str,
                 emp2tot_1 + er12tot_1);
  }
  
  double etotal = escf + mp2_corr_energy_ + r12_corr_energy_;
  so <<endl<<indent
    <<scprintf("RHF energy [au]:                           %17.12lf\n", escf);
  so <<indent
    <<scprintf("MP2 correlation energy [au]:               %17.12lf\n", mp2_corr_energy_);
  so <<indent
    <<scprintf("(MBPT2)-R12/%2s correlation energy [au]:    %17.12lf\n", SA_str, r12_corr_energy_);
  so <<indent
    <<scprintf("MBPT2-R12/%2s correlation energy [au]:      %17.12lf\n", SA_str,
               mp2_corr_energy_ + r12_corr_energy_);
  so <<indent
    <<scprintf("MBPT2-R12/%2s energy [au]:                  %17.12lf\n", SA_str, etotal) << endl;

  so.flush();

  free(SA_str);
  
  return;
}


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

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