compute_sbs_a.cc

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

	      }
	    }
	  }
	  ij_index++;
	}
      }
    }
    }

#if PRINT4Q
    if (me == 0) {
      int index = 0;
      int ij_index = 0;
      int j_offset = nfzc;
      for (int i = 0; i<ni; i++) {
	for (int j = 0; j<nocc_act; j++) {
	  if (index++ % nproc == me) {
	    double *integral_ij_offset = mo_int + num_te_types*nbasis*nbasis*ij_index;
	    for(int te_type=0; te_type<PRINT_NUM_TE_TYPES; te_type++,integral_ij_offset+=nbasis*nbasis) {
	      for (int y = 0; y < noso; y++) {
		double *integral_ijyx_ptr = integral_ij_offset + y*nbasis;
		for (int x = 0; x<noso; x++) {
		  printf("4Q: type = %d (%d %d|%d %d) = %12.8f\n",
		       te_type,i+i_offset,x,j+j_offset,y,*integral_ijyx_ptr);
		integral_ijyx_ptr++;
		}
	      }
	    }
	    ij_index++;
	  }
	}
      }
    }
#endif

    // For now compute MP2 energy to verify the transformed ERIs
    tim_enter("compute emp2");

    index = 0;
    ij_index = 0;
    for (int i=0; i<ni; i++) {
      int ii = i + i_offset - nfzc;
      for (int j=0; j<nocc_act; j++) {
	int jj = j;
        double ecorr_ij = 0.0;
	int ij_aa = (ii > jj) ? ii*(ii-1)/2 + jj : jj*(jj-1)/2 + ii;
	int ij_ab = ii*nocc_act + jj;
	double eaa = 0.0;
	double eab = 0.0;

        if (index++ % nproc == me) {
	  for (int b=0; b<nvir; b++) {
	    double* iajb_ptr = &mo_int[nocc + nbasis*(b+nocc + nbasis*num_te_types*ij_index)];
	    double* ibja_ptr = &mo_int[b+nocc + nbasis*(nocc + nbasis*num_te_types*ij_index)];
	    for (int a=0; a<nvir; a++) {
#if PRINT4Q_MP2
	      printf("4Q: (%d %d|%d %d) = %12.8f\n",
		     i,a+nocc,j+nfzc,b+nocc,*iajb_ptr);
#endif
	      double delta_ijab = evals[i_offset+i]+evals[j+nfzc]-evals[nocc+a]-evals[nocc+b];
	      // only include determinants with unique coefficients
	      if (a>=b && i_offset+i>=j+nfzc) {
		if (a>b && i_offset+i>j+nfzc) {
		  // aaaa or bbbb
		  biggest_coefs.insert(*iajb_ptr - *ibja_ptr,
				       i_offset+i,j,a,b,1111);
		  // aabb or bbaa or abba or baab
		  biggest_coefs.insert(*ibja_ptr,i_offset+i,j,b,a,1212);
		} // endif
		// aabb or bbaa or abba or baab
		biggest_coefs.insert(*iajb_ptr,i_offset+i,j,a,b,1212);
	      } // endif
	      double tmpval = (*iajb_ptr - *ibja_ptr)*(*iajb_ptr - *ibja_ptr)/delta_ijab;
	      eaa += tmpval;
	      tmpval = 0.5*(*iajb_ptr * *iajb_ptr + *ibja_ptr * *ibja_ptr)/delta_ijab;
	      eab += tmpval;
	      ecorr_ij += *iajb_ptr*(2.0 * *iajb_ptr - *ibja_ptr)/delta_ijab;
	      iajb_ptr++;
	      ibja_ptr += nbasis;
	    } // exit a loop
	  }   // exit b loop
	  ij_index++;
	  if (ii > jj)
	    emp2_aa.set_element(ij_aa,eaa);
	  emp2_ab.set_element(ij_ab,eab);
	}     // endif
	if (debug_) {
	  msg->sum(ecorr_ij);
	  ExEnv::out0() << indent
			<< scprintf("correlation energy for pair %3d %3d = %16.12f",
				    i+i_offset, j+nfzc, ecorr_ij)
			<< endl;
	}
      }         // exit j loop
    }           // exit i loop
    tim_exit("compute emp2");
    
    // debug print
    if (debug_) {
      ExEnv::out0() << indent << "End of ecorr" << endl;
    }
    // end of debug print
    
    integral_ijsq = 0;
    mem->sync(); // Make sure MO integrals are complete on all nodes before continuing

    // Push locally stored integrals to an accumulator
    // This could involve storing the data to disk or simply remembering the pointer
    tim_enter("MO ints store");
    r12intsacc->store_memorygrp(mem,ni);
    tim_exit("MO ints store");
    mem->sync();

    if (me == 0 && mole->if_to_checkpoint() && r12intsacc->can_restart()) {
      current_orbital_ += ni;
      StateOutBin stateout(mole->checkpoint_file());
      SavableState::save_state(mole,stateout);
      ExEnv::out0() << indent << "Checkpointed the wave function" << endl;
    }

  } // end of loop over passes
  tim_exit("mp2-r12/a passes");
  if (debug_)
    ExEnv::out0() << indent << "End of mp2-r12/a transformation" << endl;
  // Done storing integrals - commit the content
  // WARNING: it is not safe to use mem until deactivate has been called on the accumulator
  //          After that deactivate the size of mem will be 0 [mem->set_localsize(0)]
  r12intsacc->commit();


  /*-----------------------------------------------
    Compute dipole and quadrupole moment integrals
   -----------------------------------------------*/
  RefSymmSCMatrix MX, MY, MZ, MXX, MYY, MZZ;
  r12info()->compute_multipole_ints(MX,MY,MZ,MXX,MYY,MZZ);
  if (debug_)
    ExEnv::out0() << indent << "Computed multipole moment integrals" << endl;
  
  /*--------------------------------
    Compute MP2-R12/A intermediates
    and collect on node0
   --------------------------------*/
  tim_enter("mp2-r12a intermeds");
  int naa = (nocc_act*(nocc_act-1))/2;          // Number of alpha-alpha pairs
  int nab = nocc_act*nocc_act;                  // Number of alpha-beta pairs
  if (debug_) {
    ExEnv::out0() << indent << "naa = " << naa << endl;
    ExEnv::out0() << indent << "nab = " << nab << endl;
  }
  double *Vaa_ijkl = new double[naa*naa];
  double *Taa_ijkl = new double[naa*naa];
  double *Xaa_ijkl = new double[naa*naa];
  double *Vab_ijkl = new double[nab*nab];
  double *Tab_ijkl = new double[nab*nab];
  double *Xab_ijkl = new double[nab*nab];
  if (debug_)
    ExEnv::out0() << indent << "Allocated intermediates V, X, and T" << endl;
  bzerofast(Vaa_ijkl,naa*naa);
  bzerofast(Taa_ijkl,naa*naa);
  bzerofast(Xaa_ijkl,naa*naa);
  bzerofast(Vab_ijkl,nab*nab);
  bzerofast(Tab_ijkl,nab*nab);
  bzerofast(Xab_ijkl,nab*nab);

  // Compute intermediates
  if (debug_)
    ExEnv::out0() << indent << "Ready to compute intermediates V, X, and T" << endl;
  const int pair_block_size = num_te_types*nbasis*nbasis;
  const double oosqrt2 = 1.0/sqrt(2.0);
  // Compute the number of tasks that have full access to the integrals
  // and split the work among them
  int nproc_with_ints = 0;
  for(int proc=0;proc<nproc;proc++)
    if (r12intsacc->has_access(proc)) nproc_with_ints++;
  int *proc_with_ints = new int[nproc];
  int count = 0;
  for(int proc=0;proc<nproc;proc++)
    if (r12intsacc->has_access(proc)) {
      proc_with_ints[proc] = count;
      count++;
    }
    else
      proc_with_ints[proc] = -1;
  if (debug_)
    ExEnv::out0() << indent << "Computing intermediates on " << nproc_with_ints
		  << " processors" << endl;

  
  //////////////////////////////////////////////////////////////
  //
  // Evaluation of the intermediates proceeds as follows:
  //
  //    loop over batches of kl, k >= l,  0<=k,l<nocc_act
  //      load (kl|xy), (kl| [T1,r12] |xy), and (lk| [T1,r12] |xy)
  //           (aka kl-sets) into memory
  //
  //      loop over batches of ij, i>=j, 0<=i,j<nocc_act
  //        load (ij|r12|xy) into memory
  //           (aka ij-sets) into memory
  //        compute V[ij][kl] and T[ij][kl] for all ij and kl in
  //                the "direct product" batch
  //      end ij loop
  //    end kl loop
  //
  /////////////////////////////////////////////////////////////////////////////////

  if (r12intsacc->has_access(me)) {
    int kl = 0;
    for(int k=0;k<nocc_act;k++)
      for(int l=0;l<=k;l++,kl++) {
        int kl_proc = kl%nproc_with_ints;
        if (kl_proc != proc_with_ints[me])
          continue;
	int kl_aa = k*(k-1)/2 + l;
	int kl_ab = k*nocc_act + l;
	int lk_ab = l*nocc_act + k;
        
        // Get (|r12|) and (|[r12,T1]|) integrals only
        tim_enter("MO ints retrieve");
        double *klyx_buf_eri = r12intsacc->retrieve_pair_block(k,l,R12IntsAcc::eri);
        double *klyx_buf_r12 = r12intsacc->retrieve_pair_block(k,l,R12IntsAcc::r12);
        double *klyx_buf_r12t1 = r12intsacc->retrieve_pair_block(k,l,R12IntsAcc::r12t1);
	double *lkyx_buf_r12t1 = r12intsacc->retrieve_pair_block(l,k,R12IntsAcc::r12t1);
        tim_exit("MO ints retrieve");

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


	    int ij_aa = i*(i-1)/2 + j;
	    int ij_ab = i*nocc_act + j;
	    int ji_ab = j*nocc_act + i;
            
            tim_enter("MO ints retrieve");
            double *ijyx_buf_r12 = r12intsacc->retrieve_pair_block(i,j,R12IntsAcc::r12);
            tim_exit("MO ints retrieve");

	    double *Vaa_ij = Vaa_ijkl + ij_aa*naa;
            double *Vab_ij = Vab_ijkl + ij_ab*nab;
            double *Vab_ji = Vab_ijkl + ji_ab*nab;
            double *Taa_ij = Taa_ijkl + ij_aa*naa;
            double *Tab_ij = Tab_ijkl + ij_ab*nab;
            double *Tab_ji = Tab_ijkl + ji_ab*nab;
            double *Xaa_ij = Xaa_ijkl + ij_aa*naa;
            double *Xab_ij = Xab_ijkl + ij_ab*nab;
            double *Xab_ji = Xab_ijkl + ji_ab*nab;
            

            tim_enter("MO ints contraction");
            double Vaa_ijkl, Vab_ijkl, Vab_jikl, Vab_ijlk, Vab_jilk;
            double Xaa_ijkl, Xab_ijkl, Xab_jikl, Xab_ijlk, Xab_jilk;
	    double Taa_ijkl, Tab_ijkl, Tab_jikl, Tab_ijlk, Tab_jilk;
	    Vaa_ijkl = Vab_ijkl = Vab_jikl = Vab_ijlk = Vab_jilk = 0.0;
	    Xaa_ijkl = Xab_ijkl = Xab_jikl = Xab_ijlk = Xab_jilk = 0.0;
	    Taa_ijkl = Tab_ijkl = Tab_jikl = Tab_ijlk = Tab_jilk = 0.0;

	    /*----------------------------------
	      Compute (r12)^2 contribution to X
	     ----------------------------------*/
	    double r1r1_ik = -1.0*(MXX->get_element(i,k) + MYY->get_element(i,k) + MZZ->get_element(i,k));
	    double r1r1_il = -1.0*(MXX->get_element(i,l) + MYY->get_element(i,l) + MZZ->get_element(i,l));
	    double r1r1_jk = -1.0*(MXX->get_element(j,k) + MYY->get_element(j,k) + MZZ->get_element(j,k));
	    double r1r1_jl = -1.0*(MXX->get_element(j,l) + MYY->get_element(j,l) + MZZ->get_element(j,l));
	    double r1r2_ijkl = MX->get_element(i,k)*MX->get_element(j,l) +
	      MY->get_element(i,k)*MY->get_element(j,l) +
	      MZ->get_element(i,k)*MZ->get_element(j,l);
	    double r1r2_ijlk = MX->get_element(i,l)*MX->get_element(j,k) +
	      MY->get_element(i,l)*MY->get_element(j,k) +
	      MZ->get_element(i,l)*MZ->get_element(j,k);
	    double delta_ik = (i==k ? 1.0 : 0.0);
	    double delta_il = (i==l ? 1.0 : 0.0);
	    double delta_jk = (j==k ? 1.0 : 0.0);
	    double delta_jl = (j==l ? 1.0 : 0.0);
	    Xab_ijkl += r1r1_ik * delta_jl + r1r1_jl * delta_ik - 2.0*r1r2_ijkl;
	    if (i != j)
	      Xab_jikl += r1r1_jk * delta_il + r1r1_il * delta_jk - 2.0*r1r2_ijlk;
	    if (k != l)
	      Xab_ijlk += r1r1_il * delta_jk + r1r1_jk * delta_il - 2.0*r1r2_ijlk;
	    if (i != j && k != l) {
	      Xaa_ijkl += r1r1_ik * delta_jl + r1r1_jl * delta_ik - 2.0*r1r2_ijkl -
		r1r1_jk * delta_il - r1r1_il * delta_jk + 2.0*r1r2_ijlk;
	      Xab_jilk += r1r1_ik * delta_jl + r1r1_jl * delta_ik - 2.0*r1r2_ijkl;
	    }

            for(int y=0;y<noso;y++) {
	      double pfac_xy;
	      if (y >= nocc)
		pfac_xy = pfac_xy_1;
	      else
		pfac_xy = pfac_xy_2;
              for(int x=0;x<noso;x++) {
                int yx_offset = y*nbasis+x;
                int xy_offset = x*nbasis+y;
                double ij_r12_xy = ijyx_buf_r12[yx_offset];
                double ij_r12_yx = ijyx_buf_r12[xy_offset];
                double kl_eri_xy = klyx_buf_eri[yx_offset];
                double kl_eri_yx = klyx_buf_eri[xy_offset];
                Vab_ijkl -= pfac_xy * (ij_r12_xy * kl_eri_xy + ij_r12_yx * kl_eri_yx);
		if (i != j)
		  Vab_jikl -= pfac_xy * (ij_r12_yx * kl_eri_xy + ij_r12_xy * kl_eri_yx);
		if (k != l)
		  Vab_ijlk -= pfac_xy * (ij_r12_xy * kl_eri_yx + ij_r12_yx * kl_eri_xy);
		if (i != j && k != l) {
		  Vaa_ijkl -= pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_eri_xy - kl_eri_yx);
		  Vab_jilk -= pfac_xy * (ij_r12_yx * kl_eri_yx + ij_r12_xy * kl_eri_xy);
		}
                double kl_r12_xy = klyx_buf_r12[yx_offset];
                double kl_r12_yx = klyx_buf_r12[xy_offset];
                Xab_ijkl -= pfac_xy * (ij_r12_xy * kl_r12_xy + ij_r12_yx * kl_r12_yx);
		if (i != j)
		  Xab_jikl -= pfac_xy * (ij_r12_yx * kl_r12_xy + ij_r12_xy * kl_r12_yx);
		if (k != l)
		  Xab_ijlk -= pfac_xy * (ij_r12_xy * kl_r12_yx + ij_r12_yx * kl_r12_xy);
		if (i != j && k != l) {
		  Xaa_ijkl -= pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_r12_xy - kl_r12_yx);
		  Xab_jilk -= pfac_xy * (ij_r12_yx * kl_r12_yx + ij_r12_xy * kl_r12_xy);
		}
                double kl_r12t1_xy = klyx_buf_r12t1[yx_offset];
                double kl_r12t1_yx = klyx_buf_r12t1[xy_offset];
                double lk_r12t1_xy = lkyx_buf_r12t1[yx_offset];
                double lk_r12t1_yx = lkyx_buf_r12t1[xy_offset];
                double kl_Tr12_xy = -kl_r12t1_xy-lk_r12t1_yx;
                double kl_Tr12_yx = -kl_r12t1_yx-lk_r12t1_xy;
                Tab_ijkl += pfac_xy * (ij_r12_xy * kl_Tr12_xy + ij_r12_yx * kl_Tr12_yx);
		if (i != j)
		  Tab_jikl += pfac_xy * (ij_r12_yx * kl_Tr12_xy + ij_r12_xy * kl_Tr12_yx);
		if (k != l)
		  Tab_ijlk += pfac_xy * (ij_r12_xy * kl_Tr12_yx + ij_r12_yx * kl_Tr12_xy);
		if (i != j && k != l) {
		  Taa_ijkl += pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_Tr12_xy - kl_Tr12_yx);
		  Tab_jilk += pfac_xy * (ij_r12_yx * kl_Tr12_yx + ij_r12_xy * kl_Tr12_xy);
		}
              }
	    }
            Vab_ij[kl_ab] += Vab_ijkl;
	    if (i != j)
	      Vab_ji[kl_ab] += Vab_jikl;
	    if (k != l)
	      Vab_ij[lk_ab] += Vab_ijlk;
	    if (i != j && k != l) {
	      Vaa_ij[kl_aa] += Vaa_ijkl;
	      Vab_ji[lk_ab] += Vab_jilk;
	    }
            Xab_ij[kl_ab] += Xab_ijkl;
	    if (i != j)
	      Xab_ji[kl_ab] += Xab_jikl;
	    if (k != l)
	      Xab_ij[lk_ab] += Xab_ijlk;