Source code for pyscf.tools.cubegen

#!/usr/bin/env python
# Copyright 2014-2020 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#         Peter Koval <koval.peter@gmail.com>
#         Paul J. Robinson <pjrobinson@ucla.edu>
#

'''
Gaussian cube file format.  Reference:
http://paulbourke.net/dataformats/cube/
https://h5cube-spec.readthedocs.io/en/latest/cubeformat.html
http://gaussian.com/cubegen/

The output cube file has the following format

Comment line
Comment line
N_atom Ox Oy Oz         # number of atoms, followed by the coordinates of the origin
N1 vx1 vy1 vz1          # number of grids along each axis, followed by the step size in x/y/z direction.
N2 vx2 vy2 vz2          # ...
N3 vx3 vy3 vz3          # ...
Atom1 Z1 x y z          # Atomic number, charge, and coordinates of the atom
...                     # ...
AtomN ZN x y z          # ...
Data on grids           # (N1*N2) lines of records, each line has N3 elements
'''

import time
import numpy
import pyscf
from pyscf import lib
from pyscf import gto
from pyscf import df
from pyscf.dft import numint
from pyscf import __config__

RESOLUTION = getattr(__config__, 'cubegen_resolution', None)
BOX_MARGIN = getattr(__config__, 'cubegen_box_margin', 3.0)
ORIGIN = getattr(__config__, 'cubegen_box_origin', None)
# If given, EXTENT should be a 3-element ndarray/list/tuple to represent the
# extension in x, y, z
EXTENT = getattr(__config__, 'cubegen_box_extent', None)


[docs] def density(mol, outfile, dm, nx=80, ny=80, nz=80, resolution=RESOLUTION, margin=BOX_MARGIN): """Calculates electron density and write out in cube format. Args: mol : Mole Molecule to calculate the electron density for. outfile : str Name of Cube file to be written. dm : ndarray Density matrix of molecule. Kwargs: nx : int Number of grid point divisions in x direction. Note this is function of the molecule's size; a larger molecule will have a coarser representation than a smaller one for the same value. Conflicts to keyword resolution. ny : int Number of grid point divisions in y direction. nz : int Number of grid point divisions in z direction. resolution: float Resolution of the mesh grid in the cube box. If resolution is given in the input, the input nx/ny/nz have no effects. The value of nx/ny/nz will be determined by the resolution and the cube box size. """ from pyscf.pbc.gto import Cell cc = Cube(mol, nx, ny, nz, resolution, margin) GTOval = 'GTOval' if isinstance(mol, Cell): GTOval = 'PBC' + GTOval # Compute density on the .cube grid coords = cc.get_coords() ngrids = cc.get_ngrids() blksize = min(8000, ngrids) rho = numpy.empty(ngrids) for ip0, ip1 in lib.prange(0, ngrids, blksize): ao = mol.eval_gto(GTOval, coords[ip0:ip1]) rho[ip0:ip1] = numint.eval_rho(mol, ao, dm) rho = rho.reshape(cc.nx,cc.ny,cc.nz) # Write out density to the .cube file cc.write(rho, outfile, comment='Electron density in real space (e/Bohr^3)') return rho
[docs] def orbital(mol, outfile, coeff, nx=80, ny=80, nz=80, resolution=RESOLUTION, margin=BOX_MARGIN): """Calculate orbital value on real space grid and write out in cube format. Args: mol : Mole Molecule to calculate the electron density for. outfile : str Name of Cube file to be written. coeff : 1D array coeff coefficient. Kwargs: nx : int Number of grid point divisions in x direction. Note this is function of the molecule's size; a larger molecule will have a coarser representation than a smaller one for the same value. Conflicts to keyword resolution. ny : int Number of grid point divisions in y direction. nz : int Number of grid point divisions in z direction. resolution: float Resolution of the mesh grid in the cube box. If resolution is given in the input, the input nx/ny/nz have no effects. The value of nx/ny/nz will be determined by the resolution and the cube box size. """ from pyscf.pbc.gto import Cell cc = Cube(mol, nx, ny, nz, resolution, margin) GTOval = 'GTOval' if isinstance(mol, Cell): GTOval = 'PBC' + GTOval # Compute density on the .cube grid coords = cc.get_coords() ngrids = cc.get_ngrids() blksize = min(8000, ngrids) orb_on_grid = numpy.empty(ngrids) for ip0, ip1 in lib.prange(0, ngrids, blksize): ao = mol.eval_gto(GTOval, coords[ip0:ip1]) orb_on_grid[ip0:ip1] = numpy.dot(ao, coeff) orb_on_grid = orb_on_grid.reshape(cc.nx,cc.ny,cc.nz) # Write out orbital to the .cube file cc.write(orb_on_grid, outfile, comment='Orbital value in real space (1/Bohr^3)') return orb_on_grid
[docs] def mep(mol, outfile, dm, nx=80, ny=80, nz=80, resolution=RESOLUTION, margin=BOX_MARGIN): """Calculates the molecular electrostatic potential (MEP) and write out in cube format. Args: mol : Mole Molecule to calculate the electron density for. outfile : str Name of Cube file to be written. dm : ndarray Density matrix of molecule. Kwargs: nx : int Number of grid point divisions in x direction. Note this is function of the molecule's size; a larger molecule will have a coarser representation than a smaller one for the same value. Conflicts to keyword resolution. ny : int Number of grid point divisions in y direction. nz : int Number of grid point divisions in z direction. resolution: float Resolution of the mesh grid in the cube box. If resolution is given in the input, the input nx/ny/nz have no effects. The value of nx/ny/nz will be determined by the resolution and the cube box size. """ cc = Cube(mol, nx, ny, nz, resolution, margin) coords = cc.get_coords() # Nuclear potential at given points Vnuc = 0 for i in range(mol.natm): r = mol.atom_coord(i) Z = mol.atom_charge(i) rp = r - coords Vnuc += Z / numpy.einsum('xi,xi->x', rp, rp)**.5 # Potential of electron density Vele = numpy.empty_like(Vnuc) for p0, p1 in lib.prange(0, Vele.size, 600): fakemol = gto.fakemol_for_charges(coords[p0:p1]) ints = df.incore.aux_e2(mol, fakemol) Vele[p0:p1] = numpy.einsum('ijp,ij->p', ints, dm) MEP = Vnuc - Vele # MEP at each point MEP = MEP.reshape(cc.nx,cc.ny,cc.nz) # Write the potential cc.write(MEP, outfile, 'Molecular electrostatic potential in real space') return MEP
[docs] class Cube: ''' Read-write of the Gaussian CUBE files Attributes: nx : int Number of grid point divisions in x direction. Note this is function of the molecule's size; a larger molecule will have a coarser representation than a smaller one for the same value. Conflicts to keyword resolution. ny : int Number of grid point divisions in y direction. nz : int Number of grid point divisions in z direction. resolution: float Resolution of the mesh grid in the cube box. If resolution is given in the input, the input nx/ny/nz have no effects. The value of nx/ny/nz will be determined by the resolution and the cube box size. The unit is Bohr. ''' def __init__(self, mol, nx=80, ny=80, nz=80, resolution=RESOLUTION, margin=BOX_MARGIN, origin=ORIGIN, extent=EXTENT): from pyscf.pbc.gto import Cell self.mol = mol coord = mol.atom_coords() # If the molecule is periodic, use lattice vectors as the box # and automatically determine margin, origin, and extent if isinstance(mol, Cell): self.box = mol.lattice_vectors() atom_center = (numpy.max(coord, axis=0) + numpy.min(coord, axis=0))/2 box_center = (self.box[0] + self.box[1] + self.box[2])/2 self.boxorig = atom_center - box_center else: # Create a box based on its origin and extents/lengths (rectangular cuboid only) # If extent is not supplied, use the coordinates plus a margin on both sides if extent is None: extent = numpy.max(coord, axis=0) - numpy.min(coord, axis=0) + 2*margin self.box = numpy.diag(extent) # if origin is not supplied, set it as the minimum coordinate minus a margin. if origin is None: origin = numpy.min(coord, axis=0) - margin self.boxorig = numpy.asarray(origin) if resolution is not None: nx, ny, nz = numpy.ceil(numpy.diag(self.box) / resolution).astype(int) self.nx = nx self.ny = ny self.nz = nz if isinstance(mol, Cell): # Use an asymmetric mesh for tiling unit cells self.xs = numpy.linspace(0, 1, nx, endpoint=False) self.ys = numpy.linspace(0, 1, ny, endpoint=False) self.zs = numpy.linspace(0, 1, nz, endpoint=False) else: # Use endpoint=True to get a symmetric mesh # see also the discussion https://github.com/sunqm/pyscf/issues/154 self.xs = numpy.linspace(0, 1, nx, endpoint=True) self.ys = numpy.linspace(0, 1, ny, endpoint=True) self.zs = numpy.linspace(0, 1, nz, endpoint=True)
[docs] def get_coords(self) : """ Result: set of coordinates to compute a field which is to be stored in the file. """ frac_coords = lib.cartesian_prod([self.xs, self.ys, self.zs]) return frac_coords @ self.box + self.boxorig # Convert fractional coordinates to real-space coordinates
[docs] def get_ngrids(self): return self.nx * self.ny * self.nz
[docs] def get_volume_element(self): return (self.xs[1]-self.xs[0])*(self.ys[1]-self.ys[0])*(self.zs[1]-self.zs[0])
[docs] def write(self, field, fname, comment=None): """ Result: .cube file with the field in the file fname. """ assert (field.ndim == 3) assert (field.shape == (self.nx, self.ny, self.nz)) if comment is None: comment = 'Generic field? Supply the optional argument "comment" to define this line' mol = self.mol coord = mol.atom_coords() with open(fname, 'w') as f: f.write(comment+'\n') f.write(f'PySCF Version: {pyscf.__version__} Date: {time.ctime()}\n') f.write(f'{mol.natm:5d}') f.write('%12.6f%12.6f%12.6f\n' % tuple(self.boxorig.tolist())) dx = self.xs[-1] if len(self.xs) == 1 else self.xs[1] dy = self.ys[-1] if len(self.ys) == 1 else self.ys[1] dz = self.zs[-1] if len(self.zs) == 1 else self.zs[1] delta = (self.box.T * [dx,dy,dz]).T f.write(f'{self.nx:5d}{delta[0,0]:12.6f}{delta[0,1]:12.6f}{delta[0,2]:12.6f}\n') f.write(f'{self.ny:5d}{delta[1,0]:12.6f}{delta[1,1]:12.6f}{delta[1,2]:12.6f}\n') f.write(f'{self.nz:5d}{delta[2,0]:12.6f}{delta[2,1]:12.6f}{delta[2,2]:12.6f}\n') for ia in range(mol.natm): atmsymb = mol.atom_symbol(ia) f.write('%5d%12.6f'% (gto.charge(atmsymb), 0.)) f.write('%12.6f%12.6f%12.6f\n' % tuple(coord[ia])) for ix in range(self.nx): for iy in range(self.ny): for iz0, iz1 in lib.prange(0, self.nz, 6): fmt = '%13.5E' * (iz1-iz0) + '\n' f.write(fmt % tuple(field[ix,iy,iz0:iz1].tolist()))
[docs] def read(self, cube_file): with open(cube_file, 'r') as f: f.readline() f.readline() data = f.readline().split() natm = int(data[0]) self.boxorig = numpy.array([float(x) for x in data[1:]]) def parse_nx(data): from pyscf.pbc.gto import Cell d = data.split() nx = int(d[0]) x_vec = numpy.array([float(x) for x in d[1:]]) * nx if isinstance(self.mol, Cell): # Use an asymmetric mesh for tiling unit cells xs = numpy.linspace(0, 1, nx, endpoint=False) else: # Use endpoint=True to get a symmetric mesh # see also the discussion https://github.com/sunqm/pyscf/issues/154 xs = numpy.linspace(0, 1, nx, endpoint=True) return x_vec, nx, xs self.box = numpy.zeros((3,3)) self.box[0], self.nx, self.xs = parse_nx(f.readline()) self.box[1], self.ny, self.ys = parse_nx(f.readline()) self.box[2], self.nz, self.zs = parse_nx(f.readline()) atoms = [] for ia in range(natm): d = f.readline().split() atoms.append([int(d[0]), [float(x) for x in d[2:]]]) self.mol = gto.M(atom=atoms, unit='Bohr') data = f.read() cube_data = numpy.array([float(x) for x in data.split()]) return cube_data.reshape([self.nx, self.ny, self.nz])
if __name__ == '__main__': from pyscf import scf from pyscf.tools import cubegen mol = gto.M(atom='''O 0.00000000, 0.000000, 0.000000 H 0.761561, 0.478993, 0.00000000 H -0.761561, 0.478993, 0.00000000''', basis='6-31g*') mf = scf.RHF(mol).run() cubegen.density(mol, 'h2o_den.cube', mf.make_rdm1()) #makes total density cubegen.mep(mol, 'h2o_pot.cube', mf.make_rdm1()) cubegen.orbital(mol, 'h2o_mo1.cube', mf.mo_coeff[:,0])