Skip to content
This repository has been archived by the owner on Apr 14, 2023. It is now read-only.

Commit

Permalink
Faster and cleaner code for bond orders
Browse files Browse the repository at this point in the history
  • Loading branch information
@tovrstra
tovrstra committed Aug 3, 2010
1 parent 3045230 commit 3ca13a0
Show file tree
Hide file tree
Showing 13 changed files with 445 additions and 94 deletions.
2 changes: 1 addition & 1 deletion cleanfiles.sh
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/bin/bash
echo Removing redundant files
for i in `find hipart scripts | egrep "\.pyc$|\.py~$|\.pyc~$|\.bak|\.so$"` ; do rm -v ${i}; done
for i in `find hipart scripts | egrep "\.pyc$|\.py~$|\.pyc~$|\.bak$|\.so$|extmodule\.c$"` ; do rm -v ${i}; done
rm -vr python-build-stamp-*
rm -vr HiPart.egg-info

Expand Down
120 changes: 69 additions & 51 deletions hipart/cache.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
# --


# TODO: also evaluate overlap matrices in the basis of contracted gaussians
# TODO: Support for Gaussian/GAMESS wfn files.
# TODO: Extend hi-atomdb.py to work with GAMESS
# TODO: Implement the potential generated by an atomic density, and evaluate it
Expand All @@ -35,12 +34,13 @@

from hipart.log import log
from hipart.io import dump_atom_scalars, dump_atom_vectors, dump_atom_matrix, \
dump_atom_fields
dump_atom_fields, dump_overlap_matrices
from hipart.fit import ESPCostFunction
from hipart.lebedev_laikov import get_grid, get_atom_grid, integrate_lebedev
from hipart.atoms import AtomTable
from hipart.grids import Grid
from hipart.spline import get_radial_weights_log, CubicSpline
from hipart.gint import dmat_to_full

from molmod import Rotation, angstrom
from molmod.periodic import periodic
Expand Down Expand Up @@ -142,7 +142,7 @@ def do_atgrids(self):
prefix = os.path.join(workdir, "atom%05i" % i)
atgrid = Grid.from_prefix(prefix)
rs = self.get_rs(i)
num_points = len(self.context.lebedev_xyz.shape)*len(rs)
num_points = len(self.context.lebedev_xyz)*len(rs)
if atgrid is None or atgrid.points.shape != (num_points,3):
center = molecule.coordinates[i]
points = get_atom_grid(self.context.lebedev_xyz, center, rs)
Expand Down Expand Up @@ -582,7 +582,7 @@ def do_atgrids_orbitals(self):
self.context.wavefn.compute_orbitals(self.atgrids[i])
pb()

@OnlyOnce("Atomic overlap matrices")
@OnlyOnce("Atomic overlap matrices (in orbital basis)")
def do_atgrids_overlap_matrix_orb(self):
# Note that the overlap matrices are computed in the basis of the
# orbitals. Each kind of overlap matrix is thus computed in the basis
Expand Down Expand Up @@ -634,21 +634,44 @@ def do_one_kind(kind):
pb()

filename = os.path.join(self.context.outdir, "%s_%s_overlap_matrices_orb.txt" % (self.prefix, kind))
f = file(filename, "w")
print >> f, "number of orbitals:", num_orbitals
print >> f, "number of atoms: ", molecule.size
for i, number_i in enumerate(molecule.numbers):
print >> f, "Atom %i: %s" % (i, periodic[number_i].symbol)
matrix = getattr(self.atgrids[i], "%s_overlap_matrix_orb" % kind)
for row in matrix:
print >> f, " ".join("% 15.10e" % value for value in row)
f.close()
overlap_matrices = [
getattr(grid, "%s_overlap_matrix_orb" % kind)
for grid in self.atgrids
]
dump_overlap_matrices(filename, overlap_matrices, molecule.numbers)
log("Written %s" % filename)

do_one_kind("alpha")
do_one_kind("beta")
do_one_kind("natural")

@OnlyOnce("Atomic overlap matrices (in basis of contracted Gaussians)")
def do_atgrids_overlap_matrix(self):
self.do_integration_weights()
self.do_atgrids()
self.do_atgrids_atweights()

molecule = self.context.wavefn.molecule
num_orbitals = self.context.wavefn.num_orbitals
pb = log.pb("Computing matrices", molecule.size)
for i in xrange(molecule.size):
pb()
atgrid = self.atgrids[i]
rs = self.get_rs(i)
rw = self.radial_weights_map[len(rs)].copy()
rw *= 4*numpy.pi
rw *= rs
rw *= rs
weights = numpy.outer(rw, self.context.lebedev_weights).ravel()
weights *= atgrid.atweights
# = numpy.zeros((num_orbitals, num_orbitals), float)
self.context.wavefn.compute_atomic_overlap(atgrid, weights, self.prefix)
pb()

filename = os.path.join(self.context.outdir, "%s_overlap_matrices.txt" % self.prefix)
overlap_matrices = [atgrid.overlap_matrix for atgrid in self.atgrids]
dump_overlap_matrices(filename, overlap_matrices, molecule.numbers)

@OnlyOnce("Bond orders and valences")
def do_bond_orders(self):
# first try to load the results from the work dir
Expand All @@ -660,56 +683,51 @@ def do_bond_orders(self):
self.valences = numpy.fromfile(valences_fn_bin)
else:
self.do_charges()
self.do_atgrids_overlap_matrix_orb()
self.do_atgrids_overlap_matrix()

self.bond_orders = numpy.zeros((molecule.size, molecule.size))
self.valences = numpy.zeros(molecule.size)
num_orbitals = self.context.wavefn.num_orbitals
n_alpha = self.context.wavefn.num_alpha
n_beta = self.context.wavefn.num_beta
n_natural = self.context.wavefn.num_natural
o_alpha = self.context.wavefn.alpha_occupations[:n_alpha]
o_beta = self.context.wavefn.beta_occupations[:n_beta]
o_natural = self.context.wavefn.natural_occupations[:n_natural]

def get_pm(om, num, occ):
pm = om[:num,:num].copy()
tmp = numpy.sqrt(occ)
pm *= tmp
pm *= tmp.reshape((-1,1))
return pm

def get_correlation(om1, om2, num, occ):
pm1 = get_pm(om1, num, occ)
pm2 = get_pm(om2, num, occ)
return numpy.dot(pm1.ravel(),pm2.ravel())
num_dof = self.context.wavefn.num_orbitals

full = numpy.zeros((num_dof, num_dof), float)
dmat_to_full(self.context.wavefn.density_matrix, full)
if self.context.wavefn.spin_density_matrix is None:
full_alpha = 0.5*full
full_beta = full_alpha
else:
full_alpha = numpy.zeros((num_dof, num_dof), float)
full_beta = numpy.zeros((num_dof, num_dof), float)
dmat_to_full(
0.5*(self.context.wavefn.density_matrix +
self.context.wavefn.spin_density_matrix), full_alpha
)
dmat_to_full(
0.5*(self.context.wavefn.density_matrix -
self.context.wavefn.spin_density_matrix), full_beta
)

pb = log.pb("Computing bond orders", (molecule.size*(molecule.size+1))/2)
for i in xrange(molecule.size):
for j in xrange(i+1):
pb()
if i==j:
# compute valence
bo = get_correlation(
self.atgrids[i].natural_overlap_matrix_orb,
self.atgrids[j].natural_overlap_matrix_orb,
n_natural, o_natural
)
self.valences[i] = 2*self.populations[i] - bo
tmp = numpy.dot(full, self.atgrids[i].overlap_matrix)
self.valences[i] = 2*self.populations[i] - (tmp*tmp.transpose()).sum()
else:
# compute bond order
bo = 2*(
get_correlation(
self.atgrids[i].alpha_overlap_matrix_orb,
self.atgrids[j].alpha_overlap_matrix_orb,
n_alpha, o_alpha
)+
get_correlation(
self.atgrids[i].beta_overlap_matrix_orb,
self.atgrids[j].beta_overlap_matrix_orb,
n_beta, o_beta
)
)
bo = (
numpy.dot(full_alpha, self.atgrids[i].overlap_matrix)*
numpy.dot(full_alpha, self.atgrids[j].overlap_matrix).transpose()
).sum()
if full_alpha is full_beta:
bo *= 2
else:
bo += (
numpy.dot(full_beta, self.atgrids[i].overlap_matrix)*
numpy.dot(full_beta, self.atgrids[j].overlap_matrix).transpose()
).sum()
bo *= 2
self.bond_orders[i,j] = bo
self.bond_orders[j,i] = bo
pb()
Expand Down
30 changes: 22 additions & 8 deletions hipart/gint/basis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,16 +160,30 @@ def from_fchk(cls, fchk):

return result

def call_gint(self, gint1_fn, weights, points):
fns = numpy.zeros(len(points), float)
retcode = gint1_fn(
weights, fns, points, self.molecule.coordinates, self.shell_types,
def call_gint_grid(self, gint_c_routine, weights, points):
result = numpy.zeros(len(points), float)
retcode = gint_c_routine(
weights, result, points, self.molecule.coordinates, self.shell_types,
self.shell_map, self.num_primitives, self.ccoeffs, self.exponents
)
if retcode == -1:
raise MemoryError("Out of memory while calling %s." % gint1_fn)
raise MemoryError("Out of memory while calling %s." % gint_c_routine)
elif retcode == -2:
raise ValueError("Unsuported shell type when calling %s." % gint1_fn)
raise ValueError("Unsuported shell type when calling %s." % gint_c_routine)
elif retcode != 0:
raise RuntimeError("Something went wrong when calling %s. Got retcode=%i." % (gint1_fn, retcode))
return fns
raise RuntimeError("Something went wrong when calling %s. Got retcode=%i." % (gint_c_routine, retcode))
return result

def call_gint_atomic_operator(self, gint_c_routine, points, weights):
result = numpy.zeros((self.num_dof, self.num_dof), float)
retcode = gint_c_routine(
result, points, weights, self.molecule.coordinates, self.shell_types,
self.shell_map, self.num_primitives, self.ccoeffs, self.exponents
)
if retcode == -1:
raise MemoryError("Out of memory while calling %s." % gint_c_routine)
elif retcode == -2:
raise ValueError("Unsuported shell type when calling %s." % gint_c_routine)
elif retcode != 0:
raise RuntimeError("Something went wrong when calling %s. Got retcode=%i." % (gint_c_routine, retcode))
return result
96 changes: 96 additions & 0 deletions hipart/gint/gint1_fn.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -341,3 +341,99 @@ int gint1_fn_dmat(double* dmat, double* density, double* points,
free(basis_fns);
return result;
}

int gint1_fn_overlap(double* overlap, double* points, double* weights,
double* centers, int* shell_types, int* shell_map, int* num_primitives,
double* ccoeffs, double* exponents, int num_overlap, int num_points,
int num_centers, int num_shells, int num_ccoeffs, int num_exponents)
{
int shell, shell_type, primitive, dof1, dof2, num_dof, num_shell_dof, result, i_point;
double *work, *out, *center, *ccoeff, *exponent, *basis_fns, *shell_fns, *fn, *overlap_element;

result = 0;

work = malloc(MAX_SHELL_DOF*sizeof(double));
CHECK_ALLOC(work);
num_dof = (int)(sqrt(num_overlap));
basis_fns = malloc(num_dof*sizeof(double));
CHECK_ALLOC(basis_fns);

// Clear the output
out = overlap;
for (dof1=0; dof1<num_overlap; dof1++) {
*out = 0.0;
out++;
}

for (i_point=0; i_point<num_points; i_point++) {
// A) clear the basis functions.
fn = basis_fns;
for (dof1=0; dof1<num_dof; dof1++) {
*fn = 0.0;
fn++;
}
// B) evaluate the basis functions in the current point.
ccoeff = ccoeffs;
exponent = exponents;
shell_fns = basis_fns;
for (shell=0; shell<num_shells; shell++) {
center = centers + (3*shell_map[shell]);
shell_type = shell_types[shell];
CHECK_SHELL(shell_type);
for (primitive=0; primitive<num_primitives[shell]; primitive++) {
gint1_fn_dispatch(shell_type, center, *exponent, points, work);
//printf("shell_type=%d primitive=%d exponent=%f\n", shell_type, primitive, *exponent);
out = work;
exponent++;
fn = shell_fns;
if (shell_type==-1) {
*fn += (*out)*(*ccoeff);
//printf("out=%f ccoeff=%f fn=%f\n", *out, *ccoeff, (*out)*(*ccoeff));
fn++;
out++;
ccoeff++;
num_shell_dof = 3;
} else if (shell_type > 0) {
num_shell_dof = ((shell_type+1)*(shell_type+2))/2;
} else {
num_shell_dof = -2*shell_type+1;
}
for (dof1=0; dof1<num_shell_dof; dof1++) {
*fn += (*out)*(*ccoeff);
//printf("out=%f ccoeff=%f fn=%f\n", *out, *ccoeff, (*out)*(*ccoeff));
fn++;
out++;
}
ccoeff++;
}
if (shell_type==-1) {
num_shell_dof = 4;
} else if (shell_type > 0) {
num_shell_dof = ((shell_type+1)*(shell_type+2))/2;
} else {
num_shell_dof = -2*shell_type+1;
}
shell_fns += num_shell_dof;
}
//printf("\n");
// C) Multiply overlap elements with the wieght and add to the atomic
// overlap matrix
overlap_element = overlap;
for (dof1=0; dof1<num_dof; dof1++) {
for (dof2=0; dof2<num_dof; dof2++) {
*overlap_element += (*weights)*basis_fns[dof1]*basis_fns[dof2];
//printf("dof1=%d dof2=%d basis1=%f basis2=%f overlap_element=%f weight=%f\n", dof1, dof2, basis_fns[dof1], basis_fns[dof2], *overlap_element, *weight);
overlap_element++;
}
}
// D) Prepare for next iteration
points += 3;
weights++;
//printf("\n");
}

EXIT:
free(work);
free(basis_fns);
return result;
}
20 changes: 20 additions & 0 deletions hipart/gint/gint1_fn.pyf.inc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,3 +107,23 @@
integer intent(hide), depend(ccoeffs) :: num_ccoeffs=len(ccoeffs)
integer intent(hide), depend(exponents) :: num_exponents=len(exponents)
end function gint1_fn_dmat

integer function gint1_fn_overlap(overlap, points, weights, centers, shell_types, shell_map, num_primitives, ccoeffs, exponents, num_overlap, num_points, num_centers, num_shells, num_ccoeffs, num_exponents)
intent(c) gint1_fn_overlap
intent(c)
double precision intent(inout) :: overlap(num_overlap)
double precision intent(in) :: points(num_points,3)
double precision intent(in) :: weights(num_points)
double precision intent(in) :: centers(num_centers,3)
integer intent(int) :: shell_types(num_shells)
integer intent(int) :: shell_map(num_shells)
integer intent(int) :: num_primitives(num_shells)
double precision intent(in) :: ccoeffs(num_ccoeffs)
double precision intent(in) :: exponents(num_exponents)
integer intent(hide), depend(overlap) :: num_overlap=len(overlap)
integer intent(hide), depend(points) :: num_points=len(points)
integer intent(hide), depend(centers) :: num_centers=len(centers)
integer intent(hide), depend(shell_types) :: num_shells=len(shell_types)
integer intent(hide), depend(ccoeffs) :: num_ccoeffs=len(ccoeffs)
integer intent(hide), depend(exponents) :: num_exponents=len(exponents)
end function gint1_fn_overlap
Loading

0 comments on commit 3ca13a0

Please sign in to comment.