[0b990d] | 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 |
|
---|
| 80 | ~LocalTBGrad() {}
|
---|
| 81 |
|
---|
| 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:
|
---|