solvent.cc

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

      efdn_mat->assign(0.0);
      efdn_mat->element_op(efdn_op);
      sp->init();
      efdn_mat->element_op(generic_sp, ao_density);
      efield_dot_normals_[i] = sp->result();
    }
  RefSCDimension aodim = ao_density.dim();
  Ref<SCMatrixKit> aokit = ao_density.kit();
  ao_density = 0;
  efdn_mat = 0;
  tim_exit("efield");

  // compute a new set of charges
  tim_enter("charges");
  // electron contrib
  solvent_->compute_charges(efield_dot_normals_, charges_);
  double qeenc = solvent_->computed_enclosed_charge();
  // nuclear contrib
  for (i=0; i<ncharge; i++) {
      double nuc_efield[3];
      wfn_->molecule()->nuclear_efield(charge_positions_[i], nuc_efield);
      double tmp = 0.0;
      for (j=0; j<3; j++) {
          tmp += nuc_efield[j] * normals_[i][j];
        }
      efield_dot_normals_[i] = tmp;
    }
  solvent_->compute_charges(efield_dot_normals_, charges_n_);
  double qnenc = solvent_->computed_enclosed_charge();
  tim_exit("charges");

  // normalize the charges
  // e and n are independently normalized since the nature of the
  // errors in e and n are different: n error is just numerical and
  // e error is numerical plus diffuseness of electron distribution
  if (normalize_q_) {
      tim_enter("norm");
      // electron contrib
      solvent_->normalize_charge(-wfn_->nelectron(), charges_);
      // nuclear contrib
      solvent_->normalize_charge(wfn_->molecule()->nuclear_charge(),
                                 charges_n_);
      tim_exit("norm");
    }
  // sum the nuclear and electron contrib
  for (i=0; i<ncharge; i++) charges_[i] += charges_n_[i];

  //// compute scalar contributions
  double A = solvent_->area();

  // the cavitation energy
  ecavitation_ = A * gamma_;

  // compute the nuclear-surface interaction energy
  tim_enter("n-s");
  enucsurf_
      = solvent_->nuclear_interaction_energy(charge_positions_, charges_);
  tim_exit("n-s");

  double enqn = 0.0, enqe = 0.0;
  if (y_equals_j_ || separate_surf_charges_) {
      tim_enter("n-qn");
      enqn = solvent_->nuclear_interaction_energy(charge_positions_,
                                                  charges_n_);
      enqe = enucsurf_ - enqn;
      tim_exit("n-qn");
    }

  //// compute one body contributions

  // compute the electron-surface interaction matrix elements
  tim_enter("e-s");
  Ref<PointChargeData> pc_dat = new PointChargeData(ncharge,
                                                  charge_positions_, charges_);
  Ref<OneBodyInt> pc = wfn_->integral()->point_charge(pc_dat);
  Ref<SCElementOp> pc_op = new OneBodyIntOp(pc);

  // compute matrix elements in the ao basis
  RefSymmSCMatrix h_ao(aodim, aokit);
  h_ao.assign(0.0);
  h_ao.element_op(pc_op);
  // transform to the so basis and add to h
  RefSymmSCMatrix h_so = wfn_->integral()->petite_list()->to_SO_basis(h_ao);
  if (onebody_) h->accumulate(h_so);
  // compute the contribution to the energy
  sp->init();
  RefSymmSCMatrix so_density = wfn_->density()->copy();
  // for the scalar products, scale the density's off-diagonals by two
  so_density->scale(2.0);
  so_density->scale_diagonal(0.5);
  h_so->element_op(generic_sp, so_density);
  eelecsurf_ = sp->result();
  tim_exit("e-s");

  double eeqn = 0.0, eeqe = 0.0;
  if (y_equals_j_ || separate_surf_charges_) {
      tim_enter("e-qn");
      pc_dat = new PointChargeData(ncharge, charge_positions_, charges_n_);
      pc = wfn_->integral()->point_charge(pc_dat);
      pc_op = new OneBodyIntOp(pc);

      // compute matrix elements in the ao basis
      h_ao.assign(0.0);
      h_ao.element_op(pc_op);
      // transform to the so basis
      h_so = wfn_->integral()->petite_list()->to_SO_basis(h_ao);
      // compute the contribution to the energy
      sp->init();
      h_so->element_op(generic_sp, so_density);
      eeqn = sp->result();
      eeqe = eelecsurf_ - eeqn;
      tim_exit("e-qn");
    }

  if (y_equals_j_) {
      // Remove the y term (enqe) and add the j term (eeqn).  Formally,
      // they are equal, but they are not because some e-density is outside
      // the surface and because of the numerical approximations.
      enucsurf_ += eeqn - enqe;
    }

  // compute the surface-surface interaction energy
  esurfsurf_ = -0.5*(eelecsurf_+enucsurf_);
  // (this can also be computed as below, but is much more expensive)
  //tim_enter("s-s");
  //double esurfsurf_;
  //esurfsurf_ = solvent_->self_interaction_energy(charge_positions_, charges_);
  //tim_exit("s-s");

  escalar_ = enucsurf_ + esurfsurf_ + ecavitation_ + edisprep_;
  // NOTE: SCF currently only adds h_so to the Fock matrix
  // so a term is missing in the energy.  This term is added here
  // and when SCF is fixed, should no longer be included.
  if (onebody_) escalar_ += 0.5 * eelecsurf_;

  if (!onebody_) escalar_ += eelecsurf_;

  ExEnv::out0() << incindent;
  ExEnv::out0() << indent
       << "Solvent: "
       << scprintf("q(e-enc)=%12.10f q(n-enc)=%12.10f", qeenc, qnenc)
       << endl;
  ExEnv::out0() << incindent;
  if (separate_surf_charges_) {
      ExEnv::out0() << indent
           << scprintf("E(n-qn)=%10.8f ", enqn)
           << scprintf("E(n-qe)=%10.8f", enqe)
           << endl;
      ExEnv::out0() << indent
           << scprintf("E(e-qn)=%10.8f ", eeqn)
           << scprintf("E(e-qe)=%10.8f", eeqe)
           << endl;
      //ExEnv::out0() << indent
      //     << scprintf("DG = %12.8f ", 0.5*627.51*(enqn+enqe+eeqn+eeqe))
      //     << scprintf("DG(Y=J) = %12.8f", 0.5*627.51*(enqn+2*eeqn+eeqe))
      //     << endl;
    }
  ExEnv::out0() << indent
       << scprintf("E(c)=%10.8f ", ecavitation_)
       << scprintf("E(disp-rep)=%10.8f", edisprep_)
       << endl;
  ExEnv::out0() << indent
       << scprintf("E(n-s)=%10.8f ", enucsurf_)
       << scprintf("E(e-s)=%10.8f ", eelecsurf_)
       << scprintf("E(s-s)=%10.8f ", esurfsurf_)
       << endl;
  ExEnv::out0() << decindent;
  ExEnv::out0() << decindent;

  tim_exit("accum");
  tim_exit("solvent");
}

void
BEMSolventH::done()
{
  solvent_->free_normals(normals_);
  normals_ = 0;
  solvent_->free_efield_dot_normals(efield_dot_normals_);
  efield_dot_normals_ = 0;
  solvent_->free_charges(charges_);
  solvent_->free_charges(charges_n_);
  charges_ = 0;
  charges_n_ = 0;
  solvent_->free_charge_positions(charge_positions_);
  charge_positions_ = 0;
  solvent_->done();
}

void
BEMSolventH::print_summary()
{
  Ref<Units> unit = new Units("kcal/mol");
  ExEnv::out0() << endl;
  ExEnv::out0() << "Summary of solvation calculation:" << endl;
  ExEnv::out0() << "_______________________________________________" << endl;
  ExEnv::out0() << endl;
  ExEnv::out0().setf(ios::scientific,ios::floatfield); // use scientific format
  ExEnv::out0().precision(5);
  ExEnv::out0() << indent << "E(nuc-surf):              " 
       << setw(12) << setfill(' ')
       << enucsurf_*unit->from_atomic_units() << " kcal/mol" << endl; 
  ExEnv::out0() << indent << "E(elec-surf):             " 
       << setw(12) << setfill(' ')
       << eelecsurf_*unit->from_atomic_units() << " kcal/mol" << endl; 
  ExEnv::out0() << indent << "E(surf-surf):             " 
       << setw(12) << setfill(' ')
       << esurfsurf_*unit->from_atomic_units() << " kcal/mol" << endl; 
  ExEnv::out0() << indent << "Electrostatic energy:     " 
       << setw(12) << setfill(' ')
       << (enucsurf_+eelecsurf_+esurfsurf_)*unit->from_atomic_units()
       << " kcal/mol" << endl; 
  ExEnv::out0() << "_______________________________________________" << endl;
  ExEnv::out0() << endl;
  ExEnv::out0() << indent << "E(cav):                   " 
       << setw(12) << setfill(' ')
       << ecavitation_*unit->from_atomic_units() << " kcal/mol" << endl; 
  ExEnv::out0() << indent << "E(disp):                  " 
       << setw(12) << setfill(' ')
       << solvent_->disp()*unit->from_atomic_units() << " kcal/mol" << endl; 
  ExEnv::out0() << indent << "E(rep):                   " 
       << setw(12) << setfill(' ')
       << solvent_->rep()*unit->from_atomic_units() << " kcal/mol" << endl; 
  ExEnv::out0() << indent << "Non-electrostatic energy: "
       << setw(12) << setfill(' ')
       << (ecavitation_+solvent_->disp()+solvent_->rep())
          *unit->from_atomic_units() << " kcal/mol" << endl; 
  ExEnv::out0() << "_______________________________________________" << endl;

}

double
BEMSolventH::e()
{
  return escalar_;
}

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

}

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