From 64c724372879bcd6bcd8df25a5c4cd53c6eb8492 Mon Sep 17 00:00:00 2001 From: Ivo Filot Date: Tue, 7 Nov 2023 11:55:52 +0100 Subject: [PATCH] Develop (#23) * Initial commit docs * Adding deployment docs * Adding deployment docs * Expanding documentation * Building docker image for docs * Adding Dockerfile * Adding example script * Expanding documentation * Expanding documentation * Redesigning README.md * Expanding documentation * Initial commit Foster-Boys procedure * Fixing requirements * Adding scipy in build and host * Adding scipy to environment.yml * Adding Foster-Boys script * Adding testing * Initial commit COHP and adding RNG seed for FB localization * Adding seed * Add multiple runners * Decreasing accuracy to two decimals for FB * Putting RNG as class member * Fixing nuclear derivatives * Fixing meta.yaml * Adding derivatives tests * Adding error catching in diis * Adding LiH example for COHP * Adding molecule builder * Adding molecules * Fixing manifest file * Fixing type * Making test more lenient * Update cohp_lih.py * Adding bounds for Foster-Boys * Refactoring HF derivatives * Pre-seeding calculation * Initial commit steepest descent * Initial commit geometric derivatives parallel parser * Update cohp_lih.py * Initial commit OpenMP calculation of derivatives * Adding steepest descent * More stringent conversion criterium * Reducing number of digits for MacOS * Adding conjugate gradient algorithm * Making test more lenient * Fixing issue #17 * Tracking history in geometry optimization class * Fixing bugs in HF algo * Fixing bugs in GeometryOptimization and in MoleculeBuilder * Fixing molecule builder test * Forcing CH4 * Adding .xyz file to unit test * Adding example scripts * Adding nelec() method * Adding charges * Adding molecule rendering * Incrementing version number * Adding safeguard for PyTessel * Removing tqdm and adding mendeleev in gaurds * Changing versions of dependencies * Adding tqdm and mendeleev as dependencies * Removing mendeleev as package * Fixing test environment * Removing potential troublesome mendeleev package * Testing compilation Python 3.11 * Specifying versions in conda_build_config.yml --- COMPILATION.md | 4 +- conda_build_config.yaml | 9 + environment.yml | 1 + examples/foster_boys.py | 81 ++++--- examples/molecular_orbitals_co.py | 148 ++++++------- examples/molecular_orbitals_hco.py | 74 +++++++ meta.yaml | 3 +- pyqint/__init__.py | 1 + pyqint/_version.py | 2 +- pyqint/blender/blender_render_molecule.py | 245 +++++++++++++++++++++ pyqint/blender/blender_render_molecules.py | 245 +++++++++++++++++++++ pyqint/blenderrender.py | 162 ++++++++++++++ pyqint/geometry_optimization.py | 1 + pyqint/hf.py | 8 +- pyqint/molecule.py | 18 ++ tests/test_hf_molecules_charged.py | 30 +++ tests/test_molecule_builder.py | 6 + 17 files changed, 918 insertions(+), 120 deletions(-) create mode 100644 examples/molecular_orbitals_hco.py create mode 100644 pyqint/blender/blender_render_molecule.py create mode 100644 pyqint/blender/blender_render_molecules.py create mode 100644 pyqint/blenderrender.py create mode 100644 tests/test_hf_molecules_charged.py diff --git a/COMPILATION.md b/COMPILATION.md index 7351055..501782b 100644 --- a/COMPILATION.md +++ b/COMPILATION.md @@ -99,11 +99,11 @@ python3 setup.py build and install it locally ``` -pip install -e . +pip3 install -e . ``` and finally test it ``` -pytest tests/* +pytest-3 tests/* ``` diff --git a/conda_build_config.yaml b/conda_build_config.yaml index da95e80..6034800 100644 --- a/conda_build_config.yaml +++ b/conda_build_config.yaml @@ -2,3 +2,12 @@ python: - 3.8 - 3.9 - 3.10 + - 3.11 +numpy: + - 1.22 + - 1.22 + - 1.22 + - 1.24 +zip_keys: + - python + - numpy diff --git a/environment.yml b/environment.yml index 131cd06..6779298 100644 --- a/environment.yml +++ b/environment.yml @@ -5,3 +5,4 @@ dependencies: - conda-build - conda-verify - anaconda-client + - tqdm diff --git a/examples/foster_boys.py b/examples/foster_boys.py index b0d2849..513d3a9 100644 --- a/examples/foster_boys.py +++ b/examples/foster_boys.py @@ -1,50 +1,69 @@ -from pyqint import Molecule, HF, PyQInt, FosterBoys -import pyqint +from pyqint import Molecule, PyQInt, FosterBoys, GeometryOptimization import numpy as np -from pytessel import PyTessel +import matplotlib.pyplot as plt # # Plot the isosurfaces for the CO molecule # def main(): - res = calculate_co(1.145414) + res = optimize_co() resfb = FosterBoys(res).run() + energies = resfb['orbe'] + coeff = resfb['orbc'] + cgfs = res['cgfs'] - for i in range(len(res['cgfs'])): - build_isosurface('MO_%03i' % (i+1), - res['cgfs'], - resfb['orbc'][:,i], - 0.1) + fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5)) + for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)): + res, x, y = build_contourplot(cgfs, chi, sz=3, plane='xz') + vmax = np.max(np.abs(res)) + vmin = -vmax + ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG', + vmin=vmin, vmax=vmax) + ax[i//5,i%5].contour(x, y, res, levels=15, colors='black', + vmin=vmin, vmax=vmax) + ax[i//5,i%5].set_xlabel('x [Bohr]') + ax[i//5,i%5].set_ylabel('z [Bohr]') + ax[i//5,i%5].set_title('Energy = %.4f Ht' % e) + ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5) + + plt.tight_layout() + plt.savefig('co_fb.jpg') -def calculate_co(d): +def optimize_co(): """ - Full function for evaluation + Optimization function for scipy.optimize.minimize """ mol = Molecule() - mol.add_atom('C', 0.0, 0.0, -d/2, unit='angstrom') - mol.add_atom('O', 0.0, 0.0, d/2, unit='angstrom') - - result = HF().rhf(mol, 'sto3g') - - return result + mol.add_atom('C', 0.0, 0.0, -0.6, unit='angstrom') + mol.add_atom('O', 0.0, 0.0, 0.6, unit='angstrom') + + res = GeometryOptimization().run(mol, 'sto3g') + + return res['data'] -def build_isosurface(filename, cgfs, coeff, isovalue, sz=5, npts=100): - # generate some data - isovalue = np.abs(isovalue) +def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'): integrator = PyQInt() - grid = integrator.build_rectgrid3d(-sz, sz, npts) - scalarfield = np.reshape(integrator.plot_wavefunction(grid, coeff, cgfs), (npts, npts, npts)) - unitcell = np.diag(np.ones(3) * 2 * sz) - - pytessel = PyTessel() - vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), isovalue) - fname = filename + '_pos.ply' - pytessel.write_ply(fname, vertices, normals, indices) - vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), -isovalue) - fname = filename + '_neg.ply' - pytessel.write_ply(fname, vertices, normals, indices) + # build grid + x = np.linspace(-sz, sz, 50) + y = np.linspace(-sz, sz, 50) + xx, yy = np.meshgrid(x,y) + zz = np.zeros(len(x) * len(y)) + + if plane == 'xy': + points = [xx.flatten(), yy.flatten(), zz] + elif plane == 'xz': + points = [xx.flatten(), zz, yy.flatten()] + elif plane == 'yz': + points = [zz, xx.flatten(), yy.flatten()] + else: + raise Exception('Unknown plane: %s' % plane) + + grid = np.vstack(points).reshape(3,-1).T + res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x))) + + return res, x, y if __name__ == '__main__': main() diff --git a/examples/molecular_orbitals_co.py b/examples/molecular_orbitals_co.py index f9bea40..3252e1e 100644 --- a/examples/molecular_orbitals_co.py +++ b/examples/molecular_orbitals_co.py @@ -1,80 +1,68 @@ -from pyqint import Molecule, HF, PyQInt -import numpy as np -import matplotlib.pyplot as plt - -# -# Plot the isosurfaces for the CO molecule -# - -def main(): - # calculate sto-3g coefficients for co - result = calculate_co(1.145414) - energies = result['orbe'] - coeff = result['orbc'] - - fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5)) - for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)): - res, x, y = build_contourplot(result['cgfs'], chi, sz=3, plane='xz') - vmax = np.max(np.abs(res)) - vmin = -vmax - ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG', - vmin=vmin, vmax=vmax) - ax[i//5,i%5].contour(x, y, res, levels=15, colors='black', - vmin=vmin, vmax=vmax) - ax[i//5,i%5].set_xlabel('x [Bohr]') - ax[i//5,i%5].set_ylabel('z [Bohr]') - ax[i//5,i%5].set_title('Energy = %.4f Ht' % e) - ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5) - - plt.tight_layout() - plt.savefig('co.jpg') - -def optimize_co(d): - """ - Optimization function for scipy.optimize.minimize - """ - mol = Molecule() - mol.add_atom('C', 0.0, 0.0, -d[0]/2, unit='angstrom') - mol.add_atom('O', 0.0, 0.0, d[0]/2, unit='angstrom') - - result = HF().rhf(mol, 'sto3g') - - return result['energy'] - -def calculate_co(d): - """ - Full function for evaluation - """ - mol = Molecule() - mol.add_atom('C', 0.0, 0.0, -d/2, unit='angstrom') - mol.add_atom('O', 0.0, 0.0, d/2, unit='angstrom') - - result = HF().rhf(mol, 'sto3g') - - return result - -def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'): - integrator = PyQInt() - - # build grid - x = np.linspace(-sz, sz, 50) - y = np.linspace(-sz, sz, 50) - xx, yy = np.meshgrid(x,y) - zz = np.zeros(len(x) * len(y)) - - if plane == 'xy': - points = [xx.flatten(), yy.flatten(), zz] - elif plane == 'xz': - points = [xx.flatten(), zz, yy.flatten()] - elif plane == 'yz': - points = [zz, xx.flatten(), yy.flatten()] - else: - raise Exception('Unknown plane: %s' % plane) - - grid = np.vstack(points).reshape(3,-1).T - res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x))) - - return res, x, y - -if __name__ == '__main__': - main() \ No newline at end of file +from pyqint import Molecule, HF, PyQInt, GeometryOptimization +import numpy as np +import matplotlib.pyplot as plt + +# +# Plot the isosurfaces for the CO molecule +# + +def main(): + # calculate sto-3g coefficients for co + result = optimize_co() + energies = result['orbe'] + coeff = result['orbc'] + + fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5)) + for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)): + res, x, y = build_contourplot(result['cgfs'], chi, sz=3, plane='xz') + vmax = np.max(np.abs(res)) + vmin = -vmax + ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG', + vmin=vmin, vmax=vmax) + ax[i//5,i%5].contour(x, y, res, levels=15, colors='black', + vmin=vmin, vmax=vmax) + ax[i//5,i%5].set_xlabel('x [Bohr]') + ax[i//5,i%5].set_ylabel('z [Bohr]') + ax[i//5,i%5].set_title('Energy = %.4f Ht' % e) + ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5) + + plt.tight_layout() + plt.savefig('co.jpg') + +def optimize_co(): + """ + Optimization function for scipy.optimize.minimize + """ + mol = Molecule() + mol.add_atom('C', 0.0, 0.0, -0.6, unit='angstrom') + mol.add_atom('O', 0.0, 0.0, 0.6, unit='angstrom') + + res = GeometryOptimization().run(mol, 'sto3g') + + return res['data'] + +def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'): + integrator = PyQInt() + + # build grid + x = np.linspace(-sz, sz, 50) + y = np.linspace(-sz, sz, 50) + xx, yy = np.meshgrid(x,y) + zz = np.zeros(len(x) * len(y)) + + if plane == 'xy': + points = [xx.flatten(), yy.flatten(), zz] + elif plane == 'xz': + points = [xx.flatten(), zz, yy.flatten()] + elif plane == 'yz': + points = [zz, xx.flatten(), yy.flatten()] + else: + raise Exception('Unknown plane: %s' % plane) + + grid = np.vstack(points).reshape(3,-1).T + res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x))) + + return res, x, y + +if __name__ == '__main__': + main() diff --git a/examples/molecular_orbitals_hco.py b/examples/molecular_orbitals_hco.py new file mode 100644 index 0000000..e3b1f44 --- /dev/null +++ b/examples/molecular_orbitals_hco.py @@ -0,0 +1,74 @@ +from pyqint import Molecule, HF, PyQInt, GeometryOptimization +import numpy as np +import matplotlib.pyplot as plt + +# +# Plot the isosurfaces for the CO molecule +# + +def main(): + # calculate sto-3g coefficients for co + result = optimize_hco() + energies = result['orbe'] + coeff = result['orbc'] + + fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5)) + for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)): + + if i >= 10: + break + + res, x, y = build_contourplot(result['cgfs'], chi, sz=3, plane='xz') + vmax = np.max(np.abs(res)) + vmin = -vmax + ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG', + vmin=vmin, vmax=vmax) + ax[i//5,i%5].contour(x, y, res, levels=15, colors='black', + vmin=vmin, vmax=vmax) + ax[i//5,i%5].set_xlabel('x [Bohr]') + ax[i//5,i%5].set_ylabel('z [Bohr]') + ax[i//5,i%5].set_title('Energy = %.4f Ht' % e) + ax[i//5,i%5].grid(linestyle='--', color='black', alpha=0.5) + + plt.tight_layout() + plt.savefig('hco.jpg') + +def optimize_hco(): + """ + Optimization function for scipy.optimize.minimize + """ + mol = Molecule() + mol.add_atom('O', -0.12770367, -0.00000000, 1.23419284) + mol.add_atom('C', -0.50625665, 0.00000000, -1.09149431) + mol.add_atom('H', 1.57882331, -0.00000000, -2.06681794) + mol.set_charge(-1) + + res = GeometryOptimization().run(mol, 'p321') + + return res['data'] + +def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'): + integrator = PyQInt() + + # build grid + x = np.linspace(-sz, sz, 50) + y = np.linspace(-sz, sz, 50) + xx, yy = np.meshgrid(x,y) + zz = np.zeros(len(x) * len(y)) + + if plane == 'xy': + points = [xx.flatten(), yy.flatten(), zz] + elif plane == 'xz': + points = [xx.flatten(), zz, yy.flatten()] + elif plane == 'yz': + points = [zz, xx.flatten(), yy.flatten()] + else: + raise Exception('Unknown plane: %s' % plane) + + grid = np.vstack(points).reshape(3,-1).T + res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x))) + + return res, x, y + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/meta.yaml b/meta.yaml index 0b47242..8fbc2da 100644 --- a/meta.yaml +++ b/meta.yaml @@ -1,6 +1,6 @@ package: name: "pyqint" - version: "0.12.2.0" + version: "0.13.0.0" source: path: . @@ -23,6 +23,7 @@ requirements: - python - numpy - scipy + - tqdm test: requires: diff --git a/pyqint/__init__.py b/pyqint/__init__.py index c201978..4e1a87b 100644 --- a/pyqint/__init__.py +++ b/pyqint/__init__.py @@ -7,5 +7,6 @@ from .cohp import COHP from .molecule_builder import MoleculeBuilder from .geometry_optimization import GeometryOptimization +from .blenderrender import BlenderRender from ._version import __version__ diff --git a/pyqint/_version.py b/pyqint/_version.py index d19ae08..be00ee7 100644 --- a/pyqint/_version.py +++ b/pyqint/_version.py @@ -1,2 +1,2 @@ -__version__ = "0.12.2.0" +__version__ = "0.14.0.0" diff --git a/pyqint/blender/blender_render_molecule.py b/pyqint/blender/blender_render_molecule.py new file mode 100644 index 0000000..c986444 --- /dev/null +++ b/pyqint/blender/blender_render_molecule.py @@ -0,0 +1,245 @@ +import bpy +import numpy as np +import os +import time +import json + +# +# IMPORTANT +# +# Do not run this script natively. This script is meant to be run in Blender +# via one of the call routines +# + +with open(os.path.join(os.path.dirname(__file__), 'manifest.json')) as f: + manifest = json.load(f) + +def main(): + # set the scene + settings = { + 'resolution': 512, + 'camera_location': (-10,0,0), + 'camera_rotation': (np.pi/2,0,-np.pi/2), + 'camera_scale' : 10 + } + set_environment(settings) + + # read molecule file and load it + mol = read_xyz(manifest['xyzfile']) + create_atoms(mol) + create_bonds(mol) + + add_isosurface(manifest['mo_name'], + manifest['mo_neg_path'], + manifest['mo_pos_path']) + + render_scene(manifest['png_output']) + +def add_isosurface(label, filename_neg, filename_pos): + """ + Add the two isosurfaces + """ + bpy.ops.import_mesh.ply( + filepath=filename_neg + ) + obj = bpy.context.object + bpy.ops.object.shade_smooth() + obj.data.materials.append(create_material('matneg', manifest['mo_colors']['neg'], alpha=0.5)) + obj.name = 'isosurface ' + label + '_neg' + + obj = bpy.ops.import_mesh.ply( + filepath=filename_pos + ) + bpy.ops.object.shade_smooth() + obj = bpy.context.object + obj.data.materials.append(create_material('matpos', manifest['mo_colors']['pos'], alpha=0.5)) + obj.name = 'isosurface ' + label + '_pos' + +def create_atoms(mol): + """ + Create atoms + """ + for i,at in enumerate(mol): + scale = manifest['atom_radii'][at[0]] + bpy.ops.surface.primitive_nurbs_surface_sphere_add( + radius=scale, + enter_editmode=False, + align='WORLD', + location=at[1]) + obj = bpy.context.view_layer.objects.active + obj.name = "atom-%s-%03i" % (at[0],i) + bpy.ops.object.shade_smooth() + + # set a material + mat = create_material(at[0], manifest['atom_colors'][at[0]]) + print(mat) + obj.data.materials.append(mat) + +def create_bonds(mol): + """ + Create bonds between atoms + """ + # set default orientation of bonds (fixed!) + z = np.array([0,0,1]) + + # add new bonds material if it does not yet exist + matbond = create_material('bond', '222222') + + for i,at1 in enumerate(mol): + r1 = np.array(at1[1]) + for j,at2 in enumerate(mol[i+1:]): + r2 = np.array(at2[1]) + dist = np.linalg.norm(r2 - r1) + + # only create a bond if the distance is less than 1.5 A + if dist < 2.5: + axis = np.cross(z,r2-r1) + if np.linalg.norm(axis) < 1e-5: + axis = np.array([0,0,1]) + angle = 0.0 + else: + axis /= np.linalg.norm(axis) + angle = np.arccos(np.dot(r2-r1,z)/dist) + + bpy.ops.surface.primitive_nurbs_surface_cylinder_add( + enter_editmode=False, + align='WORLD', + location=tuple((r1 + r2) / 2) + ) + + obj = bpy.context.view_layer.objects.active + obj.scale = (manifest['bond_thickness'], manifest['bond_thickness'], dist/2) + obj.rotation_mode = 'AXIS_ANGLE' + obj.rotation_axis_angle = (angle, axis[0], axis[1], axis[2]) + + obj.name = "bond-%s-%03i-%s-%03i" % (at1[0],i,at2[0],j) + bpy.ops.object.shade_smooth() + obj.data.materials.append(matbond) + +def set_environment(settings): + """ + Specify canvas size, remove default objects, reset positions of + camera and light, define film and set background + """ + bpy.context.scene.render.engine = 'CYCLES' + bpy.context.scene.cycles.device = 'GPU' + bpy.context.scene.render.resolution_x = settings['resolution'] + bpy.context.scene.render.resolution_y = settings['resolution'] + bpy.context.scene.cycles.samples = 1024 + bpy.context.scene.cycles.tile_size = 2048 + + # remove cube + if 'Cube' in bpy.data.objects: + o = bpy.data.objects['Cube'] + bpy.data.objects.remove(o, do_unlink=True) + + # set camera into default position + bpy.data.objects['Camera'].location = tuple(settings['camera_location']) + bpy.data.objects['Camera'].rotation_euler = tuple(settings['camera_rotation']) + bpy.data.objects['Camera'].data.clip_end = 1000 + bpy.data.objects['Camera'].data.type = 'ORTHO' + bpy.data.objects['Camera'].data.ortho_scale = settings['camera_scale'] + + # set lights + bpy.data.objects['Light'].data.type = 'AREA' + bpy.data.objects['Light'].data.energy = 1e4 + bpy.data.objects['Light'].location = (-10,10,10) + bpy.data.objects['Light'].rotation_euler = tuple(np.radians([55, 0, 225])) + bpy.data.objects['Light'].data.shape = 'DISK' + bpy.data.objects['Light'].data.size = 10 + + # set film + bpy.context.scene.render.film_transparent = True + + # set background + bpy.data.worlds["World"].node_tree.nodes["Background"].inputs[0].default_value = (1,1,1,1) + bpy.data.worlds["World"].node_tree.nodes["Background"].inputs[1].default_value = 1 + +def create_material(name, color, alpha=1.0): + """ + Build a new material + """ + # early exit if material already exists + if name in bpy.data.materials: + return bpy.data.materials[name] + + mat = bpy.data.materials.new(name) + mat.use_nodes = True + + # set base color + mat.node_tree.nodes["Principled BSDF"].inputs[0].default_value = hex2rgbtuple(color) + + # subsurface modifier + mat.node_tree.nodes["Principled BSDF"].inputs[1].default_value = 0.2 + + # set subsurface radii + mat.node_tree.nodes["Principled BSDF"].inputs[2].default_value = (0.3,0.3,0.3) + + # set subsurface color + mat.node_tree.nodes["Principled BSDF"].inputs[3].default_value = hex2rgbtuple('000000') + + # metallic + mat.node_tree.nodes["Principled BSDF"].inputs[4].default_value = 0.3 + + # roughness + mat.node_tree.nodes["Principled BSDF"].inputs[7].default_value = 0.05 + + # alpha + mat.node_tree.nodes["Principled BSDF"].inputs[21].default_value = alpha + + return mat + +def render_scene(outputfile, samples=512): + """ + Render the scene + """ + bpy.context.scene.cycles.samples = samples + + print('Start render') + start = time.time() + bpy.data.scenes['Scene'].render.filepath = outputfile + bpy.ops.render.render(write_still=True) + end = time.time() + print('Finished rendering frame in %.1f seconds' % (end - start)) + +def read_xyz(filename): + f = open(filename) + nratoms = int(f.readline()) + f.readline() # skip line + + angtobohr = 1.8897259886 + + mol = [] + for i in range(0,nratoms): + pieces = f.readline().split() + mol.append( + [pieces[0], + ( + float(pieces[1]) * angtobohr, + float(pieces[2]) * angtobohr, + float(pieces[3]) * angtobohr + )] + ) + + return mol + +def hex2rgbtuple(hexcode): + """ + Convert 6-digit color hexcode to a tuple of floats + """ + hexcode += "FF" + hextuple = tuple([int(hexcode[i:i+2], 16)/255.0 for i in [0,2,4,6]]) + + return tuple([color_srgb_to_scene_linear(c) for c in hextuple]) + +def color_srgb_to_scene_linear(c): + """ + Convert RGB to sRGB + """ + if c < 0.04045: + return 0.0 if c < 0.0 else c * (1.0 / 12.92) + else: + return ((c + 0.055) * (1.0 / 1.055)) ** 2.4 + +if __name__ == '__main__': + main() diff --git a/pyqint/blender/blender_render_molecules.py b/pyqint/blender/blender_render_molecules.py new file mode 100644 index 0000000..c986444 --- /dev/null +++ b/pyqint/blender/blender_render_molecules.py @@ -0,0 +1,245 @@ +import bpy +import numpy as np +import os +import time +import json + +# +# IMPORTANT +# +# Do not run this script natively. This script is meant to be run in Blender +# via one of the call routines +# + +with open(os.path.join(os.path.dirname(__file__), 'manifest.json')) as f: + manifest = json.load(f) + +def main(): + # set the scene + settings = { + 'resolution': 512, + 'camera_location': (-10,0,0), + 'camera_rotation': (np.pi/2,0,-np.pi/2), + 'camera_scale' : 10 + } + set_environment(settings) + + # read molecule file and load it + mol = read_xyz(manifest['xyzfile']) + create_atoms(mol) + create_bonds(mol) + + add_isosurface(manifest['mo_name'], + manifest['mo_neg_path'], + manifest['mo_pos_path']) + + render_scene(manifest['png_output']) + +def add_isosurface(label, filename_neg, filename_pos): + """ + Add the two isosurfaces + """ + bpy.ops.import_mesh.ply( + filepath=filename_neg + ) + obj = bpy.context.object + bpy.ops.object.shade_smooth() + obj.data.materials.append(create_material('matneg', manifest['mo_colors']['neg'], alpha=0.5)) + obj.name = 'isosurface ' + label + '_neg' + + obj = bpy.ops.import_mesh.ply( + filepath=filename_pos + ) + bpy.ops.object.shade_smooth() + obj = bpy.context.object + obj.data.materials.append(create_material('matpos', manifest['mo_colors']['pos'], alpha=0.5)) + obj.name = 'isosurface ' + label + '_pos' + +def create_atoms(mol): + """ + Create atoms + """ + for i,at in enumerate(mol): + scale = manifest['atom_radii'][at[0]] + bpy.ops.surface.primitive_nurbs_surface_sphere_add( + radius=scale, + enter_editmode=False, + align='WORLD', + location=at[1]) + obj = bpy.context.view_layer.objects.active + obj.name = "atom-%s-%03i" % (at[0],i) + bpy.ops.object.shade_smooth() + + # set a material + mat = create_material(at[0], manifest['atom_colors'][at[0]]) + print(mat) + obj.data.materials.append(mat) + +def create_bonds(mol): + """ + Create bonds between atoms + """ + # set default orientation of bonds (fixed!) + z = np.array([0,0,1]) + + # add new bonds material if it does not yet exist + matbond = create_material('bond', '222222') + + for i,at1 in enumerate(mol): + r1 = np.array(at1[1]) + for j,at2 in enumerate(mol[i+1:]): + r2 = np.array(at2[1]) + dist = np.linalg.norm(r2 - r1) + + # only create a bond if the distance is less than 1.5 A + if dist < 2.5: + axis = np.cross(z,r2-r1) + if np.linalg.norm(axis) < 1e-5: + axis = np.array([0,0,1]) + angle = 0.0 + else: + axis /= np.linalg.norm(axis) + angle = np.arccos(np.dot(r2-r1,z)/dist) + + bpy.ops.surface.primitive_nurbs_surface_cylinder_add( + enter_editmode=False, + align='WORLD', + location=tuple((r1 + r2) / 2) + ) + + obj = bpy.context.view_layer.objects.active + obj.scale = (manifest['bond_thickness'], manifest['bond_thickness'], dist/2) + obj.rotation_mode = 'AXIS_ANGLE' + obj.rotation_axis_angle = (angle, axis[0], axis[1], axis[2]) + + obj.name = "bond-%s-%03i-%s-%03i" % (at1[0],i,at2[0],j) + bpy.ops.object.shade_smooth() + obj.data.materials.append(matbond) + +def set_environment(settings): + """ + Specify canvas size, remove default objects, reset positions of + camera and light, define film and set background + """ + bpy.context.scene.render.engine = 'CYCLES' + bpy.context.scene.cycles.device = 'GPU' + bpy.context.scene.render.resolution_x = settings['resolution'] + bpy.context.scene.render.resolution_y = settings['resolution'] + bpy.context.scene.cycles.samples = 1024 + bpy.context.scene.cycles.tile_size = 2048 + + # remove cube + if 'Cube' in bpy.data.objects: + o = bpy.data.objects['Cube'] + bpy.data.objects.remove(o, do_unlink=True) + + # set camera into default position + bpy.data.objects['Camera'].location = tuple(settings['camera_location']) + bpy.data.objects['Camera'].rotation_euler = tuple(settings['camera_rotation']) + bpy.data.objects['Camera'].data.clip_end = 1000 + bpy.data.objects['Camera'].data.type = 'ORTHO' + bpy.data.objects['Camera'].data.ortho_scale = settings['camera_scale'] + + # set lights + bpy.data.objects['Light'].data.type = 'AREA' + bpy.data.objects['Light'].data.energy = 1e4 + bpy.data.objects['Light'].location = (-10,10,10) + bpy.data.objects['Light'].rotation_euler = tuple(np.radians([55, 0, 225])) + bpy.data.objects['Light'].data.shape = 'DISK' + bpy.data.objects['Light'].data.size = 10 + + # set film + bpy.context.scene.render.film_transparent = True + + # set background + bpy.data.worlds["World"].node_tree.nodes["Background"].inputs[0].default_value = (1,1,1,1) + bpy.data.worlds["World"].node_tree.nodes["Background"].inputs[1].default_value = 1 + +def create_material(name, color, alpha=1.0): + """ + Build a new material + """ + # early exit if material already exists + if name in bpy.data.materials: + return bpy.data.materials[name] + + mat = bpy.data.materials.new(name) + mat.use_nodes = True + + # set base color + mat.node_tree.nodes["Principled BSDF"].inputs[0].default_value = hex2rgbtuple(color) + + # subsurface modifier + mat.node_tree.nodes["Principled BSDF"].inputs[1].default_value = 0.2 + + # set subsurface radii + mat.node_tree.nodes["Principled BSDF"].inputs[2].default_value = (0.3,0.3,0.3) + + # set subsurface color + mat.node_tree.nodes["Principled BSDF"].inputs[3].default_value = hex2rgbtuple('000000') + + # metallic + mat.node_tree.nodes["Principled BSDF"].inputs[4].default_value = 0.3 + + # roughness + mat.node_tree.nodes["Principled BSDF"].inputs[7].default_value = 0.05 + + # alpha + mat.node_tree.nodes["Principled BSDF"].inputs[21].default_value = alpha + + return mat + +def render_scene(outputfile, samples=512): + """ + Render the scene + """ + bpy.context.scene.cycles.samples = samples + + print('Start render') + start = time.time() + bpy.data.scenes['Scene'].render.filepath = outputfile + bpy.ops.render.render(write_still=True) + end = time.time() + print('Finished rendering frame in %.1f seconds' % (end - start)) + +def read_xyz(filename): + f = open(filename) + nratoms = int(f.readline()) + f.readline() # skip line + + angtobohr = 1.8897259886 + + mol = [] + for i in range(0,nratoms): + pieces = f.readline().split() + mol.append( + [pieces[0], + ( + float(pieces[1]) * angtobohr, + float(pieces[2]) * angtobohr, + float(pieces[3]) * angtobohr + )] + ) + + return mol + +def hex2rgbtuple(hexcode): + """ + Convert 6-digit color hexcode to a tuple of floats + """ + hexcode += "FF" + hextuple = tuple([int(hexcode[i:i+2], 16)/255.0 for i in [0,2,4,6]]) + + return tuple([color_srgb_to_scene_linear(c) for c in hextuple]) + +def color_srgb_to_scene_linear(c): + """ + Convert RGB to sRGB + """ + if c < 0.04045: + return 0.0 if c < 0.0 else c * (1.0 / 12.92) + else: + return ((c + 0.055) * (1.0 / 1.055)) ** 2.4 + +if __name__ == '__main__': + main() diff --git a/pyqint/blenderrender.py b/pyqint/blenderrender.py new file mode 100644 index 0000000..1a2fecb --- /dev/null +++ b/pyqint/blenderrender.py @@ -0,0 +1,162 @@ +import os +import tempfile +from pyqint import PyQInt +import numpy as np +import json +import subprocess +import shutil +import tqdm + +# try to import PyTessel but do not throw an error if it cannot be loaded +try: + from pytessel import PyTessel +except ModuleNotFoundError: + print('Cannot find PyTessel') + +try: + from mendeleev import element +except ModuleNotFoundError: + print('Cannot find mendeleev') + +class BlenderRender: + """ + This class leverages blender for rendering molecular orbitals + """ + def __init__(self): + self.executable = self.__find_blender() + if self.executable is None: + raise Exception('Cannot find Blender executable') + else: + print('Found executable: %s' % self.executable) + + def render_molecular_orbitals(self, molecule, cgfs, orbc, outpath, + mo_indices=None, sz=5, isovalue=0.03, + npts=100): + if mo_indices is None: # render all orbitals + mo_indices = np.arange(0, len(orbc)) + + # build a temporary folder + tempdir = tempfile.mkdtemp() + # store .xyz file + xyzfile = os.path.join(tempdir, 'mol.xyz') + self.__store_xyz(molecule, xyzfile) + + for idx in mo_indices: + # build isosurfaces + plyfile = os.path.join(tempdir, 'MO_%04i' % (idx+1)) + plypos, plyneg = self.__build_isosurface(plyfile, cgfs, orbc[:,idx], isovalue, sz, npts) + + # execute blender + outfile = os.path.join(outpath, 'MO_%04i.png' % (idx+1)) + self.__run_blender(plypos, plyneg, xyzfile, outfile, tempdir) + + shutil.rmtree(tempdir) + + def get_executable(self): + """ + Get the Blender executable + """ + return self.executable + + def __find_blender(self): + """ + Find the Blender executable + """ + searchpath = os.path.join('C:\\','Program Files','Blender Foundation') + name = 'blender.exe' + results = [] + for root, dirs, files in os.walk(searchpath): + if name in files: + results.append(os.path.join(root, name)) + + for res in results: + if '3.3' in res: + return res + + return None + + def __run_blender(self, negfile, posfile, xyzfile, pngfile, cwd): + # set path to xyz file + blendpysrc = os.path.join(os.path.dirname(__file__), 'blender', 'blender_render_molecule.py') + blendpydst = os.path.join(cwd, 'blender_render_molecule.py') + shutil.copyfile(blendpysrc, blendpydst) + + manifest = { + 'mo_colors' : { + 'neg': 'E72F65', + 'pos': '3F9EE7', + }, + 'xyzfile' : xyzfile, + 'mo_name' : 'isosurface', + 'mo_neg_path' : negfile, + 'mo_pos_path' : posfile, + 'png_output': pngfile, + 'bond_thickness': 0.2, + 'atom_radii' : { + 'H': 0.4, + 'N': 0.6, + 'C': 0.6, + 'O': 0.6, + }, + 'atom_colors' : { + 'H': 'FFFFFF', + 'N': '0000FF', + 'C': '000000', + 'O': 'DD0000', + } + } + + with open(os.path.join(cwd, 'manifest.json'), 'w') as f: + f.write(json.dumps(manifest)) + + # run blender + out = subprocess.check_output( + [self.executable, '-b', '-P', blendpydst], + cwd=cwd + ) + + def __store_xyz(self, mol, filename): + """ + Convert final result of GeometryOptimization class to an .xyz file + which can be relayed to Blender for rendering + + Note that we have to convert Bohr units back to Angstrom as .xyz files + are in units of Angstrom. Blender actually converts them from Angstrom back + to Bohr, which looks like wasted effort but this way all .xyz files + are consistent. + """ + f = open(filename, 'w') + f.write(str(len(mol.atoms)) + '\n') + f.write('\n') + + angtobohr = 1.8897259886 + + for a in mol.atoms: + elname = element(a[0]).symbol + f.write('%s %12.6f %12.6f %12.6f\n' % (elname, a[1][0] / angtobohr, + a[1][1] / angtobohr, + a[1][2] / angtobohr)) + + f.close() + + def __build_isosurface(self, filename, cgfs, coeff, isovalue, sz=5, npts=100): + """ + Construct isosurfaces from PyQInt output + """ + # generate some data + isovalue = np.abs(isovalue) + integrator = PyQInt() + grid = integrator.build_rectgrid3d(-sz, sz, npts) + scalarfield = np.reshape(integrator.plot_wavefunction(grid, coeff, cgfs), (npts, npts, npts)) + unitcell = np.diag(np.ones(3) * 2 * sz) + + pytessel = PyTessel() + vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), isovalue) + posfile = filename + '_pos.ply' + pytessel.write_ply(posfile, vertices, normals, indices) + + vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), -isovalue) + negfile = filename + '_neg.ply' + pytessel.write_ply(negfile, vertices, normals, indices) + + return posfile, negfile diff --git a/pyqint/geometry_optimization.py b/pyqint/geometry_optimization.py index e3dd0ad..52a231d 100644 --- a/pyqint/geometry_optimization.py +++ b/pyqint/geometry_optimization.py @@ -127,6 +127,7 @@ def __pack(self, mol, coords): """ newmol = Molecule(mol.name) + newmol.charge = mol.charge coords = coords.reshape((len(mol.atoms), 3)) for i,atom in enumerate(mol.atoms): newmol.add_atom(mol.atoms[i][0], coords[i][0], coords[i][1], coords[i][2]) diff --git a/pyqint/hf.py b/pyqint/hf.py index 474f71e..d19f4c0 100644 --- a/pyqint/hf.py +++ b/pyqint/hf.py @@ -1,7 +1,5 @@ # -*- coding: utf-8 -*- -import json -import os import numpy as np from . import PyQInt import time @@ -31,7 +29,7 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100, # build cgfs, nuclei and calculate nr of electrons cgfs, nuclei = mol.build_basis(basis) - nelec = int(np.sum([at[1] for at in nuclei])) + nelec = mol.get_nelec() N = len(cgfs) occ = [2 if i < nelec//2 else 0 for i in range(N)] @@ -129,7 +127,7 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100, # calculate DIIS coefficients e = (F.dot(P.dot(S)) - S.dot(P.dot(F))).flatten() # calculate error vector - enorm = np.linalg.norm(e) # store error vector norm + #enorm = np.linalg.norm(e) # store error vector norm fmats_diis.append(F) # add Fock matrix to list pmat_diis.append(P) # add density matrix to list evs_diis.append(e) @@ -251,7 +249,7 @@ def rhf_forces(self, mol, basis, C, P, e): # intialization integrator = PyQInt() cgfs, nuclei = mol.build_basis(basis) - nelec = np.sum([nucleus[1] for nucleus in nuclei]) + nelec = mol.get_nelec() forces = np.zeros((len(nuclei),3)) N = len(cgfs) occ = [2 if i < nelec//2 else 0 for i in range(N)] diff --git a/pyqint/molecule.py b/pyqint/molecule.py index d3caccb..9046ba0 100644 --- a/pyqint/molecule.py +++ b/pyqint/molecule.py @@ -13,6 +13,8 @@ def __init__(self, _name='unknown'): self.atoms = [] self.charges = [] self.name = _name + self.charge = 0 + self.nelec = None def __str__(self): res = "Molecule: %s\n" % self.name @@ -21,6 +23,20 @@ def __str__(self): return res + def get_nelec(self): + """ + Get the number of electrons + """ + if self.nelec == None: + raise Exception('You need to use build_basis() before using this function.') + return self.nelec - self.charge + + def set_charge(self, charge): + """ + Set the charge of the molecule + """ + self.charge = charge + def add_atom(self, atom, x, y, z, unit='bohr'): ang2bohr = 1.8897259886 @@ -100,4 +116,6 @@ def build_basis(self, name): for gto in cgf_t['gtos']: self.cgfs[-1].add_gto(gto['coeff'], gto['alpha'], 0, 1, 1) + self.nelec = np.sum(self.charges) + return self.cgfs, self.nuclei diff --git a/tests/test_hf_molecules_charged.py b/tests/test_hf_molecules_charged.py new file mode 100644 index 0000000..bb39324 --- /dev/null +++ b/tests/test_hf_molecules_charged.py @@ -0,0 +1,30 @@ +import unittest +from pyqint import PyQInt, cgf, gto, Molecule, HF, GeometryOptimization +from copy import deepcopy +import numpy as np +import multiprocessing +import os +from nose.tools import nottest + +class TestHFMoleculeCharged(unittest.TestCase): + + def testHCO(self): + """ + Test Hartree-Fock calculation for CO + + Data is compared to results obtained from Gaussian + """ + mol = Molecule() + mol.add_atom('O', -0.12770367, -0.00000000, 1.23419284) + mol.add_atom('C', -0.50625665, 0.00000000, -1.09149431) + mol.add_atom('H', 1.57882331, -0.00000000, -2.06681794) + mol.set_charge(-1) + + results = GeometryOptimization().run(mol, 'sto3g') + + # check that energy matches + self.assertEqual(results['data']['nelec'], 16) + np.testing.assert_almost_equal(results['energies'][-1], -111.523147758727, 5) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_molecule_builder.py b/tests/test_molecule_builder.py index 4495e72..a9dbbdb 100644 --- a/tests/test_molecule_builder.py +++ b/tests/test_molecule_builder.py @@ -23,6 +23,9 @@ def test_load_molecule_from_name(self): np.array([0.6327670,0.6327670,0.6327670]) * 1.8897259886) np.testing.assert_equal(mol.atoms[0][0], 'C') + mol.build_basis('sto3g') + np.testing.assert_equal(mol.get_nelec(), 10) + def test_load_molecule_from_file(self): """ Build a molecule from file @@ -35,6 +38,9 @@ def test_load_molecule_from_file(self): np.testing.assert_almost_equal(mol.atoms[1][1], np.array([0.6327670,0.6327670,0.6327670]) * 1.8897259886) np.testing.assert_equal(mol.atoms[0][0], 'C') + + mol.build_basis('sto3g') + np.testing.assert_equal(mol.get_nelec(), 10) if __name__ == '__main__': unittest.main()