// // mbpt.cc // // Copyright (C) 1996 Limit Point Systems, Inc. // // Author: Ida Nielsen // Maintainer: LPS // // This file is part of the SC Toolkit. // // The SC Toolkit is free software; you can redistribute it and/or modify // it under the terms of the GNU Library General Public License as published by // the Free Software Foundation; either version 2, or (at your option) // any later version. // // The SC Toolkit is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Library General Public License for more details. // // You should have received a copy of the GNU Library General Public License // along with the SC Toolkit; see the file COPYING.LIB. If not, write to // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. // // The U.S. Government is granted a limited license as per AL 91-7. // #ifdef __GNUC__ #pragma implementation #endif #include #include #include #include #include #include #include using namespace std; using namespace sc; ///////////////////////////////////////////////////////////////// // Function dquicksort performs a quick sort (smaller -> larger) // of the double data in item by the integer indices in index; // data in item remain unchanged ///////////////////////////////////////////////////////////////// static void dqs(double *item,int *index,int left,int right) { register int i,j; double x; int y; i=left; j=right; x=item[index[(left+right)/2]]; do { while(item[index[i]]left) j--; if (i<=j) { if (item[index[i]] != item[index[j]]) { y=index[i]; index[i]=index[j]; index[j]=y; } i++; j--; } } while(i<=j); if (left, create); MBPT2::MBPT2(StateIn& s): SavableState(s), Wavefunction(s) { reference_ << SavableState::restore_state(s); s.get(nfzc); s.get(nfzv); if (s.version(::class_desc()) >= 8) { double dmem_alloc; s.get(dmem_alloc); mem_alloc = size_t(dmem_alloc); } else { unsigned int imem_alloc; s.get(imem_alloc); mem_alloc = imem_alloc; } s.getstring(method_); s.getstring(algorithm_); if (s.version(::class_desc()) <= 6) { int debug_old; s.get(debug_old); } if (s.version(::class_desc()) >= 2) { s.get(do_d1_); } else { do_d1_ = 0; } if (s.version(::class_desc()) >= 3) { s.get(dynamic_); } else { dynamic_ = 0; } if (s.version(::class_desc()) >= 9) { s.get(print_percent_); } else { print_percent_ = 10.0; } if (s.version(::class_desc()) >= 4) { s.get(cphf_epsilon_); } else { cphf_epsilon_ = 1.0e-8; } if (s.version(::class_desc()) >= 5) { s.get(max_norb_); } else { max_norb_ = 0; } if (s.version(::class_desc()) >= 6) { s.get(do_d2_); } else { do_d2_ = 1; } hf_energy_ = 0.0; symorb_irrep_ = 0; symorb_num_ = 0; eliminate_in_gmat_ = 1; mem = MemoryGrp::get_default_memorygrp(); msg_ = MessageGrp::get_default_messagegrp(); thr_ = ThreadGrp::get_default_threadgrp(); restart_ecorr_ = 0.0; restart_orbital_v1_ = 0; restart_orbital_memgrp_ = 0; } MBPT2::MBPT2(const Ref& keyval): Wavefunction(keyval) { reference_ << keyval->describedclassvalue("reference"); if (reference_.null()) { ExEnv::err0() << "MBPT2::MBPT2: no reference wavefunction" << endl; abort(); } copy_orthog_info(reference_); nfzc = keyval->intvalue("nfzc"); char *nfzc_charval = keyval->pcharvalue("nfzc"); if (nfzc_charval && !strcmp(nfzc_charval, "auto")) { if (molecule()->max_z() > 30) { ExEnv::err0() << "MBPT2: cannot use \"nfzc = auto\" for Z > 30" << endl; abort(); } nfzc = molecule()->n_core_electrons()/2; ExEnv::out0() << indent << "MBPT2: auto-freezing " << nfzc << " core orbitals" << endl; } delete[] nfzc_charval; nfzv = keyval->intvalue("nfzv"); mem_alloc = keyval->sizevalue("memory"); if (keyval->error() != KeyVal::OK) { // by default, take half of the memory mem_alloc = ExEnv::memory()/2; if (mem_alloc == 0) { mem_alloc = DEFAULT_SC_MEMORY; } } mem << keyval->describedclassvalue("memorygrp"); if (mem.null()) { mem = MemoryGrp::get_default_memorygrp(); } msg_ = MessageGrp::get_default_messagegrp(); thr_ = ThreadGrp::get_default_threadgrp(); method_ = keyval->pcharvalue("method"); algorithm_ = keyval->pcharvalue("algorithm"); do_d1_ = keyval->booleanvalue("compute_d1"); do_d2_ = keyval->booleanvalue("compute_d2",KeyValValueboolean(1)); restart_ecorr_ = keyval->doublevalue("restart_ecorr"); restart_orbital_v1_ = keyval->intvalue("restart_orbital_v1"); restart_orbital_memgrp_ = keyval->intvalue("restart_orbital_memgrp"); KeyValValueint default_dynamic(0); dynamic_ = keyval->booleanvalue("dynamic", default_dynamic); KeyValValuedouble default_print_percent(10.0); print_percent_ = keyval->doublevalue("print_percent", default_print_percent); cphf_epsilon_ = keyval->doublevalue("cphf_epsilon",KeyValValuedouble(1.e-8)); max_norb_ = keyval->intvalue("max_norb",KeyValValueint(-1)); hf_energy_ = 0.0; symorb_irrep_ = 0; symorb_num_ = 0; eliminate_in_gmat_ = 1; } MBPT2::~MBPT2() { delete[] method_; delete[] algorithm_; delete[] symorb_irrep_; delete[] symorb_num_; } void MBPT2::save_data_state(StateOut& s) { Wavefunction::save_data_state(s); SavableState::save_state(reference_.pointer(),s); s.put(nfzc); s.put(nfzv); double dmem_alloc = mem_alloc; s.put(dmem_alloc); s.putstring(method_); s.putstring(algorithm_); s.put(do_d1_); s.put(dynamic_); s.put(print_percent_); s.put(cphf_epsilon_); s.put(max_norb_); s.put(do_d2_); } void MBPT2::print(ostream&o) const { o << indent << "MBPT2:" << endl; o << incindent; Wavefunction::print(o); o << indent << "Reference Wavefunction:" << endl; o << incindent; reference_->print(o); o << decindent << endl; o << decindent; } ////////////////////////////////////////////////////////////////////////////// int MBPT2::spin_polarized() { return reference_->spin_polarized(); } RefSymmSCMatrix MBPT2::density() { return 0; } ////////////////////////////////////////////////////////////////////////////// void MBPT2::compute() { if (std::string(reference_->integral()->class_name()) !=integral()->class_name()) { FeatureNotImplemented ex( "cannot use a reference with a different Integral specialization", __FILE__, __LINE__, class_desc()); try { ex.elaborate() << "reference uses " << reference_->integral()->class_name() << " but this object uses " << integral()->class_name() << std::endl; } catch (...) {} throw ex; } init_variables(); reference_->set_desired_value_accuracy(desired_value_accuracy() / ref_to_mp2_acc); if (gradient_needed()) { if (nsocc) { ExEnv::errn() << "MBPT2: cannot compute open shell gradients" << endl; abort(); } compute_cs_grad(); } else { if (nsocc && algorithm_ && !strcmp(algorithm_,"memgrp")) { ExEnv::errn() << "MBPT2: memgrp algorithm cannot compute open shell energy" << endl; abort(); } if (!nsocc && (!algorithm_ || !strcmp(algorithm_,"memgrp"))) { compute_cs_grad(); } else if ((!algorithm_ && msg_->n() <= 32) || (algorithm_ && !strcmp(algorithm_,"v1"))) { compute_hsos_v1(); } else if (!algorithm_ || !strcmp(algorithm_,"v2")) { compute_hsos_v2(); } else if (!strcmp(algorithm_,"v2lb")) { compute_hsos_v2_lb(); } else { ExEnv::errn() << "MBPT2: unknown algorithm: " << algorithm_ << endl; abort(); } } } ////////////////////////////////////////////////////////////////////////////// void MBPT2::obsolete() { // Solaris 2.7 workshop 5.0 is causing this routine to // be incorrectly called in a base class CTOR. Thus // reference_ might be null and it must be tested. if (reference_.nonnull()) reference_->obsolete(); Wavefunction::obsolete(); } ////////////////////////////////////////////////////////////////////////////// int MBPT2::gradient_implemented() const { int nb = reference_->oso_dimension()->n(); int n = 0; for (int i=0; ioccupation(i) == 1.0) n++; } if (n) return 0; return 1; } ////////////////////////////////////////////////////////////////////////////// int MBPT2::value_implemented() const { return 1; } ////////////////////////////////////////////////////////////////////////////// void MBPT2::eigen(RefDiagSCMatrix &vals, RefSCMatrix &vecs, RefDiagSCMatrix &occs) { int i, j; if (nsocc) { if (reference_->n_fock_matrices() != 2) { ExEnv::errn() << "MBPT2: given open reference with" << " wrong number of Fock matrices" << endl; abort(); } // Notation: oo = orthonormal symmetry orbital basis // ao = atomic orbital basis // so = symmetrized atomic orbital basis // mo1 = SCF molecular orbital basis // mo2 = MBPT molecular orbital basis // get the closed shell and open shell AO fock matrices RefSymmSCMatrix fock_c_so = reference_->fock(0); RefSymmSCMatrix fock_o_so = reference_->fock(1); // transform the AO fock matrices to the MO basis RefSymmSCMatrix fock_c_mo1 = basis_matrixkit()->symmmatrix(oso_dimension()); RefSymmSCMatrix fock_o_mo1 = basis_matrixkit()->symmmatrix(oso_dimension()); RefSCMatrix vecs_so_mo1 = reference_->eigenvectors(); fock_c_mo1.assign(0.0); fock_o_mo1.assign(0.0); fock_c_mo1.accumulate_transform(vecs_so_mo1.t(), fock_c_so); fock_o_mo1.accumulate_transform(vecs_so_mo1.t(), fock_o_so); fock_c_so = 0; fock_o_so = 0; /* Convert to the Guest & Saunders general form. This is the form used for an OPT2 calculation. C O V ---------- | | C | fc | | | ------------------- | | | O | 2fc-fo | fc | | | | ---------------------------- | | | | V | fc | fo | fc | | | | | ---------------------------- */ RefSymmSCMatrix fock_eff_mo1 = fock_c_mo1.clone(); fock_eff_mo1.assign(fock_c_mo1); for (i=0; in(); i++) { double occi = reference_->occupation(i); for (j=0; j<=i; j++) { double occj = reference_->occupation(j); if (occi == 2.0 && occj == 1.0 || occi == 1.0 && occj == 2.0) { fock_eff_mo1.accumulate_element(i,j, fock_c_mo1(i,j)-fock_o_mo1(i,j)); } else if (occi == 0.0 && occj == 1.0 || occi == 1.0 && occj == 0.0) { fock_eff_mo1.accumulate_element(i,j, fock_o_mo1(i,j)-fock_c_mo1(i,j)); } } } // diagonalize the new fock matrix RefDiagSCMatrix vals_mo2(fock_eff_mo1.dim(), fock_eff_mo1.kit()); RefSCMatrix vecs_mo1_mo2(fock_eff_mo1.dim(), fock_eff_mo1.dim(), fock_eff_mo1.kit()); fock_eff_mo1.diagonalize(vals_mo2, vecs_mo1_mo2); vals = vals_mo2; // compute the AO to new MO scf vector RefSCMatrix so_ao = reference_->integral()->petite_list()->sotoao(); if (debug_ > 1) { vecs_mo1_mo2.t().print("vecs_mo1_mo2.t()"); vecs_so_mo1.t().print("vecs_so_mo1.t()"); so_ao.print("so_ao"); } vecs = vecs_mo1_mo2.t() * vecs_so_mo1.t() * so_ao; } else { if (debug_) ExEnv::out0() << indent << "getting fock matrix" << endl; // get the closed shell AO fock matrices RefSymmSCMatrix fock_c_so = reference_->fock(0); // transform the AO fock matrices to the MO basis RefSymmSCMatrix fock_c_mo1 = basis_matrixkit()->symmmatrix(oso_dimension()); RefSCMatrix vecs_so_mo1 = reference_->eigenvectors(); fock_c_mo1.assign(0.0); fock_c_mo1.accumulate_transform(vecs_so_mo1.t(), fock_c_so); fock_c_so = 0; if (debug_) ExEnv::out0() << indent << "diagonalizing" << endl; // diagonalize the fock matrix vals = fock_c_mo1.eigvals(); // compute the AO to new MO scf vector if (debug_) ExEnv::out0() << indent << "AO to MO" << endl; RefSCMatrix so_ao = reference_->integral()->petite_list()->sotoao(); vecs = vecs_so_mo1.t() * so_ao; } // fill in the occupations occs = matrixkit()->diagmatrix(vals.dim()); for (i=0; in(); i++) { occs(i) = reference_->occupation(i); } // allocate storage for symmetry information if (!symorb_irrep_) symorb_irrep_ = new int[nbasis]; if (!symorb_num_) symorb_num_ = new int[nbasis]; // Check for degenerate eigenvalues. Use unsorted eigenvalues since it // only matters if the degeneracies occur within a given irrep. BlockedDiagSCMatrix *bvals = dynamic_cast(vals.pointer()); for (i=0; inblocks(); i++) { int done = 0; RefDiagSCMatrix valsi = bvals->block(i); for (j=1; jpoint_group()->char_table().order() != 1) { if (debug_) ExEnv::out0() << indent << "sorting eigenvectors" << endl; double *evals = new double[noso]; vals->convert(evals); int *indices = new int[noso]; dquicksort(evals,indices,noso); delete[] evals; // make sure all nodes see the same indices and evals msg_->bcast(indices,noso); RefSCMatrix newvecs(vecs.rowdim(), vecs.coldim(), matrixkit()); RefDiagSCMatrix newvals(vals.dim(), matrixkit()); RefDiagSCMatrix newoccs(vals.dim(), matrixkit()); for (i=0; ipoint_group()->char_table(); int orbnum = 0; int *tmp_irrep = new int[noso]; int *tmp_num = new int[noso]; for (i=0; iblocks()->nblock(); i++) { for (j=0; jblocks()->size(i); j++, orbnum++) { tmp_irrep[orbnum] = i; tmp_num[orbnum] = j; } } for (i=0; i= 0) { double split = vals(nbasis-nfzv) - vals(nbasis-nfzv-1); if (split < 0.2) { ExEnv::out0() << endl << indent << "WARNING: " << "MBPT2: gap between frozen and active virtual orbitals is " << split << " au" << endl << endl; } } if (debug_) ExEnv::out0() << indent << "eigen done" << endl; } ///////////////////////////////////////////////////////////////////////////// void MBPT2::init_variables() { nbasis = so_dimension()->n(); noso = oso_dimension()->n(); // if (nbasis != noso) { // ExEnv::outn() << "MBPT2: Noso != Nbasis: MBPT2 not checked for this case" << endl; // abort(); // } nocc = nvir = nsocc = 0; for (int i=0; ioccupation(i) == 2.0) nocc++; else if (reference_->occupation(i) == 1.0) nsocc++; else nvir++; } } ///////////////////////////////////////////////////////////////////////////// void MBPT2::symmetry_changed() { Wavefunction::symmetry_changed(); reference_->symmetry_changed(); } ///////////////////////////////////////////////////////////////////////////// int MBPT2::nelectron() { return reference_->nelectron(); } ///////////////////////////////////////////////////////////////////////////// double MBPT2::ref_energy() { return reference_->energy(); } ///////////////////////////////////////////////////////////////////////////// double MBPT2::corr_energy() { return energy() - reference_->energy(); } ///////////////////////////////////////////////////////////////////////////// RefSCVector MBPT2::ref_energy_gradient() { gradient(); return hf_gradient_; } ///////////////////////////////////////////////////////////////////////////// RefSCVector MBPT2::corr_energy_gradient() { gradient(); return get_cartesian_gradient() - hf_gradient_; } ///////////////////////////////////////////////////////////////////////////// // Local Variables: // mode: c++ // c-file-style: "CLJ" // End: