← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/lib/libmints/molecule.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/molecule.cc
 
1
#include <cmath>
 
2
#include <cstdio>
 
3
 
 
4
#include <libmints/molecule.h>
 
5
 
 
6
#include <masses.h>
 
7
#include <physconst.h>
 
8
 
 
9
using namespace std;
 
10
using namespace psi;
 
11
 
 
12
Molecule::Molecule():
 
13
    natoms_(0), nirreps_(0)
 
14
{
 
15
    
 
16
}
 
17
 
 
18
Molecule::~Molecule()
 
19
{
 
20
    clear();
 
21
}
 
22
 
 
23
void Molecule::clear()
 
24
{
 
25
    natoms_ = 0;
 
26
    nirreps_ = 0;
 
27
    atoms_.empty();
 
28
}
 
29
 
 
30
void Molecule::add_atom(int Z, double x, double y, double z,
 
31
                        const char *label, double mass,
 
32
                        int have_charge, double charge)
 
33
{
 
34
    atom_info info;
 
35
    
 
36
    info.x = x;
 
37
    info.y = y;
 
38
    info.z = z;
 
39
    info.Z = Z;
 
40
    info.charge = charge;
 
41
    info.label = label;
 
42
    info.mass = mass;
 
43
    
 
44
    natoms_++;
 
45
    atoms_.push_back(info);
 
46
}
 
47
 
 
48
double Molecule::mass(int atom) const
 
49
{
 
50
    if (atoms_[atom].mass != 0.0)
 
51
        return atoms_[atom].mass;
 
52
        
 
53
    return an2masses[atoms_[atom].Z];
 
54
}
 
55
 
 
56
const std::string Molecule::label(int atom) const
 
57
{
 
58
    return atoms_[atom].label; 
 
59
}
 
60
 
 
61
int Molecule::atom_at_position(double *xyz, double tol) const
 
62
{
 
63
    Vector3 b(xyz);
 
64
    for (int i=0; i < natom(); ++i) {
 
65
        Vector3 a(atoms_[i].x, atoms_[i].y, atoms_[i].z);
 
66
        if (b.distance(a) < tol)
 
67
            return i;
 
68
    }
 
69
    return -1;
 
70
}
 
71
 
 
72
Vector3 Molecule::center_of_mass() const
 
73
{
 
74
    Vector3 ret;
 
75
    double total_m;
 
76
    
 
77
    ret = 0.0;
 
78
    total_m = 0.0;
 
79
    
 
80
    for (int i=0; i<natom(); ++i) {
 
81
        double m = mass(i);
 
82
        ret += m * xyz(i);
 
83
        total_m += m;
 
84
    }
 
85
    
 
86
    ret *= 1.0/total_m;
 
87
    
 
88
    return ret;
 
89
}
 
90
 
 
91
double Molecule::nuclear_repulsion_energy()
 
92
{
 
93
    double e=0.0;
 
94
    
 
95
    for (int i=1; i<natom(); ++i) {
 
96
        for (int j=0; j<i; ++j) {
 
97
            e += Z(i) * Z(j) / (xyz(i).distance(xyz(j)));
 
98
        }
 
99
    }
 
100
    
 
101
    return e;
 
102
}
 
103
 
 
104
double* Molecule::nuclear_repulsion_energy_deriv1()
 
105
{
 
106
    double *de = new double[3*natom()];
 
107
    
 
108
    memset(de, 0, sizeof(double)*3*natom());
 
109
    for (int i=1; i<natom(); ++i) {
 
110
        for (int j=0; j<natom(); ++j) {
 
111
            if (i != j) {
 
112
                double temp = pow((xyz(i).distance(xyz(j))), 3.0);
 
113
                de[3*i+0] -= (x(i) - x(j)) * Z(i) * Z(j) / temp;
 
114
                de[3*i+1] -= (y(i) - y(j)) * Z(i) * Z(j) / temp;
 
115
                de[3*i+2] -= (z(i) - z(j)) * Z(i) * Z(j) / temp;
 
116
            }
 
117
        }
 
118
    }
 
119
    
 
120
    return de;
 
121
}
 
122
 
 
123
void Molecule::translate(const Vector3& r)
 
124
{
 
125
    for (int i=0; i<natom(); ++i) {
 
126
        atoms_[i].x += r[0];
 
127
        atoms_[i].y += r[1];
 
128
        atoms_[i].z += r[2];
 
129
    }
 
130
}
 
131
 
 
132
void Molecule::move_to_com()
 
133
{
 
134
    Vector3 com = -center_of_mass();
 
135
    translate(com);
 
136
}
 
137
 
 
138
void Molecule::init_with_chkpt(Ref<PSIO> &psio)
 
139
{
 
140
    // User sent a psio object. Create a chkpt object based on it.
 
141
    Ref<Chkpt> chkpt(new Chkpt(psio.pointer(), PSIO_OPEN_OLD));
 
142
    init_with_chkpt(chkpt);
 
143
}
 
144
 
 
145
void Molecule::init_with_chkpt(Ref<Chkpt> &chkpt)
 
146
{
 
147
    int atoms = chkpt->rd_natom();
 
148
    double *zvals = chkpt->rd_zvals();
 
149
    double **geom = chkpt->rd_geom();
 
150
    
 
151
    for (int i=0; i<atoms; ++i) {
 
152
        add_atom((int)zvals[i], geom[i][0], geom[i][1], geom[i][2], atomic_labels[(int)zvals[i]], an2masses[(int)zvals[i]]);
 
153
    }
 
154
    
 
155
    nirreps_ = chkpt->rd_nirreps();
 
156
    
 
157
    Chkpt::free(zvals);
 
158
    Chkpt::free(geom);
 
159
}
 
160
 
 
161
void Molecule::print(FILE *out)
 
162
{
 
163
    if (natom()) {
 
164
        fprintf(out,"       Center              X                  Y                   Z\n");
 
165
        fprintf(out,"    ------------   -----------------  -----------------  -----------------\n");
 
166
        
 
167
        for(int i = 0; i < natom(); ++i){
 
168
            Vector3 geom = xyz(i);
 
169
            fprintf(out, "    %12s ",label(i).c_str()); fflush(out);
 
170
            for(int j = 0; j < 3; j++)
 
171
                fprintf(out, "  %17.12f", geom[j]);
 
172
            fprintf(out,"\n");
 
173
        }
 
174
        fprintf(out,"\n");
 
175
        fflush(out);
 
176
    }
 
177
}
 
178
 
 
179
SimpleVector Molecule::nuclear_dipole_contribution()
 
180
{
 
181
    SimpleVector ret(3);
 
182
    
 
183
    for(int i=0; i<natom(); ++i) {
 
184
        Vector3 geom = xyz(i);
 
185
        ret[0] += Z(i) * geom[0];
 
186
        ret[1] += Z(i) * geom[1];
 
187
        ret[2] += Z(i) * geom[2];
 
188
    }
 
189
        
 
190
    return ret;
 
191
}
 
192
 
 
193
SimpleVector Molecule::nuclear_quadrupole_contribution()
 
194
{
 
195
    SimpleVector ret(6);
 
196
    double xx, xy, xz, yy, yz, zz;
 
197
    
 
198
    xx = xy = xz = yy = yz = zz = 0.0;
 
199
    
 
200
    for (int i=0; i<natom(); ++i) {
 
201
        Vector3 geom = xyz(i);
 
202
        ret[0] += Z(i) * geom[0] * geom[0]; // xx
 
203
        ret[1] += Z(i) * geom[0] * geom[1]; // xy
 
204
        ret[2] += Z(i) * geom[0] * geom[2]; // xz
 
205
        ret[3] += Z(i) * geom[1] * geom[1]; // yy
 
206
        ret[4] += Z(i) * geom[1] * geom[2]; // yz
 
207
        ret[5] += Z(i) * geom[2] * geom[2]; // zz
 
208
    }
 
209
            
 
210
    return ret;
 
211
}