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 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:
Generated at Fri Jan 10 08:14:09 2003 for MPQC 2.1.3 using the documentation package Doxygen 1.2.14.