Annotation of /home/cljanss/src/SC/src/lib/chemistry/qc/scf/scf.cc for ./mpqc.vmon.0018

    1                                                             //
    2                                                             // scf.cc --- implementation of the SCF abstract base class
    3                                                             //
    4                                                             // Copyright (C) 1996 Limit Point Systems, Inc.
    5                                                             //
    6                                                             // Author: Edward Seidl <seidl@janed.com>
    7                                                             // Maintainer: LPS
    8                                                             //
    9                                                             // This file is part of the SC Toolkit.
   10                                                             //
   11                                                             // The SC Toolkit is free software; you can redistribute it and/or modify
   12                                                             // it under the terms of the GNU Library General Public License as published by
   13                                                             // the Free Software Foundation; either version 2, or (at your option)
   14                                                             // any later version.
   15                                                             //
   16                                                             // The SC Toolkit is distributed in the hope that it will be useful,
   17                                                             // but WITHOUT ANY WARRANTY; without even the implied warranty of
   18                                                             // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   19                                                             // GNU Library General Public License for more details.
   20                                                             //
   21                                                             // You should have received a copy of the GNU Library General Public License
   22                                                             // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
   23                                                             // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
   24                                                             //
   25                                                             // The U.S. Government is granted a limited license as per AL 91-7.
   26                                                             //
   27                                                             
   28                                                             #ifdef __GNUC__
   29                                                             #pragma implementation
   30                                                             #endif
   31                                                             
   32                                                             #include <math.h>
   33                                                             #include <limits.h>
   34                                                             #include <sys/stat.h>
   35                                                             
   36                                                             #include <util/misc/formio.h>
   37                                                             #include <util/state/stateio.h>
   38                                                             #include <util/group/mstate.h>
   39                                                             
   40                                                             #include <math/scmat/local.h>
   41                                                             #include <math/scmat/repl.h>
   42                                                             #include <math/scmat/offset.h>
   43                                                             #include <math/scmat/blocked.h>
   44                                                             
   45                                                             #include <math/optimize/diis.h>
   46                                                             
   47                                                             #include <chemistry/qc/basis/petite.h>
   48                                                             #include <chemistry/qc/scf/scf.h>
   49                                                             
   50                                                             using namespace std;
   51                                                             
   52                                                             ///////////////////////////////////////////////////////////////////////////
   53                                                             // SCF
   54                                                             
   55                                                             static ClassDesc SCF_cd(
   56                                                               typeid(SCF),"SCF",2,"public OneBodyWavefunction",
   57                                                               0, 0, 0);
   58                                                             
   59                                                             SCF::SCF(StateIn& s) :
   60                                                               SavableState(s),
   61                                                               OneBodyWavefunction(s)
   62                                                             {
   63                                                               need_vec_ = 1;
   64                                                               compute_guess_ = 0;
   65                                                             
   66                                                               s.get(maxiter_,"maxiter");
   67                                                               s.get(dens_reset_freq_);
   68                                                               s.get(reset_occ_);
   69                                                               s.get(local_dens_);
   70                                                               s.get(storage_);
   71                                                               if (s.version(::class_desc<SCF>()) >= 2) {
   72                                                                 s.get(print_all_evals_);
   73                                                                 s.get(print_occ_evals_);
   74                                                               }
   75                                                               s.get(level_shift_);
   76                                                             
   77                                                               extrap_ << SavableState::restore_state(s);
   78                                                               accumdih_ << SavableState::restore_state(s);
   79                                                               accumddh_ << SavableState::restore_state(s);
   80                                                             
   81                                                               scf_grp_ = basis()->matrixkit()->messagegrp();
   82                                                               threadgrp_ = ThreadGrp::get_default_threadgrp();
   83                                                             }
   84                                                             
   85                                                             SCF::SCF(const Ref<KeyVal>& keyval) :
   86                                                               OneBodyWavefunction(keyval),
   87                                                               need_vec_(1),
   88                                                               compute_guess_(0),
   89                                                               maxiter_(40),
   90                                                               dens_reset_freq_(10),
   91                                                               reset_occ_(0),
   92                                                               local_dens_(1),
   93                                                               storage_(0),
   94                                                               level_shift_(0)
   95                                                             {
   96                                                               if (keyval->exists("maxiter"))
   97                                                                 maxiter_ = keyval->intvalue("maxiter");
   98                                                             
   99                                                               if (keyval->exists("density_reset_frequency"))
  100                                                                 dens_reset_freq_ = keyval->intvalue("density_reset_frequency");
  101                                                             
  102                                                               if (keyval->exists("reset_occupations"))
  103                                                                 reset_occ_ = keyval->booleanvalue("reset_occupations");
  104                                                             
  105                                                               if (keyval->exists("level_shift"))
  106                                                                 level_shift_ = keyval->doublevalue("level_shift");
  107                                                             
  108                                                               extrap_ << keyval->describedclassvalue("extrap");
  109                                                               if (extrap_.null())
  110                                                                 extrap_ = new DIIS;
  111                                                             
  112                                                               accumdih_ << keyval->describedclassvalue("accumdih");
  113                                                               if (accumdih_.null())
  114                                                                 accumdih_ = new AccumHNull;
  115                                                               
  116                                                               accumddh_ << keyval->describedclassvalue("accumddh");
  117                                                               if (accumddh_.null())
  118                                                                 accumddh_ = new AccumHNull;
  119                                                               
  120                                                               KeyValValueint defaultmem(8000000);
  121                                                               storage_ = keyval->intvalue("memory",defaultmem);
  122                                                               
  123                                                               if (keyval->exists("local_density"))
  124                                                                 local_dens_ = keyval->booleanvalue("local_density");
  125                                                                 
  126                                                               print_all_evals_ = keyval->booleanvalue("print_evals");
  127                                                               print_occ_evals_ = keyval->booleanvalue("print_occupied_evals");
  128                                                               
  129                                                               scf_grp_ = basis()->matrixkit()->messagegrp();
  130                                                               threadgrp_ = ThreadGrp::get_default_threadgrp();
  131                                                               
  132                                                               // first see if guess_wavefunction is a wavefunction, then check to
  133                                                               // see if it's a string.
  134                                                               if (keyval->exists("guess_wavefunction")) {
  135                                                                 ExEnv::out() << incindent << incindent;
  136                                                                 guess_wfn_ << keyval->describedclassvalue("guess_wavefunction");
  137                                                                 compute_guess_=1;
  138                                                                 if (guess_wfn_.null()) {
  139                                                                   compute_guess_=0;
  140                                                                   char *path = keyval->pcharvalue("guess_wavefunction");
  141                                                                   struct stat sb;
  142                                                                   if (path && stat(path, &sb)==0 && sb.st_size) {
  143                                                                     BcastStateInBin s(scf_grp_, path);
  144                                                             
  145                                                                     // reset the default matrixkit so that the matrices in the guess
  146                                                                     // wavefunction will match those in this wavefunction
  147                                                                     Ref<SCMatrixKit> oldkit = SCMatrixKit::default_matrixkit();
  148                                                                     SCMatrixKit::set_default_matrixkit(basis()->matrixkit());
  149                                                             
  150                                                                     guess_wfn_ << SavableState::restore_state(s);
  151                                                             
  152                                                                     // go back to the original default matrixkit
  153                                                                     SCMatrixKit::set_default_matrixkit(oldkit);
  154                                                                     delete[] path;
  155                                                                   }
  156                                                                 }
  157                                                                 ExEnv::out() << decindent << decindent;
  158                                                               }
  159                                                             }
  160                                                             
  161                                                             SCF::~SCF()
  162                                                             {
  163                                                             }
  164                                                             
  165                                                             void
  166                                                             SCF::save_data_state(StateOut& s)
  167                                                             {
  168                                                               OneBodyWavefunction::save_data_state(s);
  169                                                               s.put(maxiter_);
  170                                                               s.put(dens_reset_freq_);
  171                                                               s.put(reset_occ_);
  172                                                               s.put(local_dens_);
  173                                                               s.put(storage_);
  174                                                               s.put(print_all_evals_);
  175                                                               s.put(print_occ_evals_);
  176                                                               s.put(level_shift_);
  177                                                               SavableState::save_state(extrap_.pointer(),s);
  178                                                               SavableState::save_state(accumdih_.pointer(),s);
  179                                                               SavableState::save_state(accumddh_.pointer(),s);
  180                                                             }
  181                                                             
  182                                                             RefSCMatrix
  183                                                             SCF::oso_eigenvectors()
  184                                                             {
  185                                                               return oso_eigenvectors_.result();
  186                                                             }
  187                                                             
  188                                                             RefDiagSCMatrix
  189                                                             SCF::eigenvalues()
  190                                                             {
  191                                                               return eigenvalues_.result();
  192                                                             }
  193                                                             
  194                                                             int
  195                                                             SCF::spin_unrestricted()
  196                                                             {
  197                                                               return 0;
  198                                                             }
  199                                                             
  200                                                             void
  201                                                             SCF::print(ostream&o) const
  202                                                             {
  203                                                               OneBodyWavefunction::print(o);
  204                                                               o << node0 << indent << "SCF Parameters:\n" << incindent
  205                                                                 << indent << "maxiter = " << maxiter_ << endl
  206                                                                 << indent << "density_reset_frequency = " << dens_reset_freq_ << endl
  207                                                                 << indent << scprintf("level_shift = %f\n",level_shift_)
  208                                                                 << decindent << endl;
  209                                                             }
  210                                                             
  211                                                             //////////////////////////////////////////////////////////////////////////////
  212                                                             
  213                                                             void
  214                                                             SCF::compute()
  215                                                             {
  216                                                               local_ = (dynamic_cast<LocalSCMatrixKit*>(basis()->matrixkit().pointer()) ||
  217                                                                         dynamic_cast<ReplSCMatrixKit*>(basis()->matrixkit().pointer())) ? 1:0;
  218                                                               
  219                                                               const double hess_to_grad_acc = 1.0/100.0;
  220                                                               if (hessian_needed())
  221                                                                 set_desired_gradient_accuracy(desired_hessian_accuracy()*hess_to_grad_acc);
  222                                                             
  223                                                               const double grad_to_val_acc = 1.0/100.0;
  224                                                               if (gradient_needed())
  225                                                                 set_desired_value_accuracy(desired_gradient_accuracy()*grad_to_val_acc);
  226                                                             
  227                                                               double delta;
  228                                                               if (value_needed()) {
  229                                                                 ExEnv::out() << node0 << endl << indent
  230                                                                      << scprintf("SCF::compute: energy accuracy = %10.7e\n",
  231                                                                                  desired_value_accuracy())
  232                                                                      << endl;
  233                                                             
  234                                                                 double eelec;
  235                                                                 delta = compute_vector(eelec);
  236                                                                   
  237                                                                 // this will be done elsewhere eventually
  238                                                                 double nucrep = molecule()->nuclear_repulsion_energy();
  239                                                                 double eother = 0.0;
  240                                                                 if (accumddh_.nonnull()) eother = accumddh_->e();
  241                                                                 ExEnv::out() << node0 << endl << indent
  242                                                                      << scprintf("total scf energy = %15.10f", eelec+eother+nucrep)
  243                                                                      << endl;
  244                                                             
  245                                                                 set_energy(eelec+eother+nucrep);
  246                                                                 set_actual_value_accuracy(delta);
  247                                                               }
  248                                                               else {
  249                                                                 delta = actual_value_accuracy();
  250                                                               }
  251                                                             
  252                                                               if (gradient_needed()) {
  253                                                                 RefSCVector gradient = matrixkit()->vector(moldim());
  254                                                             
  255                                                                 ExEnv::out() << node0 << endl << indent
  256                                                                      << scprintf("SCF::compute: gradient accuracy = %10.7e\n",
  257                                                                                  desired_gradient_accuracy())
  258                                                                      << endl;
  259                                                             
  260                                                                 compute_gradient(gradient);
  261                                                                 print_natom_3(gradient,"Total Gradient:");
  262                                                                 set_gradient(gradient);
  263                                                             
  264                                                                 set_actual_gradient_accuracy(delta/grad_to_val_acc);
  265                                                               }
  266                                                               
  267                                                               if (hessian_needed()) {
  268                                                                 RefSymmSCMatrix hessian = matrixkit()->symmmatrix(moldim());
  269                                                                 
  270                                                                 ExEnv::out() << node0 << endl << indent
  271                                                                      << scprintf("SCF::compute: hessian accuracy = %10.7e\n",
  272                                                                                  desired_hessian_accuracy())
  273                                                                      << endl;
  274                                                             
  275                                                                 compute_hessian(hessian);
  276                                                                 set_hessian(hessian);
  277                                                             
  278                                                                 set_actual_hessian_accuracy(delta/grad_to_val_acc/hess_to_grad_acc);
  279                                                               }
  280                                                             }
  281                                                             
  282                                                             //////////////////////////////////////////////////////////////////////////////
  283                                                             
  284                                                             signed char *
  285                                                             SCF::init_pmax(double *pmat_data)
  286                                                             {
  287                                                               double l2inv = 1.0/log(2.0);
  288                                                               double tol = pow(2.0,-126.0);
  289                                                               
  290                                                               GaussianBasisSet& gbs = *basis().pointer();
  291                                                               
  292                                                               signed char * pmax = new signed char[i_offset(gbs.nshell())];
  293                                                             
  294                                                               int ish, jsh, ij;
  295                                                               for (ish=ij=0; ish < gbs.nshell(); ish++) {
  296                                                                 int istart = gbs.shell_to_function(ish);
  297                                                                 int iend = istart + gbs(ish).nfunction();
  298                                                                 
  299     1     0     0     0     1     0     0     0     0     0     for (jsh=0; jsh <= ish; jsh++,ij++) {
  300                                                                   int jstart = gbs.shell_to_function(jsh);
  301                                                                   int jend = jstart + gbs(jsh).nfunction();
  302                                                                   
  303                                                                   double maxp=0, tmp;
  304                                                             
  305                                                                   for (int i=istart; i < iend; i++) {
  306     3     0     0     0     0     0     0     0     0     0         int ijoff = i_offset(i) + jstart;
  307     2     0     0     0     1     0     0     0     0     0         for (int j=jstart; j < ((ish==jsh) ? i+1 : jend); j++,ijoff++)
  308                                                                       if ((tmp=fabs(pmat_data[ijoff])) > maxp)
  309     2     0     0     0     3     1     0     0     0     0             maxp=tmp;
  310    17     0     9     2     9     3     0     0     0     0       }
  311                                                             
  312     2     0     2     0     1     0     0     0     0     0       if (maxp <= tol)
  313     1     0     0     0     0     0     0     0     0     0         maxp=tol;
  314                                                             
  315                                                                   long power = long(ceil(log(maxp)*l2inv));
  316     1     0     0     0     0     0     0     0     0     0       if (power < SCHAR_MIN) pmax[ij] = SCHAR_MIN;
  317                                                                   else if (power > SCHAR_MAX) pmax[ij] = SCHAR_MAX;
  318                                                                   else pmax[ij] = (signed char) power;
  319     1     0     0     0     1     0     0     0     0     0     }
  320                                                               }
  321                                                             
  322                                                               return pmax;
  323                                                             }
  324                                                             
  325                                                             //////////////////////////////////////////////////////////////////////////////
  326                                                             
  327                                                             RefSymmSCMatrix
  328                                                             SCF::get_local_data(const RefSymmSCMatrix& m, double*& p, Access access)
  329     1     0     1     0     0     0     0     0     0     0 {
  330                                                               RefSymmSCMatrix l = m;
  331                                                               
  332                                                               if (!dynamic_cast<LocalSymmSCMatrix*>(l.pointer())
  333                                                                   && !dynamic_cast<ReplSymmSCMatrix*>(l.pointer())) {
  334                                                                 Ref<SCMatrixKit> k = new ReplSCMatrixKit;
  335                                                                 l = k->symmmatrix(m.dim());
  336                                                                 l->convert(m);
  337                                                             
  338                                                                 if (access == Accum)
  339                                                                   l->assign(0.0);
  340                                                               } else if (scf_grp_->n() > 1 && access==Accum) {
  341                                                                 l = m.clone();
  342                                                                 l.assign(0.0);
  343                                                               }
  344                                                             
  345                                                               if (dynamic_cast<ReplSymmSCMatrix*>(l.pointer()))
  346                                                                 p = dynamic_cast<ReplSymmSCMatrix*>(l.pointer())->get_data();
  347                                                               else
  348                                                                 p = dynamic_cast<LocalSymmSCMatrix*>(l.pointer())->get_data();
  349                                                             
  350                                                               return l;
  351                                                             }
  352                                                             
  353                                                             //////////////////////////////////////////////////////////////////////////////
  354                                                             
  355                                                             void
  356                                                             SCF::initial_vector(int needv)
  357                                                             {
  358     1     0     0     0     0     0     0     0     0     0   if (need_vec_) {
  359                                                                 if (oso_eigenvectors_.result_noupdate().null()) {
  360                                                                   // if guess_wfn_ is non-null then try to get a guess vector from it.
  361                                                                   // First check that the same basis is used...if not, then project the
  362                                                                   // guess vector into the present basis.
  363                                                                   if (guess_wfn_.nonnull()) {
  364                                                                     if (basis()->equiv(guess_wfn_->basis())
  365                                                                         &&orthog_method() == guess_wfn_->orthog_method()) {
  366                                                                       ExEnv::out() << node0 << indent
  367                                                                            << "Using guess wavefunction as starting vector" << endl;
  368                                                             
  369                                                                       // indent output of eigenvectors() call if there is any
  370                                                                       ExEnv::out() << incindent << incindent;
  371                                                                       SCF *g = dynamic_cast<SCF*>(guess_wfn_.pointer());
  372                                                                       if (!g || compute_guess_) {
  373                                                                         oso_eigenvectors_ = guess_wfn_->oso_eigenvectors().copy();
  374                                                                         eigenvalues_ = guess_wfn_->eigenvalues().copy();
  375                                                                       } else {
  376                                                                         oso_eigenvectors_ = g->oso_eigenvectors_.result_noupdate().copy();
  377                                                                         eigenvalues_ = g->eigenvalues_.result_noupdate().copy();
  378                                                                       }
  379                                                                       ExEnv::out() << decindent << decindent;
  380                                                                     } else {
  381                                                                       ExEnv::out() << node0 << indent
  382                                                                            << "Projecting guess wavefunction into the present basis set"
  383                                                                            << endl;
  384                                                             
  385                                                                       // indent output of projected_eigenvectors() call if there is any
  386                                                                       ExEnv::out() << incindent << incindent;
  387                                                                       oso_eigenvectors_ = projected_eigenvectors(guess_wfn_);
  388                                                                       eigenvalues_ = projected_eigenvalues(guess_wfn_);
  389                                                                       ExEnv::out() << decindent << decindent;
  390                                                                     }
  391                                                             
  392                                                                     // we should only have to do this once, so free up memory used
  393                                                                     // for the old wavefunction
  394                                                                     guess_wfn_=0;
  395                                                             
  396                                                                     ExEnv::out() << node0 << endl;
  397                                                                   
  398                                                                   } else {
  399                                                                     ExEnv::out() << node0 << indent << "Starting from core Hamiltonian guess\n"
  400                                                                          << endl;
  401                                                                     oso_eigenvectors_ = hcore_guess(eigenvalues_.result_noupdate());
  402                                                                   }
  403                                                                 } else {
  404                                                                   // this is just an old vector
  405                                                                 }
  406                                                               }
  407                                                             
  408                                                               need_vec_=needv;
  409                                                             }
  410                                                             
  411                                                             //////////////////////////////////////////////////////////////////////////////
  412                                                             
  413                                                             void
  414                                                             SCF::init_mem(int nm)
  415                                                             {
  416                                                               // if local_den_ is already 0, then that means it was set to zero by
  417                                                               // the user.
  418                                                               if (!local_dens_) {
  419                                                                 integral()->set_storage(storage_);
  420                                                                 return;
  421                                                               }
  422                                                               
  423                                                               int nmem = i_offset(basis()->nbasis())*nm*sizeof(double);
  424                                                             
  425                                                               // if we're actually using local matrices, then there's no choice
  426                                                               if (dynamic_cast<LocalSCMatrixKit*>(basis()->matrixkit().pointer())
  427                                                                   ||dynamic_cast<ReplSCMatrixKit*>(basis()->matrixkit().pointer())) {
  428                                                                 if (nmem > storage_)
  429                                                                   return;
  430                                                               } else {
  431                                                                 if (nmem > storage_) {
  432                                                                   local_dens_=0;
  433                                                                   integral()->set_storage(storage_);
  434                                                                   return;
  435                                                                 }
  436                                                               }
  437                                                             
  438                                                               integral()->set_storage(storage_-nmem);
  439                                                             }
  440                                                             
  441                                                             /////////////////////////////////////////////////////////////////////////////
  442                                                             
  443                                                             void
  444                                                             SCF::so_density(const RefSymmSCMatrix& d, double occ, int alp)
  445                                                             {
  446                                                               int i,j,k;
  447                                                               int me=scf_grp_->me();
  448                                                               int nproc=scf_grp_->n();
  449                                                               int uhf = spin_unrestricted();
  450                                                               
  451                                                               d->assign(0.0);
  452                                                               
  453                                                               RefSCMatrix oso_vector;
  454                                                               if (alp || !uhf) {
  455                                                                 if (oso_scf_vector_.nonnull())
  456                                                                   oso_vector = oso_scf_vector_;
  457                                                               }
  458                                                               else {
  459                                                                 if (oso_scf_vector_beta_.nonnull())
  460                                                                   oso_vector = oso_scf_vector_beta_;
  461                                                               }
  462                                                                   
  463                                                               if (oso_vector.null()) {
  464                                                                 if (uhf) {
  465                                                                   if (alp)
  466                                                                     oso_vector = oso_alpha_eigenvectors();
  467                                                                   else
  468                                                                     oso_vector = oso_beta_eigenvectors();
  469                                                                 } else
  470                                                                   oso_vector = oso_eigenvectors();
  471                                                               }
  472                                                             
  473     1     0     0     0     0     0     0     0     0     0   if (debug_ > 1) oso_vector.print("ortho SO vector");
  474                                                             
  475                                                               RefSCMatrix vector = so_to_orthog_so().t() * oso_vector;
  476                                                               oso_vector = 0;
  477                                                             
  478                                                               if (debug_ > 1) vector.print("SO vector");
  479                                                               
  480                                                               BlockedSCMatrix *bvec = require_dynamic_cast<BlockedSCMatrix*>(
  481                                                                 vector, "SCF::so_density: blocked vector");
  482                                                             
  483                                                               BlockedSymmSCMatrix *bd = require_dynamic_cast<BlockedSymmSCMatrix*>(
  484                                                                 d, "SCF::so_density: blocked density");
  485                                                               
  486                                                               for (int ir=0; ir < oso_dimension()->blocks()->nblock(); ir++) {
  487                                                                 RefSCMatrix vir = bvec->block(ir);
  488                                                                 RefSymmSCMatrix dir = bd->block(ir);
  489                                                                 
  490                                                                 if (vir.null() || vir.ncol()==0)
  491                                                                   continue;
  492                                                                 
  493                                                                 int n_orthoSO = oso_dimension()->blocks()->size(ir);
  494                                                                 int n_SO = so_dimension()->blocks()->size(ir);
  495                                                                 
  496                                                                 // figure out which columns of the scf vector we'll need
  497                                                                 int col0 = -1, coln = -1;
  498                                                                 for (i=0; i < n_orthoSO; i++) {
  499                                                                   double occi;
  500     1     0     0     0     0     0     0     0     0     0       if (!uhf)
  501     0     0     1     0     0     0     0     0     0     0         occi = occupation(ir, i);
  502                                                                   else if (alp)
  503                                                                     occi = alpha_occupation(ir, i);
  504                                                                   else
  505                                                                     occi = beta_occupation(ir, i);
  506                                                             
  507                                                                   if (fabs(occi-occ) < 1.0e-8) {
  508                                                                     if (col0 == -1)
  509                                                                       col0 = i;
  510                                                                     continue;
  511     0     0     1     0     0     0     0     0     0     0       } else if (col0 != -1) {
  512                                                                     coln = i-1;
  513                                                                     break;
  514                                                                   }
  515                                                                 }
  516                                                             
  517                                                                 if (col0 == -1)
  518                                                                   continue;
  519                                                             
  520     1     0     0     0     0     0     0     0     0     0     if (coln == -1)
  521                                                                   coln = n_orthoSO-1;
  522                                                                 
  523     1     0     0     0     0     0     0     0     0     0     if (local_ || local_dens_) {
  524                                                                   RefSymmSCMatrix ldir = dir;
  525                                                             
  526                                                                   RefSCMatrix occbits; // holds the occupied bits of the scf vector
  527                                                             
  528                                                                   // get local copies of vector and density matrix
  529                                                                   if (!local_) {
  530                                                                     Ref<SCMatrixKit> rk = new ReplSCMatrixKit;
  531                                                                     RefSCMatrix lvir = rk->matrix(vir.rowdim(), vir.coldim());
  532                                                                     lvir->convert(vir);
  533                                                                     occbits = lvir->get_subblock(0, n_SO-1, col0, coln);
  534                                                                     lvir = 0;
  535                                                             
  536                                                                     ldir = rk->symmmatrix(dir.dim());
  537                                                                     ldir->convert(dir);
  538                                                             
  539                                                                   } else {
  540                                                                     occbits = vir->get_subblock(0, n_SO-1, col0, coln);
  541                                                                   }
  542                                                                 
  543                                                                   double **c;
  544                                                                   double *dens;
  545                                                                   
  546                                                                   if (dynamic_cast<LocalSCMatrix*>(occbits.pointer()))
  547                                                                     c = dynamic_cast<LocalSCMatrix*>(occbits.pointer())->get_rows();
  548                                                                   else if (dynamic_cast<ReplSCMatrix*>(occbits.pointer()))
  549                                                                     c = dynamic_cast<ReplSCMatrix*>(occbits.pointer())->get_rows();
  550                                                                   else
  551                                                                     abort();
  552                                                             
  553                                                                   if (dynamic_cast<LocalSymmSCMatrix*>(ldir.pointer()))
  554                                                                     dens = dynamic_cast<LocalSymmSCMatrix*>(ldir.pointer())->get_data();
  555                                                                   else if (dynamic_cast<ReplSymmSCMatrix*>(ldir.pointer()))
  556                                                                     dens = dynamic_cast<ReplSymmSCMatrix*>(ldir.pointer())->get_data();
  557                                                                   else
  558                                                                     abort();
  559                                                             
  560                                                                   int ij=0;
  561                                                                   for (i=0; i < n_SO; i++) {
  562                                                                     for (j=0; j <= i; j++, ij++) {
  563    17     2    11     8     9     0     0     0     0     0           if (ij%nproc != me)
  564                                                                         continue;
  565                                                                       
  566                                                                       double dv = 0;
  567                                                             
  568     0     0     0     1     1     0     0     0     0     0           int kk=0;
  569     1     0     3     0     0     0     0     0     1     0           for (k=col0; k <= coln; k++, kk++)
  570    11     1    10     6     6     1     0     0     0     0             dv += c[i][kk]*c[j][kk];
  571                                                             
  572     1     0     0     0     0     0     0     0     0     0           dens[ij] = dv;
  573     3     0     0     1     0     0     0     0     0     0         }
  574                                                                   }
  575                                                             
  576                                                                   if (nproc > 1)
  577                                                                     scf_grp_->sum(dens, i_offset(n_SO));
  578                                                             
  579                                                                   if (!local_) {
  580                                                                     dir->convert(ldir);
  581                                                                   }
  582                                                                 }
  583                                                             
  584                                                                 // for now quit
  585                                                                 else {
  586                                                                   ExEnv::err() << node0 << indent
  587                                                                        << "Cannot yet use anything but Local matrices"
  588                                                                        << endl;
  589                                                                   abort();
  590                                                                 }
  591                                                               }
  592                                                             
  593                                                               if (debug_ > 0) {
  594                                                                 ExEnv::out() << node0 << indent
  595                                                                      << "Nelectron = " << 2.0 * (d * overlap()).trace() << endl;
  596                                                               }
  597                                                             
  598                                                               int use_alternate_density = 0;
  599                                                               if (use_alternate_density || debug_ > 2) {
  600                                                                 // double check the density with this simpler, slower way to compute
  601                                                                 // the density matrix
  602                                                                 RefSymmSCMatrix occ(oso_dimension(), basis_matrixkit());
  603                                                                 occ.assign(0.0);
  604                                                                 for (i=0; i<oso_dimension()->n(); i++) {
  605                                                                   occ(i,i) = occupation(i);
  606                                                                 }
  607                                                                 occ.scale(0.5);
  608                                                                 RefSymmSCMatrix d2(so_dimension(), basis_matrixkit());
  609                                                                 d2.assign(0.0);
  610                                                                 d2.accumulate_transform(vector, occ);
  611                                                                 if (debug_ > 2) {
  612                                                                   d2.print("d2 density");
  613                                                                   ExEnv::out() << node0 << indent << "d2 Nelectron = "
  614                                                                                << 2.0 * (d2 * overlap()).trace() << endl;
  615                                                                 }
  616                                                                 if (use_alternate_density) {
  617                                                                   d.assign(d2);
  618                                                                   ExEnv::out() << node0 << indent
  619                                                                                << "using alternate density algorithm" << endl;
  620                                                                 }
  621                                                               }
  622                                                             
  623     1     0     0     0     0     0     0     0     0     0   if (debug_ > 1) {
  624                                                                 d.print("SO Density");
  625                                                                 RefSCMatrix rd(d.dim(), d.dim(), basis_matrixkit());
  626                                                                 rd.assign(0.0);
  627                                                                 rd.accumulate(d);
  628                                                                 (d*overlap()*d-rd).print("SO Density idempotency error");
  629                                                               }
  630     0     0     0     1     0     0     0     0     0     0 }
  631                                                             
  632                                                             double
  633                                                             SCF::one_body_energy()
  634                                                             {
  635                                                               RefSymmSCMatrix dens = ao_density().copy();
  636                                                               RefSymmSCMatrix hcore = dens->clone();
  637                                                               hcore.assign(0.0);
  638                                                               Ref<SCElementOp> hcore_op = new OneBodyIntOp(integral()->hcore());
  639                                                               hcore.element_op(hcore_op);
  640                                                             
  641                                                               dens->scale_diagonal(0.5);
  642                                                               SCElementScalarProduct *prod = new SCElementScalarProduct;
  643                                                               prod->reference();
  644                                                               Ref<SCElementOp2> op = prod;
  645                                                               hcore->element_op(prod, dens);
  646                                                               double e = prod->result();
  647                                                               op = 0;
  648                                                               prod->dereference();
  649                                                               delete prod;
  650                                                               return 2.0 * e;
  651                                                             }
  652                                                             
  653                                                             void
  654                                                             SCF::two_body_energy(double &ec, double &ex)
  655                                                             {
  656                                                               ExEnv::err() << class_name() << ": two_body_energy not implemented" << endl;
  657                                                             }
  658                                                             
  659                                                             /////////////////////////////////////////////////////////////////////////////
  660                                                             
  661                                                             void
  662                                                             SCF::init_threads()
  663     1     0     0     0     0     0     0     0     0     0 {
  664                                                               int nthread = threadgrp_->nthread();
  665                                                               int int_store = integral()->storage_unused()/nthread;
  666                                                               
  667                                                               // initialize the two electron integral classes
  668                                                               tbis_ = new Ref<TwoBodyInt>[nthread];
  669                                                               for (int i=0; i < nthread; i++) {
  670                                                                 tbis_[i] = integral()->electron_repulsion();
  671                                                                 tbis_[i]->set_integral_storage(int_store);
  672                                                               }
  673                                                             
  674                                                             }
  675                                                             
  676                                                             void
  677                                                             SCF::done_threads()
  678                                                             {
  679                                                               for (int i=0; i < threadgrp_->nthread(); i++) tbis_[i] = 0;
  680                                                               delete[] tbis_;
  681                                                               tbis_ = 0;
  682                                                             }
  683                                                             
  684                                                             int *
  685                                                             SCF::read_occ(const Ref<KeyVal> &keyval, const char *name, int nirrep)
  686                                                             {
  687                                                               int *occ = 0;
  688                                                               if (keyval->exists(name)) {
  689                                                                 if (keyval->count(name) != nirrep) {
  690                                                                   ExEnv::err() << node0 << indent
  691                                                                                << "ERROR: SCF: have " << nirrep << " irreps but "
  692                                                                                << name << " vector is length " << keyval->count(name)
  693                                                                                << endl;
  694                                                                   abort();
  695                                                                 }
  696                                                                 occ = new int[nirrep];
  697                                                                 for (int i=0; i<nirrep; i++) {
  698                                                                   occ[i] = keyval->intvalue(name,i);
  699                                                                 }
  700                                                               }
  701                                                               return occ;
  702                                                             }
  703                                                             
  704                                                             /////////////////////////////////////////////////////////////////////////////
  705                                                             
  706                                                             // Local Variables:
  707                                                             // mode: c++
  708                                                             // c-file-style: "ETS"
  709                                                             // End:
  710