1 //
2 // densval.cc
3 //
4 // Copyright (C) 1996 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
28 #include <util/misc/formio.h>
29 #include <util/keyval/keyval.h>
30 #include <chemistry/qc/basis/basis.h>
31 #include <chemistry/qc/basis/petite.h>
32
33 #include <chemistry/qc/wfn/wfn.h>
34
35 using namespace sc;
36
37 // Function for returning electron charge density at a point
density(const SCVector3 & r)38 double Wavefunction::density(const SCVector3&r)
39 {
40 int nbasis = basis()->nbasis();
41 if (!bs_values) bs_values=new double[nbasis];
42
43 // compute the basis set values
44 GaussianBasisSet::ValueData *valdat
45 = new GaussianBasisSet::ValueData(basis(), integral_);
46 basis()->values(r,valdat,bs_values);
47 delete valdat;
48
49 //for (int i=0; i<nbasis; i++)
50 // ExEnv::out0() << indent
51 // << scprintf("bs_values[%d] = %12.8f\n",i,bs_values[i]);
52
53 // Assuming this will be called many times for the same wavefunction,
54 // it is more efficient to force the computation of the natural
55 // orbitals now and use them.
56 // get the natural orbitals and density
57 RefSCMatrix nos
58 = integral()->petite_list()->evecs_to_AO_basis(natural_orbitals());
59 RefDiagSCMatrix nd = natural_density();
60
61 // loop over natural orbitals adding contributions to elec_density
62 double elec_density=0.0;
63 for (int i=0; i<nbasis; i++) {
64 double tmp = 0.0;
65 for (int j=0; j<nbasis; j++) {
66 tmp += nos.get_element(j,i)*bs_values[j];
67 }
68 elec_density += nd.get_element(i)*tmp*tmp;
69 }
70
71 return elec_density;
72 }
73
74 // Function for returning electron charge density at a point.
75 // The grad at that point is also computed and put into double grad[3].
density_gradient(const SCVector3 & r,double * grad)76 double Wavefunction::density_gradient(const SCVector3&r,double*grad)
77 {
78 int nbasis = basis()->nbasis();
79 if (!bs_values) bs_values=new double[nbasis];
80 if (!bsg_values) bsg_values=new double[nbasis*3];
81
82 // compute the grad values and get the basis set values at the
83 // same time
84 GaussianBasisSet::ValueData *valdat
85 = new GaussianBasisSet::ValueData(basis(), integral_);
86 basis()->grad_values(r,valdat,bsg_values,bs_values);
87 delete valdat;
88
89 //for (int i=0; i<nbasis; i++)
90 // ExEnv::out0() << indent
91 // << scprintf("bs_values[%d] = % 12.8f\n",i,bs_values[i]);
92
93 // get the natural orbitals and density
94 RefSCMatrix nos
95 = integral()->petite_list()->evecs_to_AO_basis(natural_orbitals());
96 RefDiagSCMatrix nd = natural_density();
97
98 // loop over natural orbitals adding contributions to elec_density
99 double elec_density=0.0;
100 grad[0] = grad[1] = grad[2] = 0.0;
101 for (int i=0; i<nbasis; i++) {
102 double tmp = 0.0;
103 int j;
104 for (j=0; j<nbasis; j++) {
105 tmp += nos.get_element(j,i)*bs_values[j];
106 }
107 elec_density += nd.get_element(i)*tmp*tmp;
108 double tmpg[3];
109 tmpg[0] = tmpg[1] = tmpg[2] = 0.0;
110 for (j=0; j<nbasis; j++) {
111 tmpg[0] += nos.get_element(j,i)*bsg_values[j*3+0];
112 tmpg[1] += nos.get_element(j,i)*bsg_values[j*3+1];
113 tmpg[2] += nos.get_element(j,i)*bsg_values[j*3+2];
114 }
115 grad[0] += nd.get_element(i)*tmpg[0]*tmp;
116 grad[1] += nd.get_element(i)*tmpg[1]*tmp;
117 grad[2] += nd.get_element(i)*tmpg[2]*tmp;
118 }
119
120 grad[0] *= 2.0;
121 grad[1] *= 2.0;
122 grad[2] *= 2.0;
123
124 return elec_density;
125 }
126
127 /////////////////////////////////////////////////////////////////////////////
128
129 // Local Variables:
130 // mode: c++
131 // c-file-style: "CLJ"
132 // End:
133