← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/lib/libmints/basisset.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/libmints/basisset.cc
 
1
#include <stdexcept>
 
2
#include <cstdio>
 
3
#include <cstdlib>
 
4
#include <cmath>
 
5
 
 
6
#include <libciomr/libciomr.h>
 
7
#include <libchkpt/chkpt.hpp>
 
8
#include <psifiles.h>
 
9
 
 
10
#include "basisset.h"
 
11
#include "integral.h"
 
12
#include "symmetry.h"
 
13
 
 
14
using namespace psi;
 
15
 
 
16
BasisSet::BasisSet(Ref<Chkpt> &chkpt)
 
17
{
 
18
    max_nprimitives_ = 0;
 
19
    nshells_      = chkpt->rd_nshell();
 
20
    nprimitives_  = chkpt->rd_nprim();
 
21
    nao_          = chkpt->rd_nao();
 
22
    max_stability_index_ = 0;
 
23
    
 
24
    // Psi3 only allows either all Cartesian or all Spherical harmonic
 
25
    nbf_          = chkpt->rd_nso();
 
26
    puream_       = chkpt->rd_puream() ? true : false;
 
27
    max_am_       = chkpt->rd_max_am();
 
28
    uso2ao_       = chkpt->rd_usotao();
 
29
    
 
30
    // Allocate memory for the shells
 
31
    shells_       = new Ref<GaussianShell>[nshells_];
 
32
    
 
33
    // Initialize the shells
 
34
    initialize_shells(chkpt);
 
35
}
 
36
 
 
37
BasisSet::~BasisSet()
 
38
{
 
39
    Chkpt::free(shell_first_basis_function_);
 
40
    Chkpt::free(shell_first_ao_);
 
41
    Chkpt::free(shell_center_);
 
42
    Chkpt::free(uso2ao_);
 
43
}
 
44
 
 
45
void BasisSet::initialize_shells(Ref<Chkpt> &chkpt)
 
46
{
 
47
    // Retrieve angular momentum of each shell (1=s, 2=p, ...)
 
48
    int *shell_am = chkpt->rd_stype();
 
49
    
 
50
    // Retrieve number of primitives per shell
 
51
    int *shell_num_prims = chkpt->rd_snumg();
 
52
    
 
53
    // Retrieve exponents of primitive Gaussians
 
54
    double *exponents = chkpt->rd_exps();
 
55
    
 
56
    // Retrieve coefficients of primitive Gaussian
 
57
    double **ccoeffs = chkpt->rd_contr_full();
 
58
    
 
59
    // Retrieve pointer to first primitive in shell
 
60
    int *shell_fprim = chkpt->rd_sprim();
 
61
    
 
62
    // Retrieve pointer to first basis function in shell
 
63
    shell_first_basis_function_ = chkpt->rd_sloc_new();
 
64
    
 
65
    // Retrieve pointer to first AO in shell
 
66
    shell_first_ao_ = chkpt->rd_sloc();
 
67
    
 
68
    // Retrieve location of shells (which atom it's centered on)
 
69
    shell_center_ = chkpt->rd_snuc();
 
70
    
 
71
    // Initialize molecule, retrieves number of center and geometry
 
72
    molecule_ = new Molecule;
 
73
    molecule_->init_with_chkpt(chkpt);
 
74
    
 
75
    // Initialize SphericalTransform
 
76
    for (int i=0; i<=max_am_; ++i) {
 
77
        sphericaltransforms_.push_back(SphericalTransform(i));
 
78
    }
 
79
    
 
80
    // Initialize SOTransform
 
81
    sotransform_ = new SOTransform;
 
82
    sotransform_->init(nshells_);
 
83
    
 
84
    int *so2symblk = new int[nbf_];
 
85
    int *so2index  = new int[nbf_];
 
86
    int *sopi = chkpt->rd_sopi();
 
87
    int nirreps = chkpt->rd_nirreps();
 
88
    
 
89
    // Create so2symblk and so2index
 
90
    int ij = 0; int offset = 0;
 
91
    for (int h=0; h<nirreps; ++h) {
 
92
        for (int i=0; i<sopi[h]; ++i) {
 
93
            so2symblk[ij] = h;
 
94
            so2index[ij] = ij-offset;
 
95
            
 
96
            ij++;
 
97
        }
 
98
        offset += sopi[h];
 
99
    }
 
100
    
 
101
    // Currently all basis sets are treated as segmented contractions
 
102
    // even though GaussianShell is generalized (well not really).
 
103
    int ncontr = 1;
 
104
    int ao_start = 0;
 
105
    int puream_start = 0;
 
106
    int *sym_transform = new int[nirreps];
 
107
    
 
108
    // We need access to symmetry information found in the checkpoint file.
 
109
    Symmetry symmetry(chkpt);
 
110
    
 
111
    for (int i=0; i<nshells_; ++i) {
 
112
        int *am = new int[ncontr];
 
113
        am[0] = shell_am[i] - 1;
 
114
        int fprim = shell_fprim[i] - 1;
 
115
        int nprims = shell_num_prims[i];
 
116
        Vector3 center = molecule_->xyz(shell_center_[i] - 1);
 
117
        double **cc = new double*[nprims];
 
118
        for (int p=0; p<nprims; ++p) {
 
119
            cc[p] = new double[ncontr];
 
120
            cc[p][0] = ccoeffs[fprim+p][am[0]];
 
121
        }
 
122
        
 
123
        // Construct a new shell. GaussianShell copies the data to new memory
 
124
        shells_[i] = new GaussianShell(ncontr, nprims, &(exponents[fprim]), am, 
 
125
            puream_ ? GaussianShell::Pure : GaussianShell::Cartesian, cc, shell_center_[i]-1, center,
 
126
            puream_start);
 
127
            
 
128
        if (nprims > max_nprimitives_)
 
129
            max_nprimitives_ = nprims;
 
130
            
 
131
        for (int p=0; p<nprims; p++) {
 
132
            delete[] cc[p];
 
133
        }
 
134
        delete[] cc;
 
135
        
 
136
        // OK, for a given number of AO functions in a shell INT_NCART(am)
 
137
        // beginning at column ao_start go through all rows finding where this
 
138
        // AO function contributes to an SO.
 
139
        for (int ao = 0; ao < INT_NCART(am[0]); ++ao) {
 
140
            int aooffset = ao_start + ao;
 
141
            for (int so = 0; so < nbf_; ++so) {
 
142
                if (fabs(uso2ao_[so][aooffset]) >= 1.0e-14)
 
143
                    sotransform_->add_transform(i, so2symblk[so], so2index[so], uso2ao_[so][aooffset], ao, so);
 
144
            }
 
145
        }
 
146
        
 
147
        // Set the symmetry transform vector of the shell
 
148
        for (int j=0; j<nirreps; ++j) {
 
149
            sym_transform[j] = symmetry.trans_vec(i, j);
 
150
        }
 
151
        // The shell will copy the elements into a new vector.
 
152
        shells_[i]->set_sym_transform(nirreps, sym_transform);
 
153
    
 
154
        // Compute index of the stabilizer for the shell
 
155
        int count = 1;
 
156
        for (int g=1; g<nirreps; ++g) {
 
157
            if (i == symmetry.trans_vec(i, g)-1)
 
158
                count++;
 
159
        }
 
160
        int stab_index = nirreps / count;
 
161
        if (max_stability_index_ < stab_index)
 
162
            max_stability_index_ = stab_index;
 
163
        
 
164
        // Shift the ao starting index over to the next shell
 
165
        ao_start += INT_NCART(am[0]);
 
166
        puream_start += INT_NFUNC(puream_, am[0]);
 
167
    }
 
168
    
 
169
    delete[] sym_transform;
 
170
    delete[] so2symblk;
 
171
    delete[] so2index;
 
172
    free(sopi);
 
173
    free_block(ccoeffs);
 
174
    free(exponents);
 
175
    free(shell_am);
 
176
    free(shell_num_prims);
 
177
    free(shell_fprim);
 
178
}
 
179
 
 
180
void BasisSet::print(FILE *out) const
 
181
{
 
182
    fprintf(out, "  Basis Set\n");
 
183
    fprintf(out, "    Number of shells: %d\n", nshell());
 
184
    fprintf(out, "    Number of basis function: %d\n", nbf());
 
185
    fprintf(out, "    Number of Cartesian functions: %d\n", nao());
 
186
    fprintf(out, "    Spherical Harmonics?: %s\n", has_puream() ? "true" : "false");
 
187
    fprintf(out, "    Max angular momentum: %d\n\n", max_am());
 
188
    
 
189
    fprintf(out, "    Shells:\n\n");
 
190
    for (int s=0; s<nshell(); ++s)
 
191
        shells_[s]->print(out);
 
192
}
 
193
 
 
194
Ref<GaussianShell>& BasisSet::shell(int si) const
 
195
{
 
196
    #ifdef DEBUG
 
197
    assert(si < nshell());
 
198
    #endif
 
199
    return shells_[si];
 
200
}