MPQC 2.3.1
ltbgrad.h
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
45namespace sc {
46
47template<class T>
48class 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();
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:
The CharacterTable class provides a workable character table for all of the non-cubic point groups.
Definition: pointgrp.h:325
int order() const
Returns the order of the point group.
Definition: pointgrp.h:372
SymmetryOperation & symm_operation(int i)
Returns the i'th symmetry operation.
Definition: pointgrp.h:378
DerivCenters keeps track the centers that derivatives are taken with respect to.
Definition: dercent.h:41
int omitted_atom() const
Definition: dercent.h:100
int atom(int i) const
Definition: dercent.h:91
int n() const
The number of centers for which derivatives have been computed.
Definition: dercent.h:83
The GaussianBasisSet class is used describe a basis set composed of atomic gaussian orbitals.
Definition: gaussbas.h:146
int nshell() const
Return the number of shells.
Definition: gaussbas.h:417
int shell_to_function(int i) const
Return the number of the first function in the given shell.
Definition: gaussbas.h:454
Ref< Molecule > molecule() const
Return the Molecule object.
Definition: gaussbas.h:406
Definition: ltbgrad.h:48
void run()
This is called with the Thread is run from a ThreadGrp.
Definition: ltbgrad.h:82
The MessageGrp abstract class provides a mechanism for moving data and objects between nodes in a par...
Definition: message.h:109
int me()
Returns my processor number. In the range [0,n()).
Definition: message.h:138
int n()
Returns the number of processors.
Definition: message.h:136
The Molecule class contains information about molecules.
Definition: molecule.h:128
Ref< PointGroup > point_group() const
Returns the PointGroup of the molecule.
int natom() const
Returns the number of atoms in the molcule.
Definition: molecule.h:288
Definition: petite.h:119
A template class that maintains references counts.
Definition: ref.h:332
T * pointer() const
Returns a pointer the reference counted object.
Definition: ref.h:383
The SymmetryOperation class provides a 3 by 3 matrix representation of a symmetry operation,...
Definition: pointgrp.h:70
Definition: tbgrad.h:40
This is an abstract base type for classes that compute integrals involving two electrons.
Definition: tbint.h:422
virtual void compute_shell(int, int, int, int, DerivCenters &)=0
Given for shell indices, this will cause the integral buffer to be filled in.
virtual int log2_shell_bound(int=-1, int=-1, int=-1, int=-1)=0
Return log base 2 of the maximum magnitude of any integral in a shell block.
const double * buffer() const
The computed shell integrals will be put in the buffer returned by this member.

Generated at Sun Nov 27 2022 00:00:00 for MPQC 2.3.1 using the documentation package Doxygen 1.9.5.