uscf.cc

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

    UBExtrapErrorOp(UnrestrictedSCF *s) : scf_(s) {}
    ~UBExtrapErrorOp() {}

    int has_side_effects() { return 1; }

    void process(SCMatrixBlockIter& bi) {
      int ir=current_block();
      
      for (bi.reset(); bi; bi++) {
        int i=bi.i();
        int j=bi.j();
        if (scf_->beta_occupation(ir,i) == scf_->beta_occupation(ir,j))
          bi.set(0.0);
      }
    }
};

Ref<SCExtrapData>
UnrestrictedSCF::extrap_data()
{
  Ref<SCExtrapData> data =
    new SymmSCMatrix2SCExtrapData(focka_.result_noupdate(),
                                  fockb_.result_noupdate());
  return data;
}

Ref<SCExtrapError>
UnrestrictedSCF::extrap_error()
{
  RefSCMatrix so_to_ortho_so_tr = so_to_orthog_so().t();

  // form Error_a
  RefSymmSCMatrix moa(oso_dimension(), basis_matrixkit());
  moa.assign(0.0);
  moa.accumulate_transform(so_to_ortho_so_tr * oso_scf_vector_,
                           focka_.result_noupdate(),
                           SCMatrix::TransposeTransform);
  
  Ref<SCElementOp> op = new UAExtrapErrorOp(this);
  moa.element_op(op.pointer());
  
  // form Error_b
  RefSymmSCMatrix mob(oso_dimension(), basis_matrixkit());
  mob.assign(0.0);
  mob.accumulate_transform(so_to_ortho_so_tr * oso_scf_vector_beta_,
                           fockb_.result_noupdate(),
                           SCMatrix::TransposeTransform);
  
  op = new UBExtrapErrorOp(this);
  mob.element_op(op);

  RefSymmSCMatrix aoa(so_dimension(), basis_matrixkit());
  aoa.assign(0.0);
  aoa.accumulate_transform(so_to_ortho_so_tr * oso_scf_vector_, moa);
  moa = 0;

  RefSymmSCMatrix aob(so_dimension(), basis_matrixkit());
  aob.assign(0.0);
  aob.accumulate_transform(so_to_ortho_so_tr * oso_scf_vector_beta_,mob);
  mob=0;

  aoa.accumulate(aob);
  aob=0;

  Ref<SCExtrapError> error = new SymmSCMatrixSCExtrapError(aoa);
  aoa=0;

  return error;
}

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

double
UnrestrictedSCF::compute_vector(double& eelec)
{
  tim_enter("vector");
  int i;

    // reinitialize the extrapolation object
  extrap_->reinitialize();
  
  // create level shifter
  ALevelShift *alevel_shift = new ALevelShift(this);
  alevel_shift->reference();
  BLevelShift *blevel_shift = new BLevelShift(this);
  blevel_shift->reference();
  
  // calculate the core Hamiltonian
  hcore_ = core_hamiltonian();

  // add density independant contributions to Hcore
  accumdih_->accum(hcore_);

  // set up subclass for vector calculation
  init_vector();
  
  // calculate the nuclear repulsion energy
  double nucrep = molecule()->nuclear_repulsion_energy();
  ExEnv::out0() << indent
       << scprintf("nuclear repulsion energy = %15.10f", nucrep)
       << endl << endl;

  RefDiagSCMatrix evalsa(oso_dimension(), basis_matrixkit());
  RefDiagSCMatrix evalsb(oso_dimension(), basis_matrixkit());

  double delta = 1.0;
  int iter;
  
  for (iter=0; iter < maxiter_; iter++) {
    // form the density from the current vector 
    tim_enter("density");
    delta = new_density();
    tim_exit("density");
    
    // check convergence
    if (delta < desired_value_accuracy())
      break;

    // reset the density from time to time
    if (iter && !(iter%dens_reset_freq_))
      reset_density();
      
    // form the AO basis fock matrix
    tim_enter("fock");
    double accuracy = 0.01 * delta;
    if (accuracy > 0.0001) accuracy = 0.0001;
    ao_fock(accuracy);
    tim_exit("fock");

    // calculate the electronic energy
    eelec = scf_energy();
    ExEnv::out0() << indent
         << scprintf("iter %5d energy = %15.10f delta = %10.5e",
                     iter+1, eelec+nucrep, delta)
         << endl;

    // now extrapolate the fock matrix
    tim_enter("extrap");
    Ref<SCExtrapData> data = extrap_data();
    Ref<SCExtrapError> error = extrap_error();
    extrap_->extrapolate(data,error);
    data=0;
    error=0;
    tim_exit("extrap");

    // diagonalize effective MO fock to get MO vector
    tim_enter("evals");

    RefSCMatrix so_to_oso_tr = so_to_orthog_so().t();

    RefSymmSCMatrix moa(oso_dimension(), basis_matrixkit());
    moa.assign(0.0);
    moa.accumulate_transform(so_to_oso_tr * oso_scf_vector_,
                             focka_.result_noupdate(),
                             SCMatrix::TransposeTransform);
    
    RefSymmSCMatrix mob(oso_dimension(), basis_matrixkit());
    mob.assign(0.0);
    mob.accumulate_transform(so_to_oso_tr * oso_scf_vector_beta_,
                             fockb_.result_noupdate(),
                             SCMatrix::TransposeTransform);
    
    RefSCMatrix nvectora(oso_dimension(), oso_dimension(), basis_matrixkit());
    RefSCMatrix nvectorb(oso_dimension(), oso_dimension(), basis_matrixkit());
  
    // level shift effective fock in the mo basis
    alevel_shift->set_shift(level_shift_);
    moa.element_op(alevel_shift);
    blevel_shift->set_shift(level_shift_);
    mob.element_op(blevel_shift);

    // transform back to the oso basis to do the diagonalization
    RefSymmSCMatrix osoa(oso_dimension(), basis_matrixkit());
    osoa.assign(0.0);
    osoa.accumulate_transform(oso_scf_vector_,moa);
    moa = 0;
    osoa.diagonalize(evalsa,oso_scf_vector_);
    osoa = 0;

    RefSymmSCMatrix osob(oso_dimension(), basis_matrixkit());
    osob.assign(0.0);
    osob.accumulate_transform(oso_scf_vector_beta_,mob);
    mob = 0;
    osob.diagonalize(evalsb,oso_scf_vector_beta_);
    osob = 0;

    tim_exit("evals");

    // now un-level shift eigenvalues
    alevel_shift->set_shift(-level_shift_);
    evalsa.element_op(alevel_shift);
    blevel_shift->set_shift(-level_shift_);
    evalsb.element_op(blevel_shift);
    
    if (reset_occ_)
      set_occupations(evalsa, evalsb);

    savestate_iter(iter);
    }
      
  eigenvalues_ = evalsa;
  eigenvalues_.computed() = 1;
  eigenvalues_.set_actual_accuracy(delta);
  evalsa = 0;
  
  oso_eigenvectors_ = oso_scf_vector_;
  oso_eigenvectors_.computed() = 1;
  oso_eigenvectors_.set_actual_accuracy(delta);
  
  oso_eigenvectors_beta_ = oso_scf_vector_beta_;
  oso_eigenvectors_beta_.computed() = 1;
  oso_eigenvectors_beta_.set_actual_accuracy(delta);

  eigenvalues_beta_ = evalsb;
  eigenvalues_beta_.computed() = 1;
  eigenvalues_beta_.set_actual_accuracy(delta);
  evalsb = 0;
  
  {
    // compute spin contamination
    RefSCMatrix so_to_oso_tr = so_to_orthog_so().t();
    RefSCMatrix Sab
      = (so_to_oso_tr * oso_scf_vector_).t()
      * overlap()
      * (so_to_oso_tr * oso_scf_vector_beta_);
    //Sab.print("Sab");
    BlockedSCMatrix *pSab = dynamic_cast<BlockedSCMatrix*>(Sab.pointer());
    double s2=0;
    for (int ir=0; ir < nirrep_; ir++) {
      RefSCMatrix Sab_ir=pSab->block(0);
      if (Sab_ir.nonnull()) {
        for (i=0; i < nalpha_[ir]; i++)
          for (int j=0; j < nbeta_[ir]; j++)
            s2 += Sab_ir.get_element(i,j)*Sab_ir.get_element(i,j);
      }
    }

    double S2real = (double)(tnalpha_-tnbeta_)/2.;
    S2real = S2real*(S2real+1);
    double S2 = S2real + tnbeta_ - s2;

    ExEnv::out0() << endl
         << indent << scprintf("<S^2>exact = %f", S2real) << endl
         << indent << scprintf("<S^2>      = %f", S2) << endl;
  }
  
  // now clean up
  done_vector();

  alevel_shift->dereference();
  delete alevel_shift;
  blevel_shift->dereference();
  delete blevel_shift;

  tim_exit("vector");
  //tim_print(0);

  return delta;
}

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

void
UnrestrictedSCF::init_gradient()
{
  // presumably the eigenvectors have already been computed by the time
  // we get here
  oso_scf_vector_ = oso_eigenvectors_.result_noupdate();
  oso_scf_vector_beta_ = oso_eigenvectors_beta_.result_noupdate();
}

void
UnrestrictedSCF::done_gradient()
{
  densa_=0;
  densb_=0;
  oso_scf_vector_ = 0;
  oso_scf_vector_beta_ = 0;
}

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

RefSymmSCMatrix
UnrestrictedSCF::lagrangian()
{
  RefSCMatrix so_to_oso_tr = so_to_orthog_so().t();

  RefDiagSCMatrix ea = eigenvalues_.result_noupdate().copy();
  RefDiagSCMatrix eb = eigenvalues_beta_.result_noupdate().copy();
  
  BlockedDiagSCMatrix *eab = dynamic_cast<BlockedDiagSCMatrix*>(ea.pointer());
  BlockedDiagSCMatrix *ebb = dynamic_cast<BlockedDiagSCMatrix*>(eb.pointer());

  Ref<PetiteList> pl = integral()->petite_list(basis());

  for (int ir=0; ir < nirrep_; ir++) {
    RefDiagSCMatrix eair = eab->block(ir);
    RefDiagSCMatrix ebir = ebb->block(ir);

    if (eair.null())
      continue;

    int i;
    for (i=nalpha_[ir]; i < eair.dim().n(); i++)
      eair.set_element(i,0.0);
    for (i=nbeta_[ir]; i < ebir.dim().n(); i++)
      ebir.set_element(i,0.0);
  }
  
  RefSymmSCMatrix la = basis_matrixkit()->symmmatrix(so_dimension());
  la.assign(0.0);
  la.accumulate_transform(so_to_oso_tr * oso_scf_vector_, ea);

  RefSymmSCMatrix lb = la.clone();
  lb.assign(0.0);
  lb.accumulate_transform(so_to_oso_tr * oso_scf_vector_beta_, eb);

  la.accumulate(lb);
  
  la = pl->to_AO_basis(la);
  la->scale(-1.0);

  return la;
}

RefSymmSCMatrix
UnrestrictedSCF::gradient_density()
{
  densa_ = basis_matrixkit()->symmmatrix(so_dimension());
  densb_ = densa_.clone();
  
  so_density(densa_, 1.0, 1);
  so_density(densb_, 1.0, 0);
  
  Ref<PetiteList> pl = integral()->petite_list(basis());
  
  densa_ = pl->to_AO_basis(densa_);
  densb_ = pl->to_AO_basis(densb_);

  RefSymmSCMatrix tdens = densa_.copy();
  tdens.accumulate(densb_);
  return tdens;
}

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

void
UnrestrictedSCF::init_hessian()
{
}

void
UnrestrictedSCF::done_hessian()
{
}

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

void
UnrestrictedSCF::two_body_deriv_hf(double * tbgrad, double exchange_fraction)
{
  Ref<SCElementMaxAbs> m = new SCElementMaxAbs;
  densa_.element_op(m.pointer());
  double pmax = m->result();
  m=0;

  // now try to figure out the matrix specialization we're dealing with.
  // if we're using Local matrices, then there's just one subblock, or
  // see if we can convert P to local matrices

  if (local_ || local_dens_) {
    // grab the data pointers from the P matrices
    double *pmata, *pmatb;
    RefSymmSCMatrix ptmpa = get_local_data(densa_, pmata, SCF::Read);
    RefSymmSCMatrix ptmpb = get_local_data(densb_, pmatb, SCF::Read);
  
    Ref<PetiteList> pl = integral()->petite_list();
    LocalUHFGradContribution l(pmata,pmatb);

    int i;
    int na3 = molecule()->natom()*3;
    int nthread = threadgrp_->nthread();
    double **grads = new double*[nthread];
    Ref<TwoBodyDerivInt> *tbis = new Ref<TwoBodyDerivInt>[nthread];
    for (i=0; i < nthread; i++) { 
      tbis[i] = integral()->electron_repulsion_deriv();
      grads[i] = new double[na3];
      memset(grads[i], 0, sizeof(double)*na3);
    }

    LocalTBGrad<LocalUHFGradContribution> **tblds =
      new LocalTBGrad<LocalUHFGradContribution>*[nthread];

    for (i=0; i < nthread; i++) {
      tblds[i] = new LocalTBGrad<LocalUHFGradContribution>(
        l, tbis[i], pl, basis(), scf_grp_, grads[i], pmax,
        desired_gradient_accuracy(), nthread, i, exchange_fraction);
      threadgrp_->add_thread(i, tblds[i]);
    }

    if (threadgrp_->start_threads() < 0
        ||threadgrp_->wait_threads() < 0) {
      ExEnv::err0() << indent
           << "USCF: error running threads" << endl;
      abort();
    }

    for (i=0; i < nthread; i++) {
      for (int j=0; j < na3; j++)
        tbgrad[j] += grads[i][j];

      delete[] grads[i];
      delete tblds[i];
    }

    scf_grp_->sum(tbgrad,3 * basis()->molecule()->natom());
  }

  // for now quit
  else {
    ExEnv::err0() << indent
         << "USCF::two_body_deriv_hf: can't do gradient yet\n";
    abort();
  }
}

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

}

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