aotoso.cc

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

            // find the linearly independent SO's for this irrep/component
            memcpy(linorbcop[0], linorb[0], nfuncuniq*nfuncall*sizeof(double));
            double *singval = new double[nfuncuniq];
            double djunk=0; int ijunk=1;
            int lwork = 5*nfuncall;
            double *work = new double[lwork];
            int info;
            // solves At = V SIGMA Ut (since FORTRAN array layout is used)
            F77_DGESVD("N","A",&nfuncall,&nfuncuniq,linorb[0],&nfuncall,
                       singval, &djunk, &ijunk, u[0], &nfuncuniq,
                       work, &lwork, &info);
            // put the lin indep symm orbs into the so array
            for (j=0; j<nfuncuniq; j++) {
              if (singval[j] > 1.0e-6) {
                SO tso;
                tso.set_length(nfuncall);
                int ll = 0, llnonzero = 0;
                for (int k=0; k<nequiv; k++) {
                  for (int l=0; l<nfuncuniq; l++,ll++) {
                    double tmp = 0.0;
                    for (int m=0; m<nfuncuniq; m++) {
                      tmp += u[m][j] * linorbcop[m][ll] / singval[j];
                    }
                    if (fabs(tmp) > DBL_EPSILON) {
                      int equivatom = gbs.molecule()->equivalent(iuniq,k);
                      tso.cont[llnonzero].bfn
                        = l
                        + bfn_offset_in_shell
                        + gbs.shell_to_function(gbs.shell_on_center(equivatom,
                                                                    s));
                      tso.cont[llnonzero].coef = tmp;
                      llnonzero++;
                    }
                  }
                }
                tso.reset_length(llnonzero);
                if (llnonzero == 0) {
                  ExEnv::errn() << "aotoso: internal error: no bfns in SO"
                               << endl;
                  abort();
                }
                if (SOs[irnum+comp].add(tso,saoelem[irnum+comp])) {
                  saoelem[irnum+comp]++;
                }
                else {
                  ExEnv::errn() << "aotoso: internal error: "
                               << "impossible duplicate SO"
                               << endl;
                  abort();
                }
              }
            }
            delete[] singval;
            delete[] work;
          }
          irnum += ct.gamma(ir).degeneracy();
        }
        bfn_offset_in_shell += gbs(i,s).nfunction(c);

        delete[] linorb[0];
        delete[] linorb;
        delete[] linorbcop[0];
        delete[] linorbcop;
        delete[] u[0];
        delete[] u;
      }
    }
  }

  // Make sure all the nodes agree on what the symmetry orbitals are.
  // (All the above work for me > 0 is ignored.)
  Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
  for (i=0; i<ncomp; i++) {
    int len = SOs[i].len;
    grp->bcast(len);
    SOs[i].reset_length(len);
    for (j=0; j<len; j++) {
      int solen = SOs[i].so[j].length;
      grp->bcast(solen);
      SOs[i].so[j].reset_length(solen);
      for (int k=0; k<solen; k++) {
        grp->bcast(SOs[i].so[j].cont[k].bfn);
        grp->bcast(SOs[i].so[j].cont[k].coef);
      }
    }
  }

  for (i=0; i < ncomp; i++) {
    ir = whichir[i];
    int scal = ct.gamma(ir).complex() ? 2 : 1;

    if (saoelem[i] < nbf_in_ir_[ir]/scal) {
      // if we found too few, there are big problems
      
      ExEnv::err0() << indent
           << scprintf("trouble making SO's for irrep %s\n",
                       ct.gamma(ir).symbol());
      ExEnv::err0() << indent
           << scprintf("  only found %d out of %d SO's\n",
                       saoelem[i], nbf_in_ir_[ir]/scal);
      SOs[i].print("");
      abort();

    } else if (saoelem[i] > nbf_in_ir_[ir]/scal) {
      // there are some redundant so's left...need to do something to get
      // the elements we want
      
      ExEnv::err0() << indent
           << scprintf("trouble making SO's for irrep %s\n",
                       ct.gamma(ir).symbol());
      ExEnv::err0() << indent
           << scprintf("  found %d SO's, but there should only be %d\n",
                       saoelem[i], nbf_in_ir_[ir]/scal);
      SOs[i].print("");
      abort();
    }
  }

  if (ct.complex()) {
    SO_block *nSOs = new SO_block[nblocks_];

    int in=0;

    for (i=ir=0; ir < nirrep_; ir++) {
      if (ct.gamma(ir).complex()) {
        nSOs[in].set_length(nbf_in_ir_[ir]);
        int k;
        for (k=0; k < SOs[i].len; k++)
          nSOs[in].add(SOs[i].so[k],k);
        i++;

        for (j=0; j < SOs[i].len; j++,k++)
          nSOs[in].add(SOs[i].so[j],k);

        i++;
        in++;
      } else {
        for (j=0; j < ct.gamma(ir).degeneracy(); j++,i++,in++) {
          nSOs[in].set_length(nbf_in_ir_[ir]);
          for (int k=0; k < SOs[i].len; k++)
            nSOs[in].add(SOs[i].so[k],k);
        }
      }
    }

    SO_block *tmp= SOs;
    SOs = nSOs;
    delete[] tmp;
  }

  delete[] saoelem;
  delete[] whichir;
  delete[] whichcmp;

  return SOs;
}

RefSCMatrix
PetiteList::aotoso()
{
  RefSCMatrix aoso(AO_basisdim(), SO_basisdim(), gbs_->so_matrixkit());
  aoso.assign(0.0);

  if (c1_) {
    aoso->unit();
    return aoso;
  }
  
  SO_block *sos = aotoso_info();
  
  BlockedSCMatrix *aosop = dynamic_cast<BlockedSCMatrix*>(aoso.pointer());

  for (int b=0; b < aosop->nblocks(); b++) {
    RefSCMatrix aosb = aosop->block(b);

    if (aosb.null())
      continue;
    
    SO_block& sob = sos[b];
    
    Ref<SCMatrixSubblockIter> iter =
      aosb->local_blocks(SCMatrixSubblockIter::Write);

    for (iter->begin(); iter->ready(); iter->next()) {
      if (dynamic_cast<SCMatrixRectBlock*>(iter->block())) {
        SCMatrixRectBlock *blk = dynamic_cast<SCMatrixRectBlock*>(iter->block());

        int jlen = blk->jend-blk->jstart;
    
        for (int j=0; j < sob.len; j++) {
          if (j < blk->jstart || j >= blk->jend)
            continue;
      
          SO& soj = sob.so[j];
      
          for (int i=0; i < soj.len; i++) {
            int ii=soj.cont[i].bfn;
            
            if (ii < blk->istart || ii >= blk->iend)
              continue;

            blk->data[(ii-blk->istart)*jlen+(j-blk->jstart)] =
              soj.cont[i].coef;
          }
        }
      } else {
        SCMatrixRectSubBlock *blk =
          dynamic_cast<SCMatrixRectSubBlock*>(iter->block());

        for (int j=0; j < sob.len; j++) {
          if (j < blk->jstart || j >= blk->jend)
            continue;
      
          SO& soj = sob.so[j];
      
          for (int i=0; i < soj.len; i++) {
            int ii=soj.cont[i].bfn;
        
            if (ii < blk->istart || ii >= blk->iend)
              continue;

            blk->data[ii*blk->istride+j] = soj.cont[i].coef;
          }
        }
      }
    }
  }
  
  delete[] sos;

  return aoso;
}

RefSCMatrix
PetiteList::sotoao()
{
  if (c1_)
    return aotoso();
  else if (nirrep_ == ng_) // subgroup of d2h
    return aotoso().t();
  else
    return aotoso().i();
}

RefSymmSCMatrix
PetiteList::to_SO_basis(const RefSymmSCMatrix& a)
{
  // SO basis is always blocked, so first make sure a is blocked
  RefSymmSCMatrix aomatrix = dynamic_cast<BlockedSymmSCMatrix*>(a.pointer());
  if (aomatrix.null()) {
    aomatrix = gbs_->so_matrixkit()->symmmatrix(AO_basisdim());
    aomatrix->convert(a);
  }

  // if C1, then do nothing
  if (c1_)
    return aomatrix.copy();
  
  RefSymmSCMatrix somatrix(SO_basisdim(), gbs_->so_matrixkit());
  somatrix.assign(0.0);
  somatrix->accumulate_transform(aotoso(), aomatrix,
                                 SCMatrix::TransposeTransform);

  return somatrix;
}

RefSymmSCMatrix
PetiteList::to_AO_basis(const RefSymmSCMatrix& somatrix)
{
  // if C1, then do nothing
  if (c1_)
    return somatrix.copy();
  
  RefSymmSCMatrix aomatrix(AO_basisdim(), gbs_->so_matrixkit());
  aomatrix.assign(0.0);

  if (nirrep_ == ng_) // subgroup of d2h
    aomatrix->accumulate_transform(aotoso(), somatrix);
  else
    aomatrix->accumulate_transform(aotoso().i(), somatrix,
                                   SCMatrix::TransposeTransform);

  return aomatrix;
}

RefSCMatrix
PetiteList::evecs_to_SO_basis(const RefSCMatrix& aoev)
{
  ExEnv::err0() << indent
       << "PetiteList::evecs_to_SO_basis: don't work yet\n";
  abort();
  
  RefSCMatrix aoevecs = dynamic_cast<BlockedSCMatrix*>(aoev.pointer());
  if (aoevecs.null()) {
    aoevecs = gbs_->so_matrixkit()->matrix(AO_basisdim(), AO_basisdim());
    aoevecs->convert(aoev);
  }

  RefSCMatrix soev =  aotoso().t() * aoevecs;
  soev.print("soev");

  RefSCMatrix soevecs(SO_basisdim(), SO_basisdim(), gbs_->so_matrixkit());
  soevecs->convert(soev);

  return soevecs;
}

RefSCMatrix
PetiteList::evecs_to_AO_basis(const RefSCMatrix& soevecs)
{
  // if C1, then do nothing
  if (c1_)
    return soevecs.copy();
  
  RefSCMatrix aoev = aotoso() * soevecs;

  return aoev;
}

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

void
PetiteList::symmetrize(const RefSymmSCMatrix& skel,
                       const RefSymmSCMatrix& sym)
{
  GaussianBasisSet& gbs = *gbs_.pointer();

  // SO basis is always blocked, so first make sure skel is blocked
  RefSymmSCMatrix bskel = dynamic_cast<BlockedSymmSCMatrix*>(skel.pointer());
  if (bskel.null()) {
    bskel = gbs.so_matrixkit()->symmmatrix(AO_basisdim());
    bskel->convert(skel);
  }
  
  // if C1, then do nothing
  if (c1_) {
    sym.assign(bskel);
    return;
  }
  
  int b,c;

  CharacterTable ct = gbs.molecule()->point_group()->char_table();

  RefSCMatrix aoso = aotoso();
  BlockedSCMatrix *lu = dynamic_cast<BlockedSCMatrix*>(aoso.pointer());

  for (b=0; b < lu->nblocks(); b++) {
    if (lu->block(b).null())
      continue;
    
    int ir = ct.which_irrep(b);
  
    double skal = (double)ct.order()/(double)ct.gamma(ir).degeneracy();
    skal = sqrt(skal);
    lu->block(b).scale(skal);
  }

  sym.assign(0.0);
  sym.accumulate_transform(aoso,bskel,SCMatrix::TransposeTransform);
  aoso=0;

  BlockedSymmSCMatrix *la = dynamic_cast<BlockedSymmSCMatrix*>(sym.pointer());
  
  // loop through blocks and finish symmetrizing degenerate blocks
  for (b=0; b < la->nblocks(); b++) {
    if (la->block(b).null())
      continue;

    int ir=ct.which_irrep(b);

    if (ct.gamma(ir).degeneracy()==1)
      continue;

    if (ct.gamma(ir).complex()) {
      int nbf = nbf_in_ir_[ir]/2;
      
      RefSymmSCMatrix irrep = la->block(b).get_subblock(0, nbf-1);
      irrep.accumulate(la->block(b).get_subblock(nbf, 2*nbf-1));

      RefSCMatrix sub = la->block(b).get_subblock(nbf, 2*nbf-1, 0, nbf-1);
      RefSCMatrix subt = sub.t();
      subt.scale(-1.0);
      sub.accumulate(subt);
      subt=0;

      la->block(b).assign_subblock(irrep,  0, nbf-1);
      la->block(b).assign_subblock(irrep,nbf, 2*nbf-1);
      la->block(b).assign_subblock(sub, nbf, 2*nbf-1, 0, nbf-1);

    } else {
      RefSymmSCMatrix irrep = la->block(b).copy();
      for (c=1; c < ct.gamma(ir).degeneracy(); c++)
        irrep.accumulate(la->block(b+c));
      
      for (c=0; c < ct.gamma(ir).degeneracy(); c++)
        la->block(b+c).assign(irrep);

      b += ct.gamma(ir).degeneracy()-1;
    }
  }
}

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

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