1 // 2 // ltbgrad.h --- definition of the local two-electron gradient builder 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 #ifndef _chemistry_qc_scf_ltbgrad_h 29 #define _chemistry_qc_scf_ltbgrad_h 30 31 #ifdef __GNUC__ 32 #pragma interface 33 #endif 34 35 #include <math.h> 36 37 #include <util/misc/timer.h> 38 #include <math/scmat/offset.h> 39 40 #include <chemistry/qc/basis/tbint.h> 41 #include <chemistry/qc/basis/petite.h> 42 43 #include <chemistry/qc/scf/tbgrad.h> 44 45 namespace sc { 46 47 template<class T> 48 class LocalTBGrad : public TBGrad<T> { 49 public: 50 double *tbgrad; 51 52 protected: 53 MessageGrp *grp_; 54 TwoBodyDerivInt *tbi_; 55 GaussianBasisSet *gbs_; 56 PetiteList *rpl_; 57 Molecule *mol_; 58 59 double pmax_; 60 double accuracy_; 61 62 int threadno_; 63 int nthread_; 64 65 public: 66 LocalTBGrad(T& t, const Ref<TwoBodyDerivInt>& tbdi, const Ref<PetiteList>& pl, 67 const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g, 68 double *tbg, double pm, double a, int nt = 1, int tn = 0, 69 double exchange_fraction = 1.0) : 70 TBGrad<T>(t,exchange_fraction), 71 tbgrad(tbg), pmax_(pm), accuracy_(a), threadno_(tn), nthread_(nt) 72 { 73 grp_ = g.pointer(); 74 gbs_ = bs.pointer(); 75 rpl_ = pl.pointer(); 76 tbi_ = tbdi.pointer(); 77 mol_ = gbs_->molecule().pointer(); 78 } 79 ~LocalTBGrad()80 ~LocalTBGrad() {} 81 run()82 void run() { 83 int me = grp_->me(); 84 int nproc = grp_->n(); 85 86 // grab ref for convenience 87 GaussianBasisSet& gbs = *gbs_; 88 Molecule& mol = *mol_; 89 PetiteList& pl = *rpl_; 90 TwoBodyDerivInt& tbi = *tbi_; 91 92 // create vector to hold skeleton gradient 93 double *tbint = new double[mol.natom()*3]; 94 memset(tbint, 0, sizeof(double)*mol.natom()*3); 95 96 // for bounds checking 97 int PPmax = (int) (log(6.0*pmax_*pmax_)/log(2.0)); 98 int threshold = (int) (log(accuracy_)/log(2.0)); 99 100 int kindex=0; 101 int threadind=0; 102 for (int i=0; i < gbs.nshell(); i++) { 103 if (!pl.in_p1(i)) 104 continue; 105 106 int ni=gbs(i).nfunction(); 107 int fi=gbs.shell_to_function(i); 108 109 for (int j=0; j <= i; j++) { 110 int ij=i_offset(i)+j; 111 if (!pl.in_p2(ij)) 112 continue; 113 114 if (tbi.log2_shell_bound(i,j,-1,-1)+PPmax < threshold) 115 continue; 116 117 int nj=gbs(j).nfunction(); 118 int fj=gbs.shell_to_function(j); 119 120 for (int k=0; k <= i; k++,kindex++) { 121 if (kindex%nproc != me) 122 continue; 123 124 threadind++; 125 if (threadind % nthread_ != threadno_) 126 continue; 127 128 int nk=gbs(k).nfunction(); 129 int fk=gbs.shell_to_function(k); 130 131 for (int l=0; l <= ((i==k)?j:k); l++) { 132 if (tbi.log2_shell_bound(i,j,k,l)+PPmax < threshold) 133 continue; 134 135 int kl=i_offset(k)+l; 136 int qijkl; 137 if (!(qijkl=pl.in_p4(ij,kl,i,j,k,l))) 138 continue; 139 140 int nl=gbs(l).nfunction(); 141 int fl=gbs.shell_to_function(l); 142 143 DerivCenters cent; 144 tbi.compute_shell(i,j,k,l,cent); 145 146 const double * buf = tbi.buffer(); 147 148 double cscl, escl; 149 150 this->set_scale(cscl, escl, i, j, k, l); 151 152 int indijkl=0; 153 int nx=cent.n(); 154 //if (cent.has_omitted_center()) nx--; 155 for (int x=0; x < nx; x++) { 156 int ix=cent.atom(x); 157 int io=cent.omitted_atom(); 158 for (int ixyz=0; ixyz < 3; ixyz++) { 159 double tx = tbint[ixyz+ix*3]; 160 double to = tbint[ixyz+io*3]; 161 162 for (int ip=0, ii=fi; ip < ni; ip++, ii++) { 163 for (int jp=0, jj=fj; jp < nj; jp++, jj++) { 164 for (int kp=0, kk=fk; kp < nk; kp++, kk++) { 165 for (int lp=0, ll=fl; lp < nl; lp++, ll++, indijkl++) { 166 double contrib; 167 double qint = buf[indijkl]*qijkl; 168 169 contrib = cscl*qint* 170 TBGrad<T>::contribution.cont1(ij_offset(ii,jj), 171 ij_offset(kk,ll)); 172 173 tx += contrib; 174 to -= contrib; 175 176 contrib = escl*qint* 177 TBGrad<T>::contribution.cont2(ij_offset(ii,kk), 178 ij_offset(jj,ll)); 179 180 tx += contrib; 181 to -= contrib; 182 183 if (i!=j && k!=l) { 184 contrib = escl*qint* 185 TBGrad<T>::contribution.cont2(ij_offset(ii,ll), 186 ij_offset(jj,kk)); 187 188 tx += contrib; 189 to -= contrib; 190 } 191 } 192 } 193 } 194 } 195 196 tbint[ixyz+ix*3] = tx; 197 tbint[ixyz+io*3] = to; 198 } 199 } 200 } 201 } 202 } 203 } 204 205 CharacterTable ct = mol.point_group()->char_table(); 206 SymmetryOperation so; 207 208 for (int alpha=0; alpha < mol.natom(); alpha++) { 209 double tbx = tbint[alpha*3+0]; 210 double tby = tbint[alpha*3+1]; 211 double tbz = tbint[alpha*3+2]; 212 213 for (int g=1; g < ct.order(); g++) { 214 so = ct.symm_operation(g); 215 int ap = pl.atom_map(alpha,g); 216 217 tbx += tbint[ap*3+0]*so(0,0) + tbint[ap*3+1]*so(1,0) + 218 tbint[ap*3+2]*so(2,0); 219 tby += tbint[ap*3+0]*so(0,1) + tbint[ap*3+1]*so(1,1) + 220 tbint[ap*3+2]*so(2,1); 221 tbz += tbint[ap*3+0]*so(0,2) + tbint[ap*3+1]*so(1,2) + 222 tbint[ap*3+2]*so(2,2); 223 } 224 double scl = 1.0/(double)ct.order(); 225 tbgrad[alpha*3+0] += tbx*scl; 226 tbgrad[alpha*3+1] += tby*scl; 227 tbgrad[alpha*3+2] += tbz*scl; 228 } 229 230 delete[] tbint; 231 } 232 }; 233 234 } 235 236 #endif 237 238 // Local Variables: 239 // mode: c++ 240 // c-file-style: "ETS" 241 // End: 242