1 //
2 // uhf.cc --- implementation of the unrestricted Hartree-Fock class
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 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31 
32 #include <math.h>
33 
34 #include <util/misc/timer.h>
35 #include <util/misc/formio.h>
36 #include <util/state/stateio.h>
37 
38 #include <chemistry/qc/basis/petite.h>
39 
40 #include <chemistry/qc/scf/uhf.h>
41 #include <chemistry/qc/scf/lgbuild.h>
42 #include <chemistry/qc/scf/uhftmpl.h>
43 
44 using namespace std;
45 using namespace sc;
46 
47 ///////////////////////////////////////////////////////////////////////////
48 // UHF
49 
50 static ClassDesc UHF_cd(
51   typeid(UHF),"UHF",1,"public UnrestrictedSCF",
52   0, create<UHF>, create<UHF>);
53 
UHF(StateIn & s)54 UHF::UHF(StateIn& s) :
55   SavableState(s),
56   UnrestrictedSCF(s)
57 {
58 }
59 
UHF(const Ref<KeyVal> & keyval)60 UHF::UHF(const Ref<KeyVal>& keyval) :
61   UnrestrictedSCF(keyval)
62 {
63 }
64 
~UHF()65 UHF::~UHF()
66 {
67 }
68 
69 void
save_data_state(StateOut & s)70 UHF::save_data_state(StateOut& s)
71 {
72   UnrestrictedSCF::save_data_state(s);
73 }
74 
75 int
value_implemented() const76 UHF::value_implemented() const
77 {
78   return 1;
79 }
80 
81 int
gradient_implemented() const82 UHF::gradient_implemented() const
83 {
84   return 1;
85 }
86 
87 void
print(ostream & o) const88 UHF::print(ostream&o) const
89 {
90   UnrestrictedSCF::print(o);
91 }
92 
93 //////////////////////////////////////////////////////////////////////////////
94 
95 void
two_body_energy(double & ec,double & ex)96 UHF::two_body_energy(double &ec, double &ex)
97 {
98   tim_enter("uhf e2");
99   ec = 0.0;
100   ex = 0.0;
101   if (local_ || local_dens_) {
102     // grab the data pointers from the G and P matrices
103     double *apmat;
104     double *bpmat;
105     tim_enter("local data");
106     RefSymmSCMatrix adens = alpha_ao_density();
107     RefSymmSCMatrix bdens = beta_ao_density();
108     adens->scale(2.0);
109     adens->scale_diagonal(0.5);
110     bdens->scale(2.0);
111     bdens->scale_diagonal(0.5);
112     RefSymmSCMatrix aptmp = get_local_data(adens, apmat, SCF::Read);
113     RefSymmSCMatrix bptmp = get_local_data(bdens, bpmat, SCF::Read);
114     tim_exit("local data");
115 
116     // initialize the two electron integral classes
117     Ref<TwoBodyInt> tbi = integral()->electron_repulsion();
118     tbi->set_integral_storage(0);
119 
120     signed char * pmax = init_pmax(apmat);
121 
122     LocalUHFEnergyContribution lclc(apmat, bpmat);
123     Ref<PetiteList> pl = integral()->petite_list();
124     LocalGBuild<LocalUHFEnergyContribution>
125       gb(lclc, tbi, pl, basis(), scf_grp_, pmax,
126          desired_value_accuracy()/100.0);
127     gb.run();
128 
129     delete[] pmax;
130 
131     ec = lclc.ec;
132     ex = lclc.ex;
133   }
134   else {
135     ExEnv::err0() << indent << "Cannot yet use anything but Local matrices\n";
136     abort();
137   }
138   tim_exit("uhf e2");
139 }
140 
141 //////////////////////////////////////////////////////////////////////////////
142 
143 void
ao_fock(double accuracy)144 UHF::ao_fock(double accuracy)
145 {
146   Ref<PetiteList> pl = integral()->petite_list(basis());
147 
148   // calculate G.  First transform diff_densa_ to the AO basis, then
149   // scale the off-diagonal elements by 2.0
150   RefSymmSCMatrix dda = diff_densa_;
151   diff_densa_ = pl->to_AO_basis(dda);
152   diff_densa_->scale(2.0);
153   diff_densa_->scale_diagonal(0.5);
154 
155   RefSymmSCMatrix ddb = diff_densb_;
156   diff_densb_ = pl->to_AO_basis(ddb);
157   diff_densb_->scale(2.0);
158   diff_densb_->scale_diagonal(0.5);
159 
160   // now try to figure out the matrix specialization we're dealing with
161   // if we're using Local matrices, then there's just one subblock, or
162   // see if we can convert G and P to local matrices
163   if (local_ || local_dens_) {
164     double *gmat, *gmato, *pmat, *pmato;
165 
166     // grab the data pointers from the G and P matrices
167     RefSymmSCMatrix gtmp = get_local_data(gmata_, gmat, SCF::Accum);
168     RefSymmSCMatrix ptmp = get_local_data(diff_densa_, pmat, SCF::Read);
169     RefSymmSCMatrix gotmp = get_local_data(gmatb_, gmato, SCF::Accum);
170     RefSymmSCMatrix potmp = get_local_data(diff_densb_, pmato, SCF::Read);
171 
172     signed char * pmax = init_pmax(pmat);
173 
174 //      LocalUHFContribution lclc(gmat, pmat, gmato, pmato);
175 //      LocalGBuild<LocalUHFContribution>
176 //        gb(lclc, tbi_, pl, basis(), scf_grp_, pmax,
177 //           desired_value_accuracy()/100.0);
178 //      gb.run();
179     int i;
180     int nthread = threadgrp_->nthread();
181     LocalGBuild<LocalUHFContribution> **gblds =
182       new LocalGBuild<LocalUHFContribution>*[nthread];
183     LocalUHFContribution **conts = new LocalUHFContribution*[nthread];
184 
185     double **gmats = new double*[nthread];
186     gmats[0] = gmat;
187     double **gmatos = new double*[nthread];
188     gmatos[0] = gmato;
189 
190     Ref<GaussianBasisSet> bs = basis();
191     int ntri = i_offset(bs->nbasis());
192 
193     double gmat_accuracy = accuracy;
194     if (min_orthog_res() < 1.0) { gmat_accuracy *= min_orthog_res(); }
195 
196     for (i=0; i < nthread; i++) {
197       if (i) {
198         gmats[i] = new double[ntri];
199         memset(gmats[i], 0, sizeof(double)*ntri);
200         gmatos[i] = new double[ntri];
201         memset(gmatos[i], 0, sizeof(double)*ntri);
202       }
203       conts[i] = new LocalUHFContribution(gmats[i], pmat, gmatos[i], pmato);
204       gblds[i] = new LocalGBuild<LocalUHFContribution>(*conts[i], tbis_[i],
205         pl, bs, scf_grp_, pmax, gmat_accuracy, nthread, i
206         );
207 
208       threadgrp_->add_thread(i, gblds[i]);
209     }
210 
211     tim_enter("start thread");
212     if (threadgrp_->start_threads() < 0) {
213       ExEnv::err0() << indent
214            << "UHF: error starting threads" << endl;
215       abort();
216     }
217     tim_exit("start thread");
218 
219     tim_enter("stop thread");
220     if (threadgrp_->wait_threads() < 0) {
221       ExEnv::err0() << indent
222            << "UHF: error waiting for threads" << endl;
223       abort();
224     }
225     tim_exit("stop thread");
226 
227     double tnint=0;
228     for (i=0; i < nthread; i++) {
229       tnint += gblds[i]->tnint;
230 
231       if (i) {
232         for (int j=0; j < ntri; j++) {
233           gmat[j] += gmats[i][j];
234           gmato[j] += gmatos[i][j];
235         }
236         delete[] gmats[i];
237         delete[] gmatos[i];
238       }
239 
240       delete gblds[i];
241       delete conts[i];
242     }
243 
244     delete[] gmats;
245     delete[] gmatos;
246     delete[] gblds;
247     delete[] conts;
248 
249     delete[] pmax;
250 
251     scf_grp_->sum(&tnint, 1, 0, 0);
252     ExEnv::out0() << indent << scprintf("%20.0f integrals\n", tnint);
253 
254     // if we're running on multiple processors, then sum the G matrices
255     if (scf_grp_->n() > 1) {
256       scf_grp_->sum(gmat, i_offset(basis()->nbasis()));
257       scf_grp_->sum(gmato, i_offset(basis()->nbasis()));
258     }
259 
260     // if we're running on multiple processors, or we don't have local
261     // matrices, then accumulate gtmp back into G
262     if (!local_ || scf_grp_->n() > 1) {
263       gmata_->convert_accumulate(gtmp);
264       gmatb_->convert_accumulate(gotmp);
265     }
266   }
267 
268   // for now quit
269   else {
270     ExEnv::err0() << indent << "Cannot yet use anything but Local matrices\n";
271     abort();
272   }
273 
274   // get rid of AO delta P
275   diff_densa_ = dda;
276   dda = diff_densa_.clone();
277 
278   diff_densb_ = ddb;
279   ddb = diff_densb_.clone();
280 
281   // now symmetrize the skeleton G matrix, placing the result in dda
282   RefSymmSCMatrix skel_gmat = gmata_.copy();
283   skel_gmat.scale(1.0/(double)pl->order());
284   pl->symmetrize(skel_gmat,dda);
285 
286   skel_gmat = gmatb_.copy();
287   skel_gmat.scale(1.0/(double)pl->order());
288   pl->symmetrize(skel_gmat,ddb);
289 
290   // Fa = H+Ga
291   focka_.result_noupdate().assign(hcore_);
292   focka_.result_noupdate().accumulate(dda);
293 
294   // Fb = H+Gb
295   fockb_.result_noupdate().assign(hcore_);
296   fockb_.result_noupdate().accumulate(ddb);
297 
298   dda.assign(0.0);
299   accumddh_->accum(dda);
300   focka_.result_noupdate().accumulate(dda);
301   fockb_.result_noupdate().accumulate(dda);
302 
303   focka_.computed()=1;
304   fockb_.computed()=1;
305 }
306 
307 /////////////////////////////////////////////////////////////////////////////
308 
309 void
two_body_deriv(double * tbgrad)310 UHF::two_body_deriv(double * tbgrad)
311 {
312   two_body_deriv_hf(tbgrad, 1.0);
313 }
314 
315 /////////////////////////////////////////////////////////////////////////////
316 
317 // Local Variables:
318 // mode: c++
319 // c-file-style: "ETS"
320 // End:
321