# ltbgrad.h

    00001 //
    00002 // ltbgrad.h --- definition of the local two-electron gradient builder
    00003 //
    00004 // Copyright (C) 1996 Limit Point Systems, Inc.
    00005 //
    00006 // Author: Edward Seidl <seidl@janed.com>
    00007 // Maintainer: LPS
    00008 //
    00009 // This file is part of the SC Toolkit.
    00010 //
    00011 // The SC Toolkit is free software; you can redistribute it and/or modify
    00012 // it under the terms of the GNU Library General Public License as published by
    00013 // the Free Software Foundation; either version 2, or (at your option)
    00014 // any later version.
    00015 //
    00016 // The SC Toolkit is distributed in the hope that it will be useful,
    00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
    00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    00019 // GNU Library General Public License for more details.
    00020 //
    00021 // You should have received a copy of the GNU Library General Public License
    00022 // along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
    00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
    00024 //
    00025 // The U.S. Government is granted a limited license as per AL 91-7.
    00026 //
    00027
    00028 #ifndef _chemistry_qc_scf_ltbgrad_h
    00029 #define _chemistry_qc_scf_ltbgrad_h
    00030
    00031 #ifdef __GNUC__
    00032 #pragma interface
    00033 #endif
    00034
    00035 #include <math.h>
    00036
    00037 #include <util/misc/timer.h>
    00038 #include <math/scmat/offset.h>
    00039
    00040 #include <chemistry/qc/basis/tbint.h>
    00041 #include <chemistry/qc/basis/petite.h>
    00042
    00043 #include <chemistry/qc/scf/tbgrad.h>
    00044
    00045 namespace sc {
    00046
    00047 template<class T>
    00048 class LocalTBGrad : public TBGrad<T> {
    00049   public:
    00050     double *tbgrad;
    00051
    00052   protected:
    00053     MessageGrp *grp_;
    00054     TwoBodyDerivInt *tbi_;
    00055     GaussianBasisSet *gbs_;
    00056     PetiteList *rpl_;
    00057     Molecule *mol_;
    00058
    00059     double pmax_;
    00060     double accuracy_;
    00061
    00062     int threadno_;
    00063     int nthread_;
    00064
    00065   public:
    00066     LocalTBGrad(T& t, const Ref<TwoBodyDerivInt>& tbdi, const Ref<PetiteList>& pl,
    00067                 const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g,
    00068                 double *tbg, double pm, double a, int nt = 1, int tn = 0,
    00069                 double exchange_fraction = 1.0) :
    00070       TBGrad<T>(t,exchange_fraction),
    00071       tbgrad(tbg), pmax_(pm), accuracy_(a), threadno_(tn), nthread_(nt)
    00072     {
    00073       grp_ = g.pointer();
    00074       gbs_ = bs.pointer();
    00075       rpl_ = pl.pointer();
    00076       tbi_ = tbdi.pointer();
    00077       mol_ = gbs_->molecule().pointer();
    00078     }
    00079
    00080     ~LocalTBGrad() {}
    00081
    00082     void run() {
    00083       int me = grp_->me();
    00084       int nproc = grp_->n();
    00085
    00086       // grab ref for convenience
    00087       GaussianBasisSet& gbs = *gbs_;
    00088       Molecule& mol = *mol_;
    00089       PetiteList& pl = *rpl_;
    00090       TwoBodyDerivInt& tbi = *tbi_;
    00091
    00092       // create vector to hold skeleton gradient
    00093       double *tbint = new double[mol.natom()*3];
    00094       memset(tbint, 0, sizeof(double)*mol.natom()*3);
    00095
    00096       // for bounds checking
    00097       int PPmax = (int) (log(6.0*pmax_*pmax_)/log(2.0));
    00098       int threshold = (int) (log(accuracy_)/log(2.0));
    00099
    00100       int kindex=0;
    00101       int threadind=0;
    00102       for (int i=0; i < gbs.nshell(); i++) {
    00103         if (!pl.in_p1(i))
    00104           continue;
    00105
    00106         int ni=gbs(i).nfunction();
    00107         int fi=gbs.shell_to_function(i);
    00108
    00109         for (int j=0; j <= i; j++) {
    00110           int ij=i_offset(i)+j;
    00111           if (!pl.in_p2(ij))
    00112             continue;
    00113
    00114           if (tbi.log2_shell_bound(i,j,-1,-1)+PPmax < threshold)
    00115             continue;
    00116
    00117           int nj=gbs(j).nfunction();
    00118           int fj=gbs.shell_to_function(j);
    00119
    00120           for (int k=0; k <= i; k++,kindex++) {
    00121             if (kindex%nproc != me)
    00122               continue;
    00123
    00124             threadind++;
    00125             if (threadind % nthread_ != threadno_)
    00126               continue;
    00127
    00128             int nk=gbs(k).nfunction();
    00129             int fk=gbs.shell_to_function(k);
    00130
    00131             for (int l=0; l <= ((i==k)?j:k); l++) {
    00132               if (tbi.log2_shell_bound(i,j,k,l)+PPmax < threshold)
    00133                 continue;
    00134
    00135               int kl=i_offset(k)+l;
    00136               int qijkl;
    00137               if (!(qijkl=pl.in_p4(ij,kl,i,j,k,l)))
    00138                 continue;
    00139
    00140               int nl=gbs(l).nfunction();
    00141               int fl=gbs.shell_to_function(l);
    00142
    00143               DerivCenters cent;
    00144               tbi.compute_shell(i,j,k,l,cent);
    00145
    00146               const double * buf = tbi.buffer();
    00147
    00148               double cscl, escl;
    00149
    00150               this->set_scale(cscl, escl, i, j, k, l);
    00151
    00152               int indijkl=0;
    00153               int nx=cent.n();
    00154               //if (cent.has_omitted_center()) nx--;
    00155               for (int x=0; x < nx; x++) {
    00156                 int ix=cent.atom(x);
    00157                 int io=cent.omitted_atom();
    00158                 for (int ixyz=0; ixyz < 3; ixyz++) {
    00159                   double tx = tbint[ixyz+ix*3];
    00160                   double to = tbint[ixyz+io*3];
    00161
    00162                   for (int ip=0, ii=fi; ip < ni; ip++, ii++) {
    00163                     for (int jp=0, jj=fj; jp < nj; jp++, jj++) {
    00164                       for (int kp=0, kk=fk; kp < nk; kp++, kk++) {
    00165                         for (int lp=0, ll=fl; lp < nl; lp++, ll++, indijkl++) {
    00166                           double contrib;
    00167                           double qint = buf[indijkl]*qijkl;
    00168
    00169                           contrib = cscl*qint*
    00170                             TBGrad<T>::contribution.cont1(ij_offset(ii,jj),
    00171                                                ij_offset(kk,ll));
    00172
    00173                           tx += contrib;
    00174                           to -= contrib;
    00175
    00176                           contrib = escl*qint*
    00177                             TBGrad<T>::contribution.cont2(ij_offset(ii,kk),
    00178                                                ij_offset(jj,ll));
    00179
    00180                           tx += contrib;
    00181                           to -= contrib;
    00182
    00183                           if (i!=j && k!=l) {
    00184                             contrib = escl*qint*
    00185                               TBGrad<T>::contribution.cont2(ij_offset(ii,ll),
    00186                                                  ij_offset(jj,kk));
    00187
    00188                             tx += contrib;
    00189                             to -= contrib;
    00190                           }
    00191                         }
    00192                       }
    00193                     }
    00194                   }
    00195
    00196                   tbint[ixyz+ix*3] = tx;
    00197                   tbint[ixyz+io*3] = to;
    00198                 }
    00199               }
    00200             }
    00201           }
    00202         }
    00203       }
    00204
    00205       CharacterTable ct = mol.point_group()->char_table();
    00206       SymmetryOperation so;
    00207
    00208       for (int alpha=0; alpha < mol.natom(); alpha++) {
    00209         double tbx = tbint[alpha*3+0];
    00210         double tby = tbint[alpha*3+1];
    00211         double tbz = tbint[alpha*3+2];
    00212
    00213         for (int g=1; g < ct.order(); g++) {
    00214           so = ct.symm_operation(g);
    00215           int ap = pl.atom_map(alpha,g);
    00216
    00217           tbx += tbint[ap*3+0]*so(0,0) + tbint[ap*3+1]*so(1,0) +
    00218                  tbint[ap*3+2]*so(2,0);
    00219           tby += tbint[ap*3+0]*so(0,1) + tbint[ap*3+1]*so(1,1) +
    00220                  tbint[ap*3+2]*so(2,1);
    00221           tbz += tbint[ap*3+0]*so(0,2) + tbint[ap*3+1]*so(1,2) +
    00222                  tbint[ap*3+2]*so(2,2);
    00223         }
    00224         double scl = 1.0/(double)ct.order();
    00225         tbgrad[alpha*3+0] += tbx*scl;
    00226         tbgrad[alpha*3+1] += tby*scl;
    00227         tbgrad[alpha*3+2] += tbz*scl;
    00228       }
    00229
    00230       delete[] tbint;
    00231     }
    00232 };
    00233
    00234 }
    00235
    00236 #endif
    00237
    00238 // Local Variables:
    00239 // mode: c++
    00240 // c-file-style: "ETS"
    00241 // End: