Coverage for cclib/method/volume.py : 78%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2#
3# Copyright (c) 2020, the cclib development team
4#
5# This file is part of cclib (http://cclib.github.io) and is distributed under
6# the terms of the BSD 3-Clause License.
8"""Calculation methods related to volume based on cclib data."""
10import copy
12import numpy
14from cclib.parser.utils import convertor
15from cclib.parser.utils import find_package
17""" In the dictionary sym2powerlist below, each element is a list that contain the combinations of
18 powers that are applied to x, y, and z in the equation for the gaussian primitives --
19 \psi (x, y, z) = x^a * y^b * z^c * exp(-\lambda * r^2)
20"""
21sym2powerlist = {
22 "S": [(0, 0, 0)],
23 "P": [(1, 0, 0), (0, 1, 0), (0, 0, 1)],
24 "D": [(2, 0, 0), (0, 2, 0), (0, 0, 2), (1, 1, 0), (0, 1, 1), (1, 0, 1)],
25 "F": [
26 (3, 0, 0),
27 (2, 1, 0),
28 (2, 0, 1),
29 (1, 2, 0),
30 (1, 1, 1),
31 (1, 0, 2),
32 (0, 3, 0),
33 (0, 2, 1),
34 (0, 1, 2),
35 (0, 0, 3),
36 ],
37}
40_found_pyquante = find_package("PyQuante")
41if _found_pyquante:
42 from PyQuante.CGBF import CGBF
44 def getbfs(ccdata):
45 from cclib.bridge import cclib2pyquante
47 pymol = cclib2pyquante.makepyquante(ccdata)
49 bfs = []
50 for i, atom in enumerate(pymol):
51 bs = ccdata.gbasis[i]
52 for sym, prims in bs:
53 for power in sym2powerlist[sym]:
54 bf = CGBF(atom.pos(), power)
55 for expnt, coef in prims:
56 bf.add_primitive(expnt, coef)
57 bf.normalize()
58 bfs.append(bf)
60 del cclib2pyquante
62 return bfs
64 def pyamp(bfs, bs, points):
65 """Wrapper for evaluating basis functions at one or more grid points.
67 Parameters
68 ----------
69 bfs : list
70 List of PyQuante 1 `CGBFs`s (contracted Gaussian basis functions).
71 bs : int
72 Index into the list of CGBFs for the basis function to evaluate.
73 points : numpy.ndarray
74 An [n, 3] array of `n` Cartesian grid points on which to evaluate the basis function.
76 Returns
77 -------
78 out : numpy.ndarray
79 An [n, ] array of the requested basis function's value on each grid point.
80 """
81 mesh_vals = numpy.zeros(len(points))
82 for i in range(len(points)):
83 mesh_vals[i] = bfs[bs].amp(points[i][0], points[i][1], points[i][2])
84 return mesh_vals
87_found_pyquante2 = find_package("pyquante2")
88if _found_pyquante2:
89 from pyquante2 import cgbf
91 def getbfs(ccdata):
92 bfs = []
93 # `atom` is instance of pyquante2.geo.atom.atom class.
94 for i, atom in enumerate(ccdata.atomcoords[-1]):
95 # `basis` is basis coefficients stored in ccData.
96 basis = ccdata.gbasis[i]
97 for sym, primitives in basis:
98 # `sym` is S, P, D, F and is used as key here.
99 for power in sym2powerlist[sym]:
100 exponentlist = []
101 coefficientlist = []
103 for exponents, coefficients in primitives:
104 exponentlist.append(exponents)
105 coefficientlist.append(coefficients)
107 basisfunction = cgbf(
108 convertor(atom, "Angstrom", "bohr"),
109 powers=power,
110 exps=exponentlist,
111 coefs=coefficientlist,
112 )
113 basisfunction.normalize()
114 bfs.append(basisfunction)
116 return bfs
118 def pyamp(bfs, bs, points):
119 """Wrapper for evaluating basis functions at one or more grid points.
121 Parameters
122 ----------
123 bfs : list
124 List of pyquante2 `cgbf`s (contracted Gaussian basis functions).
125 bs : int
126 Index into the list of CGBFs for the basis function to evaluate.
127 points : numpy.ndarray
128 An [n, 3] array of `n` Cartesian grid points on which to evaluate the basis function.
130 Returns
131 -------
132 out : numpy.ndarray
133 An [n, ] array of the requested basis function's value on each grid point.
134 """
135 return bfs[bs].mesh(points)
138_found_pyvtk = find_package("pyvtk")
139if _found_pyvtk:
140 from pyvtk import *
141 from pyvtk.DataSetAttr import *
144def _check_pyquante():
145 if (not _found_pyquante) and (not _found_pyquante2):
146 raise ImportError(
147 "You must install `pyquante2` or `PyQuante` to use this function."
148 )
151def _check_pyvtk(found_pyvtk):
152 if not found_pyvtk:
153 raise ImportError("You must install `pyvtk` to use this function.")
156class Volume(object):
157 """Represent a volume in space.
159 Required parameters:
160 origin -- the bottom left hand corner of the volume
161 topcorner -- the top right hand corner
162 spacing -- the distance between the points in the cube
164 Attributes:
165 data -- a NumPy array of values for each point in the volume
166 (set to zero at initialisation)
167 numpts -- the numbers of points in the (x,y,z) directions
168 """
170 def __init__(self, origin, topcorner, spacing):
172 self.origin = numpy.asarray(origin, dtype=float)
173 self.topcorner = numpy.asarray(topcorner, dtype=float)
174 self.spacing = numpy.asarray(spacing, dtype=float)
175 self.numpts = []
176 for i in range(3):
177 self.numpts.append(
178 int((self.topcorner[i] - self.origin[i]) / self.spacing[i] + 1)
179 )
180 self.data = numpy.zeros(tuple(self.numpts), "d")
182 def __str__(self):
183 """Return a string representation."""
184 return "Volume %s to %s (density: %s)" % (
185 self.origin,
186 self.topcorner,
187 self.spacing,
188 )
190 def write(self, filename, fformat="Cube"):
191 """Write the volume to a file."""
193 fformat = fformat.lower()
195 writers = {
196 "vtk": self.writeasvtk,
197 "cube": self.writeascube,
198 }
200 if fformat not in writers:
201 raise RuntimeError("File format must be either VTK or Cube")
203 writers[fformat](filename)
205 def writeasvtk(self, filename):
206 _check_pyvtk(_found_pyvtk)
207 ranges = (
208 numpy.arange(self.data.shape[2]),
209 numpy.arange(self.data.shape[1]),
210 numpy.arange(self.data.shape[0]),
211 )
212 v = VtkData(
213 RectilinearGrid(*ranges),
214 "Test",
215 PointData(Scalars(self.data.ravel(), "from cclib", "default")),
216 )
217 v.tofile(filename)
219 def integrate(self, weights=None):
220 weights = numpy.ones_like(self.data) if weights is None else weights
221 assert (
222 weights.shape == self.data.shape
223 ), "Shape of weights do not match with shape of Volume data."
225 boxvol = (
226 self.spacing[0]
227 * self.spacing[1]
228 * self.spacing[2]
229 * convertor(1, "Angstrom", "bohr") ** 3
230 )
232 return numpy.sum(self.data * weights) * boxvol
234 def integrate_square(self, weights=None):
235 weights = numpy.ones_like(self.data) if weights is None else weights
237 boxvol = (
238 self.spacing[0]
239 * self.spacing[1]
240 * self.spacing[2]
241 * convertor(1, "Angstrom", "bohr") ** 3
242 )
243 return numpy.sum((self.data * weights) ** 2) * boxvol
245 def writeascube(self, filename):
246 # Remember that the units are bohr, not Angstroms
247 def convert(x):
248 return convertor(x, "Angstrom", "bohr")
250 ans = []
251 ans.append("Cube file generated by cclib")
252 ans.append("")
253 format = "%4d%12.6f%12.6f%12.6f"
254 origin = [convert(x) for x in self.origin]
255 ans.append(format % (0, origin[0], origin[1], origin[2]))
256 ans.append(format % (self.data.shape[0], convert(self.spacing[0]), 0.0, 0.0))
257 ans.append(format % (self.data.shape[1], 0.0, convert(self.spacing[1]), 0.0))
258 ans.append(format % (self.data.shape[2], 0.0, 0.0, convert(self.spacing[2])))
259 line = []
260 for i in range(self.data.shape[0]):
261 for j in range(self.data.shape[1]):
262 for k in range(self.data.shape[2]):
263 line.append(scinotation(self.data[i, j, k]))
264 if len(line) == 6:
265 ans.append(" ".join(line))
266 line = []
267 if line:
268 ans.append(" ".join(line))
269 line = []
270 with open(filename, "w") as outputfile:
271 outputfile.write("\n".join(ans))
273 def coordinates(self, indices):
274 """Return [x, y, z] coordinates that correspond to input indices"""
275 return self.origin + self.spacing * indices
278def scinotation(num):
279 """Write in scientific notation."""
280 ans = "%10.5E" % num
281 broken = ans.split("E")
282 exponent = int(broken[1])
283 if exponent < -99:
284 return " 0.000E+00"
285 if exponent < 0:
286 sign = "-"
287 else:
288 sign = "+"
289 return ("%sE%s%s" % (broken[0], sign, broken[1][-2:])).rjust(12)
292def getGrid(vol):
293 """Helper function that returns (x, y, z), each of which are numpy array of the values that
294 correspond to grid points.
296 Input:
297 vol -- Volume object (will not be altered)
298 """
299 conversion = convertor(1, "bohr", "Angstrom")
300 gridendpt = vol.topcorner + 0.5 * vol.spacing
301 x = numpy.arange(vol.origin[0], gridendpt[0], vol.spacing[0]) / conversion
302 y = numpy.arange(vol.origin[1], gridendpt[1], vol.spacing[1]) / conversion
303 z = numpy.arange(vol.origin[2], gridendpt[2], vol.spacing[2]) / conversion
305 return (x, y, z)
308def wavefunction(ccdata, volume, mocoeffs):
309 """Calculate the magnitude of the wavefunction at every point in a volume.
311 Inputs:
312 ccdata -- ccData object
313 volume -- Volume object (will not be altered)
314 mocoeffs -- molecular orbital to use for calculation; i.e. ccdata.mocoeffs[0][3]
316 Output:
317 Volume object with wavefunction at each grid point stored in data attribute
318 """
319 _check_pyquante()
320 bfs = getbfs(ccdata)
322 wavefn = copy.copy(volume)
323 wavefn.data = numpy.zeros(wavefn.data.shape, "d")
325 x, y, z = getGrid(wavefn)
326 gridpoints = numpy.asanyarray(
327 tuple((xp, yp, zp) for xp in x for yp in y for zp in z)
328 )
330 # PyQuante & pyquante2
331 for bs in range(len(bfs)):
332 if abs(mocoeffs[bs]) > 0.0:
333 wavefn.data += numpy.resize(
334 pyamp(bfs, bs, gridpoints) * mocoeffs[bs], wavefn.data.shape
335 )
337 return wavefn
340def electrondensity_spin(ccdata, volume, mocoeffslist):
341 """Calculate the magnitude of the electron density at every point in a volume for either up or down spin
343 Inputs:
344 ccdata -- ccData object
345 volume -- Volume object (will not be altered)
346 mocoeffslist -- list of molecular orbital to calculate electron density from;
347 i.e. [ccdata.mocoeffs[0][1:2]]
349 Output:
350 Volume object with wavefunction at each grid point stored in data attribute
352 Attributes:
353 coords -- the coordinates of the atoms
354 mocoeffs -- mocoeffs for all of the occupied eigenvalues
355 gbasis -- gbasis from a parser object
356 volume -- a template Volume object (will not be altered)
358 Note: mocoeffs is a list of NumPy arrays. The list will be of length 1.
359 """
360 assert (
361 len(mocoeffslist) == 1
362 ), "mocoeffslist input to the function should have length of 1."
363 _check_pyquante()
364 bfs = getbfs(ccdata)
366 density = copy.copy(volume)
367 density.data = numpy.zeros(density.data.shape, "d")
369 x, y, z = getGrid(density)
370 gridpoints = numpy.asanyarray(
371 tuple((xp, yp, zp) for xp in x for yp in y for zp in z)
372 )
374 # For occupied orbitals
375 # `mocoeff` and `gbasis` in ccdata object is ordered in a way `homos` can specify which orbital
376 # is the highest lying occupied orbital in mocoeff and gbasis.
377 for mocoeffs in mocoeffslist:
378 for mocoeff in mocoeffs:
379 wavefn = numpy.zeros(density.data.shape, "d")
380 for bs in range(len(bfs)):
381 if abs(mocoeff[bs]) > 0.0:
382 wavefn += numpy.resize(
383 pyamp(bfs, bs, gridpoints) * mocoeff[bs], density.data.shape
384 )
385 density.data += wavefn ** 2
386 return density
389def electrondensity(ccdata, volume, mocoeffslist):
390 """Calculate the magnitude of the total electron density at every point in a volume.
392 Inputs:
393 ccdata -- ccData object
394 volume -- Volume object (will not be altered)
395 mocoeffslist -- list of molecular orbital to calculate electron density from;
396 i.e. [ccdata.mocoeffs[0][1:2]]
398 Output:
399 Volume object with wavefunction at each grid point stored in data attribute
401 Attributes:
402 coords -- the coordinates of the atoms
403 mocoeffs -- mocoeffs for all of the occupied eigenvalues
404 gbasis -- gbasis from a parser object
405 volume -- a template Volume object (will not be altered)
407 Note: mocoeffs is a list of NumPy arrays. The list will be of length 1
408 for restricted calculations, and length 2 for unrestricted.
409 """
411 if len(mocoeffslist) == 2:
412 return electrondensity_spin(
413 ccdata, volume, [mocoeffslist[0]]
414 ) + electrondensity_spin(ccdata, volume, [mocoeffslist[1]])
415 else:
416 edens = electrondensity_spin(ccdata, volume, [mocoeffslist[0]])
417 edens.data *= 2
418 return edens
421def read_from_cube(filepath):
422 """Read data from cube files and construct volume object
424 Specification of cube file format is described in http://paulbourke.net/dataformats/cube/
426 Input:
427 filepath -- path to the cube file
429 Output:
430 vol -- Volume object filled with data from cube file
431 """
433 with open(filepath) as f:
434 lines = f.readlines()
436 # First two lines are comments
437 # Lines 3-6 specify the grid in Cartesian coordinates
438 # Line 3 -- [Number of atoms] [Origin x] [Origin y] [Origin z]
439 natom = (int)(lines[2].split()[0])
440 originx, originy, originz = numpy.asanyarray(lines[2].split()[1:], dtype=float)
442 # Line 4, 5, 6 -- [Number of Grid Points] [Spacing x] [Spacing y], [Spacing z]
443 ngridx, spacingx = numpy.asanyarray(lines[3].split(), dtype=float)[[0, 1]]
444 ngridy, spacingy = numpy.asanyarray(lines[4].split(), dtype=float)[[0, 2]]
445 ngridz, spacingz = numpy.asanyarray(lines[5].split(), dtype=float)[[0, 3]]
447 # Line 7, 8, 9 -- Optional Information not used yet by cclib Volume
448 # [Atomic number] [Charge] [Position x] [Position y] [Position z]
449 skiplines = 6
450 while len(lines[skiplines].split()) == 5:
451 skiplines += 1
453 # Line 10+ (or 7+ if atomic coordinates data are not present)
454 tmp = []
455 for line in lines[skiplines:]:
456 tmp.extend(line.split())
457 tmp = numpy.asanyarray(tmp, dtype=float)
459 # Initialize volume object
460 vol = Volume(
461 (
462 convertor(originx, "bohr", "Angstrom"),
463 convertor(originy, "bohr", "Angstrom"),
464 convertor(originz, "bohr", "Angstrom"),
465 ),
466 (
467 convertor(originx + ngridx * spacingx - 0.5 * spacingx, "bohr", "Angstrom"),
468 convertor(originy + ngridy * spacingy - 0.5 * spacingy, "bohr", "Angstrom"),
469 convertor(originz + ngridz * spacingz - 0.5 * spacingz, "bohr", "Angstrom"),
470 ),
471 (
472 convertor(spacingx, "bohr", "Angstrom"),
473 convertor(spacingy, "bohr", "Angstrom"),
474 convertor(spacingz, "bohr", "Angstrom"),
475 ),
476 )
478 # Check grid size
479 assert vol.data.shape == (ngridx, ngridy, ngridz)
481 # Assign grid data
482 vol.data = numpy.resize(tmp, vol.data.shape)
484 return vol