source: ThirdParty/mpqc_open/src/lib/chemistry/qc/dft/dfttest.cc

Candidate_v1.6.1
Last change on this file was 860145, checked in by Frederik Heber <heber@…>, 8 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 18.6 KB
Line 
1//
2// dfttest.cc
3//
4// Copyright (C) 1997 Limit Point Systems, Inc.
5//
6// Author: Curtis Janssen <cljanss@limitpt.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#ifndef _GNU_SOURCE
28# define _GNU_SOURCE
29#endif
30
31#include <signal.h>
32
33#include <util/misc/math.h>
34
35#include <util/misc/formio.h>
36#include <util/group/pregtime.h>
37#include <chemistry/qc/dft/functional.h>
38#include <chemistry/qc/dft/integrator.h>
39
40#include <chemistry/qc/dft/linkage.h>
41#include <chemistry/qc/scf/linkage.h>
42
43#ifdef HAVE_MPI
44#include <mpi.h>
45#include <util/group/messmpi.h>
46#endif
47
48#ifdef HAVE_FENV_H
49# include <fenv.h>
50#endif
51
52using namespace std;
53using namespace sc;
54
55inline static double
56norm(double v[3])
57{
58 double x,y,z;
59 return sqrt((x=v[0])*x + (y=v[1])*y + (z=v[2])*z);
60}
61
62inline static double
63dot(double v[3], double w[3])
64{
65 return v[0]*w[0] + v[1]*w[1] + v[2]*w[2];
66}
67
68double *
69density_matrix(const Ref<Wavefunction> &wfn)
70{
71 int nbasis = wfn->basis()->nbasis();
72 RefSymmSCMatrix adens = wfn->alpha_ao_density();
73 double * alpha_dmat = new double[(nbasis*(nbasis+1))/2];
74 adens->convert(alpha_dmat);
75 return alpha_dmat;
76}
77
78void
79get_density(PointInputData::SpinData &d, const SCVector3 &r,
80 const Ref<Wavefunction> &wfn, double *pdmat = 0)
81{
82 double *dmat;
83 if (pdmat) dmat = pdmat;
84 else dmat = density_matrix(wfn);
85
86 double * bsg_values_ = new double[3*wfn->basis()->nbasis()];
87 double * bs_values_ = new double[wfn->basis()->nbasis()];
88 GaussianBasisSet::ValueData vdat(wfn->basis(), wfn->integral());
89 wfn->basis()->grad_values(r,&vdat,bsg_values_,bs_values_);
90
91 int i, j;
92
93 double tmp = 0.0;
94 double densij;
95 double bvi, bvix, bviy, bviz;
96 double bvixx, bviyx, bviyy, bvizx, bvizy, bvizz;
97
98 const int X = PointInputData::X;
99 const int Y = PointInputData::Y;
100 const int Z = PointInputData::Z;
101 const int XX = PointInputData::XX;
102 const int YX = PointInputData::YX;
103 const int YY = PointInputData::YY;
104 const int ZX = PointInputData::ZX;
105 const int ZY = PointInputData::ZY;
106 const int ZZ = PointInputData::ZZ;
107
108 double grad[3], hess[6];
109 for (i=0; i<3; i++) grad[i] = 0.0;
110 for (i=0; i<6; i++) hess[i] = 0.0;
111
112 int nbasis_ = wfn->basis()->nbasis();
113 int ij=0;
114 for (i=0; i < nbasis_; i++) {
115 bvi = bs_values_[i];
116 bvix = bsg_values_[i*3+X];
117 bviy = bsg_values_[i*3+Y];
118 bviz = bsg_values_[i*3+Z];
119 for (j=0; j < i; j++,ij++) {
120 densij = dmat[ij];
121 double bvj = bs_values_[j];
122 double bvjx = bsg_values_[j*3+X];
123 double bvjy = bsg_values_[j*3+Y];
124 double bvjz = bsg_values_[j*3+Z];
125
126 tmp += 2.0*densij*bvi*bvj;
127
128 grad[X] += densij*(bvi*bvjx + bvj*bvix);
129 grad[Y] += densij*(bvi*bvjy + bvj*bviy);
130 grad[Z] += densij*(bvi*bvjz + bvj*bviz);
131 }
132
133 densij = dmat[ij]*bvi;
134 tmp += densij*bvi;
135 grad[X] += densij*bvix;
136 grad[Y] += densij*bviy;
137 grad[Z] += densij*bviz;
138 ij++;
139 }
140
141
142 d.rho = tmp;
143 for (i=0; i<3; i++) d.del_rho[i] = 2.0 * grad[i];
144 for (i=0; i<6; i++) d.hes_rho[i] = 2.0 * hess[i];
145
146 d.lap_rho = d.hes_rho[XX] + d.hes_rho[YY] + d.hes_rho[ZZ];
147
148 d.gamma = dot(d.del_rho,d.del_rho);
149
150 delete[] bsg_values_;
151 delete[] bs_values_;
152 if (!pdmat) delete[] dmat;
153}
154
155double
156fd_test_do_point(const SCVector3 &point,
157 const Ref<DenFunctional> &func, const Ref<Wavefunction> &wfn,
158 double *frozen_dmat = 0)
159{
160 PointInputData id(point);
161 get_density(id.a, point, wfn, frozen_dmat);
162 id.compute_derived(0,func->need_density_gradient(),false);
163 PointOutputData od;
164 if ( (id.a.rho + id.b.rho) > 1e2*DBL_EPSILON) func->point(id, od);
165 else return 0.0;
166 return od.energy;
167}
168
169void
170fd_test_point(int acenter, const SCVector3 &tpoint,
171 const Ref<DenFunctional> &functional, const Ref<Wavefunction> &wfn)
172{
173 SCVector3 point(tpoint);
174 Ref<Molecule> mol = wfn->molecule();
175
176 double *fd_grad_f = new double[mol->natom()*3];
177 memset(fd_grad_f,0, 3*mol->natom() * sizeof(double));
178
179 double *dmat = density_matrix(wfn);
180
181 // frozen_dmat = dmat makes the overlap derivative contribution 0
182 double *frozen_dmat = dmat;
183
184 double delta = 0.0001;
185 for (int i=0; i<mol->natom(); i++) {
186 for (int j=0; j<3; j++) {
187 if (acenter == i) point[j] += delta;
188 mol->r(i,j) += delta;
189 wfn->obsolete();
190 double f_plus = fd_test_do_point(point, functional, wfn, frozen_dmat);
191 if (acenter == i) point[j] -= 2*delta;
192 mol->r(i,j) -= 2*delta;
193 wfn->obsolete();
194 double f_minus = fd_test_do_point(point, functional, wfn, frozen_dmat);
195 if (acenter == i) point[j] += delta;
196 mol->r(i,j) += delta;
197 wfn->obsolete();
198 fd_grad_f[i*3+j] = (f_plus-f_minus)/(2.0*delta);
199 }
200 }
201
202 double * bsh_values_ = new double[6*wfn->basis()->nbasis()];
203 double * bsg_values_ = new double[3*wfn->basis()->nbasis()];
204 double * bs_values_ = new double[wfn->basis()->nbasis()];
205 GaussianBasisSet::ValueData vdat(wfn->basis(), wfn->integral());
206 wfn->basis()->hessian_values(point,&vdat,bsh_values_,bsg_values_,bs_values_);
207
208 PointInputData id(point);
209 get_density(id.a, point, wfn);
210 id.compute_derived(0,functional->need_density_gradient(),false);
211
212 PointOutputData od;
213 functional->set_compute_potential(1);
214 if ( (id.a.rho + id.b.rho) > 1e2*DBL_EPSILON) functional->point(id, od);
215 else od.zero();
216
217 double *an_grad_f = new double[mol->natom()*3];
218 memset(an_grad_f,0, 3*mol->natom() * sizeof(double));
219
220 int ncontrib = wfn->basis()->nshell();
221 int ncontrib_bf = wfn->basis()->nbasis();
222 int *contrib = new int[ncontrib];
223 int *contrib_bf = new int[ncontrib_bf];
224 for (int i=0; i<ncontrib; i++) contrib[i] = i;
225 for (int i=0; i<ncontrib_bf; i++) contrib_bf[i] = i;
226 functional->gradient(id, od, an_grad_f, acenter, wfn->basis(),
227 dmat, dmat,
228 ncontrib, contrib,
229 ncontrib_bf, contrib_bf,
230 bs_values_, bsg_values_, bsh_values_);
231 delete[] contrib;
232 delete[] contrib_bf;
233
234 cout << " acenter = " << acenter << " point = " << point << endl;
235 cout << "FD df/dx:" << endl;
236 for (int i=0; i<mol->natom(); i++) {
237 cout << scprintf(" % 16.12f % 16.12f % 16.12f",
238 fd_grad_f[3*i+0],
239 fd_grad_f[3*i+1],
240 fd_grad_f[3*i+2])
241 << endl;
242 }
243
244 cout << "AN df/dx:" << endl;
245 for (int i=0; i<mol->natom(); i++) {
246 cout << scprintf(" % 16.12f % 16.12f % 16.12f",
247 an_grad_f[3*i+0],
248 an_grad_f[3*i+1],
249 an_grad_f[3*i+2])
250 << endl;
251 }
252
253 delete[] bsh_values_;
254 delete[] bsg_values_;
255 delete[] bs_values_;
256 delete[] dmat;
257 delete[] fd_grad_f;
258 delete[] an_grad_f;
259}
260
261void
262fd_test(const Ref<DenFunctional> &functional, const Ref<Wavefunction> &wfn)
263{
264 cout << "fd_test with functional:" << endl;
265 cout << functional;
266 for (int i=0; i<wfn->molecule()->natom(); i++) {
267 for (double x = -1.0; x <= 1.0; x++) {
268 for (double y = -1.0; y <= 1.0; y++) {
269 for (double z = -1.0; z <= 1.0; z++) {
270 fd_test_point(i, SCVector3(x,y,z), functional, wfn);
271 }
272 }
273 }
274 }
275}
276
277void
278fd_e_test(const Ref<Wavefunction> &wfn)
279{
280 Ref<Molecule> mol = wfn->molecule();
281
282 cout << "Testing dE/dx with:" << endl;
283 cout << incindent;
284 wfn->print();
285 cout << decindent;
286
287 double *fd_grad_e = new double[mol->natom()*3];
288 memset(fd_grad_e,0, 3*mol->natom() * sizeof(double));
289
290 double delta = 0.0001;
291 for (int i=0; i<mol->natom(); i++) {
292 for (int j=0; j<3; j++) {
293 mol->r(i,j) += delta;
294 wfn->obsolete();
295 double e_plus = wfn->energy();
296 mol->r(i,j) -= 2*delta;
297 wfn->obsolete();
298 double e_minus = wfn->energy();
299 mol->r(i,j) += delta;
300 wfn->obsolete();
301 fd_grad_e[i*3+j] = (e_plus-e_minus)/(2.0*delta);
302 }
303 }
304
305 double *an_grad_e = new double[mol->natom()*3];
306 memset(an_grad_e,0, 3*mol->natom() * sizeof(double));
307
308 RefSCVector grad = wfn->get_cartesian_gradient();
309 grad->convert(an_grad_e);
310
311 cout << "FD dE/dx:" << endl;
312 for (int i=0; i<mol->natom(); i++) {
313 cout << scprintf(" % 16.12f % 16.12f % 16.12f",
314 fd_grad_e[3*i+0],
315 fd_grad_e[3*i+1],
316 fd_grad_e[3*i+2])
317 << endl;
318 }
319
320 cout << "AN dE/dx:" << endl;
321 for (int i=0; i<mol->natom(); i++) {
322 cout << scprintf(" % 16.12f % 16.12f % 16.12f",
323 an_grad_e[3*i+0],
324 an_grad_e[3*i+1],
325 an_grad_e[3*i+2])
326 << endl;
327 }
328
329 delete[] fd_grad_e;
330 delete[] an_grad_e;
331}
332
333#define USE_ORIGINAL_CODE 0
334#define USE_MPQC_CODE 1
335#if USE_ORIGINAL_CODE
336extern "C" easypbe_(double *up,double *agrup,double *delgrup,double *uplap,
337 double *dn,double *agrdn,double *delgrdn,double *dnlap,
338 double *agr,double *delgr,int *lcor,int *lpot,
339 double *exlsd,double *vxuplsd,double *vxdnlsd,
340 double *eclsd,double *vcuplsd,double *vcdnlsd,
341 double *expw91,double *vxuppw91,double *vxdnpw91,
342 double *ecpw91,double *vcuppw91,double *vcdnpw91,
343 double *expbe,double *vxuppbe,double *vxdnpbe,
344 double *ecpbe,double *vcuppbe,double *vcdnpbe);
345#endif
346
347void
348do_valtest(const Ref<DenFunctional> &valtest)
349{
350 valtest->set_spin_polarized(1);
351
352 SCVector3 point = 0.0;
353 PointInputData id(point);
354 // zero out data that is never used
355 int ii=0;
356 for (ii=0; ii<3; ii++) id.a.del_rho[ii] = 0.0;
357 for (ii=0; ii<3; ii++) id.b.del_rho[ii] = 0.0;
358 id.a.lap_rho = 0.0;
359 id.b.lap_rho = 0.0;
360 for (ii=0; ii<6; ii++) id.a.hes_rho[ii] = 0.0;
361 for (ii=0; ii<6; ii++) id.b.hes_rho[ii] = 0.0;
362
363 // taken from PBE.f (PBE alpha2.1)
364 const double thrd = 1.0/3.0;
365 const double thrd2 = 2.0/3.0;
366 const double pi = M_PI;
367 double conf=pow((3.0*pi*pi),thrd);
368 double conrs=pow((3.0/(4.0*pi)),thrd);
369 cout << " Fup Fdn Zup Zdn Exc" << endl;
370 // BEGIN THE LOOP THAT SELECTS A TRIAL DENSITY
371 // spin-densities are of the form
372 // rho(r)=f*(Z**3/pi)*dexp(-2*Z*r)
373 // delzdn=small change in zdn to test potentials
374 // jdens=counter for which density
375 for (int jdens = 1; jdens <= 10; jdens++) {
376 double fup=1.0;
377 double fdn=0.2*(jdens-1);
378 double zup=1.0;
379 double zdn=0.5;
380 if (jdens > 6) {
381 fdn=1.0;
382 zup=0.5+0.5*(jdens-7);
383 zdn=zup;
384 }
385 double delzdn=1e-5;
386 double sumexc, mpqc_sumexc;
387 double sumexco;
388 // BEGIN THE LOOP THAT INCREMENTS THE DENSITY DIFFERENTIALLY
389 // kdif=1=>density as above
390 // kdif=2=>Zdn changed by DELZDN
391 for (int kdif=1; kdif<=2; kdif++) {
392 if (kdif == 2) zdn=zdn+delzdn;
393 // BEGIN THE RADIAL LOOP
394 // sumexc=integrated exchange-correlation energy
395 // chng1=integrated xc energy change, based on vxc
396 // nr=number of points in radial loop
397 // rf=final value of r in integrals
398 // dr=change in r
399 // wt=weight of r in trapezoidal rule
400 // dup=up density
401 // agrup=|grad up|
402 // delgrup=(grad up).(grad |grad up|)
403 // uplap=grad^2 up=Laplacian of up
404 // dn,agrdn,delgrdn,dnlap=corresponding down quantities
405 // d=up+dn
406 // agrad=|grad rho|
407 // delgrad=(grad rho).(grad |grad rho|)
408 sumexc=0.0;
409 mpqc_sumexc = 0.0;
410 sumexco=0.0;
411 double chng1=0.0;
412 int nr=10000;
413 double rf=20.0;
414 double dr=rf/nr;
415 for (int i=1; i<=nr; i++) {
416 double r=i*dr;
417 double wt=4.*pi*r*r*dr;
418 double dup=fup*(zup*zup*zup/pi)*exp(-2.0*zup*r);
419 double ddn=fdn*(zdn*zdn*zdn/pi)*exp(-2.0*zdn*r);
420 if (dup+ddn < DBL_EPSILON) continue;
421 double zdnnu=zdn+delzdn;
422 double delddn=fdn*(zdnnu*zdnnu*zdnnu/pi)*exp(-2.0*zdnnu*r)-ddn;
423 double agrup=2.0*zup*dup;
424 double delgrup=8.0*(zup*zup*zup)*dup*dup;
425 double uplap=4.0*zup*dup*(zup-1.0/r);
426 double agrdn=2.0*zdn*ddn;
427 double delgrdn=8.0*(zdn*zdn*zdn)*ddn*ddn;
428 double dnlap=4.0*zdn*ddn*(zdn-1.0/r);
429 double d=dup+ddn;
430 double agrad=2.0*(zup*dup+zdn*ddn);
431 double delgrad=4.0*agrad*(zup*zup*dup+zdn*zdn*ddn);
432#if USE_ORIGINAL_CODE
433 double exlsd;
434 double vxuplsd;
435 double vxdnlsd;
436 double exclsd;
437 double vxcuplsd;
438 double vxcdnlsd;
439 double expw91,vxuppw91,vxdnpw91,ecpw91;
440 double expbe,vxuppbe,vxdnpbe,ecpbe;
441 double eclsd, vcuplsd, vcdnlsd, vcuppw91, vcdnpw91, vcuppbe, vcdnpbe;
442 int ione=1;
443 easypbe_(&dup,&agrup,&delgrup,&uplap,&ddn,&agrdn,&delgrdn,
444 &dnlap,&agrad,&delgrad,&ione,&ione,
445 &exlsd,&vxuplsd,&vxdnlsd,&eclsd,&vcuplsd,&vcdnlsd,
446 &expw91,&vxuppw91,&vxdnpw91,&ecpw91,&vcuppw91,&vcdnpw91,
447 &expbe,&vxuppbe,&vxdnpbe,&ecpbe,&vcuppbe,&vcdnpbe);
448 //sumexc=sumexc+d*(expbe+ecpbe)*wt;
449 sumexc=sumexc+d*(expbe+ecpbe)*wt;
450 // CHNG1=CHNG1+(vxdnpbe+vcdnpbe)*DELDDN*WT/DELZDN
451#endif
452#if USE_MPQC_CODE
453 PointOutputData od;
454 id.a.rho = dup;
455 id.a.gamma = agrup*agrup;
456 id.b.rho = ddn;
457 id.b.gamma = agrdn*agrdn;
458 id.gamma_ab = 0.5*(agrad*agrad-id.a.gamma-id.b.gamma);
459 if (id.gamma_ab > sqrt(id.a.gamma*id.b.gamma))
460 id.gamma_ab = sqrt(id.a.gamma*id.b.gamma);
461 if (id.gamma_ab < -sqrt(id.a.gamma*id.b.gamma))
462 id.gamma_ab = -sqrt(id.a.gamma*id.b.gamma);
463 if (id.gamma_ab < -0.5*(id.a.gamma*id.b.gamma))
464 id.gamma_ab = -0.5*(id.a.gamma*id.b.gamma);
465 id.compute_derived(1, valtest->need_density_gradient(),false);
466 valtest->point(id,od);
467 mpqc_sumexc += od.energy*wt;
468#endif
469// cout << scprintf("d = %12.8f wt = %12.8f OK = %12.8f MPQC = %12.8f",
470// d, wt, expbe, od.energy/d) << endl;
471 }
472 if(kdif==1) {
473 sumexco=sumexc;
474 }
475 }
476 // CHNG: DIRECT XC ENERGY INCREMENT
477 // IF THE FUNCTIONAL DERIVATIVE IS CORRECT, THEN CHNG1=CHNG
478// CHNG=(sumEXC-sumEXCO)/DELZDN
479// PRINT 200,FUP,FDN,ZUP,ZDN,sumEXC,CHNG1,chng
480#if USE_ORIGINAL_CODE
481 cout << scprintf("orig %3.1f %3.1f %3.1f %3.1f %16.12f",
482 fup,fdn,zup,zdn,sumexc) << endl;
483#endif
484#if USE_MPQC_CODE
485 cout << scprintf("mpqc %3.1f %3.1f %3.1f %3.1f %16.12f",
486 fup,fdn,zup,zdn,mpqc_sumexc) << endl;
487#endif
488 }
489}
490
491int
492main(int argc, char**argv)
493{
494
495
496 ExEnv::init(argc, argv);
497
498 Ref<MessageGrp> grp;
499#if defined(HAVE_MPI) && defined(ALWAYS_USE_MPI)
500 grp = new MPIMessageGrp(&argc, &argv);
501 MessageGrp::set_default_messagegrp(grp);
502#endif
503
504#if 0
505#ifdef HAVE_FEENABLEEXCEPT
506// this uses a glibc extension to trap on individual exceptions
507# ifdef FE_DIVBYZERO
508 feenableexcept(FE_DIVBYZERO);
509# endif
510# ifdef FE_INVALID
511 feenableexcept(FE_INVALID);
512# endif
513# ifdef FE_OVERFLOW
514 feenableexcept(FE_OVERFLOW);
515# endif
516#endif
517#endif
518
519#ifdef HAVE_FEDISABLEEXCEPT
520// this uses a glibc extension to not trap on individual exceptions
521# ifdef FE_UNDERFLOW
522 fedisableexcept(FE_UNDERFLOW);
523# endif
524# ifdef FE_INEXACT
525 fedisableexcept(FE_INEXACT);
526# endif
527#endif
528
529 int i;
530 const char *input = (argc > 1)? argv[1] : SRCDIR "/dfttest.in";
531
532 // open keyval input
533 Ref<KeyVal> keyval(new ParsedKeyVal(input));
534
535 cout << "=========== Value f Tests ===========" << endl;
536 int nvaltest = keyval->count("valtest");
537 for (i=0; i<nvaltest; i++) {
538 Ref<DenFunctional> valtest;
539 valtest << keyval->describedclassvalue("valtest", i);
540 if (valtest.nonnull()) valtest->print();
541 do_valtest(valtest);
542 }
543
544 Ref<Wavefunction> dft; dft << keyval->describedclassvalue("dft");
545 if (dft.nonnull()) {
546 cout << "=========== FD dE/dx Tests ===========" << endl;
547 fd_e_test(dft);
548 }
549
550 cout << "=========== FD df/drho Tests ===========" << endl;
551 Ref<DenFunctional> funcs[] = {
552 new PBECFunctional,
553 new PW91CFunctional,
554 new PW91XFunctional,
555 new PBEXFunctional,
556 new PW92LCFunctional,
557 new mPW91XFunctional(mPW91XFunctional::B88),
558 new mPW91XFunctional(mPW91XFunctional::PW91),
559 new mPW91XFunctional(mPW91XFunctional::mPW91),
560 new SlaterXFunctional,
561 new Becke88XFunctional,
562 new VWN1LCFunctional(1),
563 new VWN1LCFunctional,
564 new VWN2LCFunctional,
565 new VWN3LCFunctional,
566 new VWN4LCFunctional,
567 new VWN5LCFunctional,
568 new PZ81LCFunctional,
569 new P86CFunctional,
570 new NewP86CFunctional,
571 new XalphaFunctional,
572 new LYPCFunctional,
573 new PW86XFunctional,
574 new G96XFunctional,
575 0
576 };
577 const int maxerr = 1000;
578 int errcount[maxerr];
579 for (i=0; funcs[i]; i++) {
580 cout << "-----------------"
581 << funcs[i]->class_name()
582 << "-----------------"
583 << endl;
584 int nerr = funcs[i]->test();
585 if (i<maxerr) errcount[i] = nerr;
586 else { cout << "dfttest: maxerr exceeded" << endl; abort(); }
587 }
588 cout << "-------------- ERROR RESULTS --------------" << endl;
589 for (int i=0; funcs[i]; i++) {
590 cout << funcs[i]->class_name() << ": " << errcount[i];
591 if (errcount[i] == 0) cout << " (OK)";
592 cout << endl;
593 }
594
595 Ref<Molecule> mol; mol << keyval->describedclassvalue("molecule");
596 if (mol.nonnull()) {
597 cout << "=========== FD Weight Tests ===========" << endl;
598 Ref<IntegrationWeight> weights[] = {
599 new BeckeIntegrationWeight,
600 0
601 };
602 for (i=0; weights[i]; i++) {
603 cout << "-----------------"
604 << weights[i]->class_name()
605 << "-----------------"
606 << endl;
607 weights[i]->init(mol,1.0e-8);
608 weights[i]->test();
609 }
610 }
611
612 Ref<DenFunctional> functional;
613 functional << keyval->describedclassvalue("functional");
614 Ref<Wavefunction> wfn;
615 wfn << keyval->describedclassvalue("wfn");
616 if (functional.nonnull() && wfn.nonnull()) {
617 cout << "=========== FD df/dx Tests ===========" << endl;
618 fd_test(functional, wfn);
619 }
620
621 return 0;
622}
623
624/////////////////////////////////////////////////////////////////////////////
625
626// Local Variables:
627// mode: c++
628// c-file-style: "CLJ-CONDENSED"
629// End:
Note: See TracBrowser for help on using the repository browser.