diff --git a/examples/.gitignore b/examples/.gitignore index ef28ff6..7f6150f 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -1,3 +1,5 @@ *.npy *.jpg *.ply +*.png +*.txt diff --git a/examples/blender_render.py b/examples/blender_render.py new file mode 100644 index 0000000..2e2b7ac --- /dev/null +++ b/examples/blender_render.py @@ -0,0 +1,133 @@ +from pyqint import Molecule, PyQInt, FosterBoys, GeometryOptimization, HF +from pyqint import MoleculeBuilder, BlenderRender +import pyqint +import numpy as np +import matplotlib.pyplot as plt +import os +import subprocess + +# +# Plot the isosurfaces for a number of molecules, prior and after +# orbital localization. Note that this script has been designed to be executed +# in a Linux Ubuntu 22.04 LTS (WSL) container. +# + +outpath = os.path.dirname(__file__) + +def main(): + build_orbitals_co() + build_orbitals_ch4() + build_orbitals_ethylene() + +def build_orbitals_co(): + """ + Build a montage image of the canonical and localized molecular orbitals + of CO + """ + molname = 'CO' + mol = Molecule('CO') + mol.add_atom('C', 0, 0, -1.08232106) + mol.add_atom('O', 0, 0, 1.08232106) + res = HF().rhf(mol, 'sto3g') + resfb = FosterBoys(res).run() + + build(molname, res, resfb, nrows=2) + +def build_orbitals_ch4(): + """ + Build a montage image of the canonical and localized molecular orbitals + of CH4 + """ + molname = 'CH4' + mol = MoleculeBuilder().from_name('ch4') + res = GeometryOptimization().run(mol, 'sto3g')['data'] + resfb = FosterBoys(res).run() + + build(molname, res, resfb, nrows=3) + +def build_orbitals_ethylene(): + """ + Build a montage image of the canonical and localized molecular orbitals + of CH4 + """ + molname = 'ethylene' + mol = MoleculeBuilder().from_name('ethylene') + res = GeometryOptimization().run(mol, 'sto3g')['data'] + resfb = FosterBoys(res).run() + + build(molname, res, resfb, nrows=2) + +def build(molname, res, resfb, nrows=2): + """ + Build isosurfaces, montage and print energies to a file + + :param molname: Name of the molecule + :type molname: string + :param res: Results of a Hartree-Fock calculation + :type res: dictionary + :param resfb: Results of a Foster-Boys localization + :type resfb: dictionary + :param nrows: Number of rows in the montage + :type nrows: int + """ + build_isosurfaces(molname, res, resfb) + montage(molname, nrows) + store_energies(os.path.join(os.path.dirname(__file__), 'MO_%s_energies.txt' % molname), res['orbe'], resfb['orbe']) + +def build_isosurfaces(molname, res, resfb): + """ + Builds isosurfaces. + + :param molname: Name of the molecule + :type molname: string + :param res: Result of a Hartree-Fock calculation + :type res: dictionary + :param resfb: Result of a Foster-Boys localization + :type resfb: dictionary + """ + br = BlenderRender() + br.render_molecular_orbitals(res['mol'], res['cgfs'], res['orbc'], outpath, + prefix='MO_CAN_%s' % molname) + + br.render_molecular_orbitals(resfb['mol'], res['cgfs'], resfb['orbc'], outpath, + prefix='MO_FB_%s' % molname) + +def montage(molname, nrows=2): + """ + Produce a montage of the images + + :param molname: Name of the molecule + :type molname: string + :param nrows: Number of rows in the montage + :type nrows: int + """ + out = subprocess.check_output( + ['montage', 'MO_CAN_%s_????.png' % molname, '-tile', 'x%i' % nrows, '-geometry', '128x128+2+2', 'MO_%s_CAN.png' % molname], + cwd=os.path.dirname(__file__) + ) + + out = subprocess.check_output( + ['montage', 'MO_FB_%s_????.png' % molname, '-tile', 'x%i' % nrows, '-geometry', '128x128+2+2', 'MO_%s_FB.png' % molname], + cwd=os.path.dirname(__file__) + ) + +def store_energies(filename, orbe, orbe_fb): + """ + Stores energies. + + :param filename: Name of the molecule + :type filename: string + :param orbe: Array of the canonical molecular orbital energies + :type orbe: list or numpy array + :param orbe_fb: Array of the localized molecular orbital energies + :type orbe_fb: list or numpy array + """ + f = open(filename, 'w') + f.write('id localized canonical\n') + f.write('----------------------------\n') + for i in range(len(orbe)): + f.write("%02i %12.4f %12.4f\n" % (i+1, orbe[i], orbe_fb[i])) + f.close() + +if __name__ == '__main__': + main() diff --git a/examples/foster_boys.py b/examples/foster_boys.py index 513d3a9..5dc58ca 100644 --- a/examples/foster_boys.py +++ b/examples/foster_boys.py @@ -40,6 +40,8 @@ def optimize_co(): res = GeometryOptimization().run(mol, 'sto3g') + print(res['data']['nuclei']) + return res['data'] def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'): diff --git a/meta.yaml b/meta.yaml index 9e9d9d6..4710bb6 100644 --- a/meta.yaml +++ b/meta.yaml @@ -1,6 +1,6 @@ package: name: "pyqint" - version: "0.13.1.0" + version: "0.14.0.0" source: path: . diff --git a/pyqint/_version.py b/pyqint/_version.py index 337d53e..be00ee7 100644 --- a/pyqint/_version.py +++ b/pyqint/_version.py @@ -1,2 +1,2 @@ -__version__ = "0.13.1.0" +__version__ = "0.14.0.0" diff --git a/pyqint/blender/blender_render_molecule.py b/pyqint/blender/blender_render_molecule.py index c986444..7a7b54a 100644 --- a/pyqint/blender/blender_render_molecule.py +++ b/pyqint/blender/blender_render_molecule.py @@ -121,10 +121,13 @@ def set_environment(settings): Specify canvas size, remove default objects, reset positions of camera and light, define film and set background """ + print('Set render engine to: CYCLES') bpy.context.scene.render.engine = 'CYCLES' + print('Set rendering device to GPU') bpy.context.scene.cycles.device = 'GPU' bpy.context.scene.render.resolution_x = settings['resolution'] bpy.context.scene.render.resolution_y = settings['resolution'] + print('Setting resolution to: ', settings['resolution']) bpy.context.scene.cycles.samples = 1024 bpy.context.scene.cycles.tile_size = 2048 diff --git a/pyqint/blender/blender_render_molecules.py b/pyqint/blender/blender_render_molecules.py deleted file mode 100644 index c986444..0000000 --- a/pyqint/blender/blender_render_molecules.py +++ /dev/null @@ -1,245 +0,0 @@ -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 index 1a2fecb..4cb584b 100644 --- a/pyqint/blenderrender.py +++ b/pyqint/blenderrender.py @@ -5,18 +5,23 @@ import json import subprocess import shutil -import tqdm +from sys import platform # 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') + print('Cannot find module PyTessel') try: from mendeleev import element except ModuleNotFoundError: - print('Cannot find mendeleev') + print('Cannot find module mendeleev') + +try: + from tqdm import tqdm +except ModuleNotFoundError: + print('Cannot find module tqdm') class BlenderRender: """ @@ -24,6 +29,10 @@ class BlenderRender: """ def __init__(self): self.executable = self.__find_blender() + self.log = [] + print('******************************************************') + print('WARNING: Blender rendering is an EXPERIMENTAL FEATURE.') + print('******************************************************') if self.executable is None: raise Exception('Cannot find Blender executable') else: @@ -31,7 +40,8 @@ def __init__(self): def render_molecular_orbitals(self, molecule, cgfs, orbc, outpath, mo_indices=None, sz=5, isovalue=0.03, - npts=100): + prefix='MO', npts=100, + negcol='E72F65', poscol='3F9EE7'): if mo_indices is None: # render all orbitals mo_indices = np.arange(0, len(orbc)) @@ -41,17 +51,31 @@ def render_molecular_orbitals(self, molecule, cgfs, orbc, outpath, xyzfile = os.path.join(tempdir, 'mol.xyz') self.__store_xyz(molecule, xyzfile) - for idx in mo_indices: + pbar = tqdm(mo_indices) + for idx in pbar: # build isosurfaces - plyfile = os.path.join(tempdir, 'MO_%04i' % (idx+1)) + pbar.set_description('Producing isosurfaces (#%i)' % (idx+1)) + plyfile = os.path.join(tempdir, '%s_%04i' % (prefix,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) + pbar.set_description('Producing molecular orbital (#%i)' % (idx+1)) + outfile = os.path.join(outpath, '%s_%04i.png' % (prefix,idx+1)) + logoutput = self.__run_blender(plypos, plyneg, xyzfile, outfile, tempdir, negcol, poscol) + + self.log.append("### START LOG: MOLECULAR ORBITAL %i ###" % (idx+1)) + for line in logoutput.splitlines(): + self.log.append(line.decode('utf-8')) + self.log.append("### END LOG: MOLECULAR ORBITAL %i ###" % (idx+1)) shutil.rmtree(tempdir) + # store log in same path + with open(os.path.join(outpath, 'renderlog.txt'), 'w') as f: + for line in self.log: + f.write(line + '\n') + f.close() + def get_executable(self): """ Get the Blender executable @@ -62,20 +86,33 @@ 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 + if platform == "linux" or platform == "linux2": + ex = '/opt/blender-3.3.11-linux-x64/blender' # preferred version and path + if os.path.exists(ex): + return ex + + print('Cannot find proper Blender executable. For Linux, please install Blender LTS 3.3.11 in /opt/blender-3.3.11-linux-x64/.') + print('Blender can be obtained via: https://ftp.nluug.nl/pub/graphics/blender/release/Blender3.3/blender-3.3.11-linux-x64.tar.xz') + print('For more details on how to install Blender, please consult the instructions in the manual: https://pyqint.imc-tue.nl/') + + return None + elif platform == 'win32': + 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 + else: + raise Exception('Your platform is not supported for Blender') return None - def __run_blender(self, negfile, posfile, xyzfile, pngfile, cwd): + def __run_blender(self, negfile, posfile, xyzfile, pngfile, cwd, negcol, poscol): # 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') @@ -83,8 +120,8 @@ def __run_blender(self, negfile, posfile, xyzfile, pngfile, cwd): manifest = { 'mo_colors' : { - 'neg': 'E72F65', - 'pos': '3F9EE7', + 'neg': negcol, + 'pos': poscol, }, 'xyzfile' : xyzfile, 'mo_name' : 'isosurface', @@ -115,6 +152,8 @@ def __run_blender(self, negfile, posfile, xyzfile, pngfile, cwd): cwd=cwd ) + return out + def __store_xyz(self, mol, filename): """ Convert final result of GeometryOptimization class to an .xyz file diff --git a/pyqint/foster_boys.py b/pyqint/foster_boys.py index 337b00e..21f1b19 100644 --- a/pyqint/foster_boys.py +++ b/pyqint/foster_boys.py @@ -14,6 +14,7 @@ def __init__(self, res, seed=None): # copy objects from Hartree-Fock result dictionary self.orbc_canonical = res['orbc'] self.orbe_canonical = res['orbe'] + self.mol = res['mol'] self.nelec = res['nelec'] self.H = res['fock'] self.cgfs = res['cgfs'] @@ -69,6 +70,7 @@ def __single_runner(self): 'orbe': orbe, 'orbc': orbc, 'nriter': nriter, + 'mol': self.mol, 'r2start': self.__calculate_r2(self.orbc_canonical), 'r2final': self.__calculate_r2(orbc) } diff --git a/pyqint/hf.py b/pyqint/hf.py index d19f4c0..9ac6c8b 100644 --- a/pyqint/hf.py +++ b/pyqint/hf.py @@ -195,6 +195,7 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100, "ex": -0.25 * np.einsum('iklj,ij,kl', tetensor, P, P), "enucrep": nuc_rep, "nelec": nelec, + "mol": mol, "forces": self.rhf_forces(mol, basis, C, P, orbe) if calc_forces else None } diff --git a/setup.py b/setup.py index d00e100..8a072d4 100644 --- a/setup.py +++ b/setup.py @@ -124,7 +124,7 @@ def find_windows_versions(): ext_modules=cythonize(ext_modules[0], language_level = "3", build_dir="build"), - packages=['pyqint', 'pyqint.basissets', 'pyqint.molecules'], + packages=['pyqint', 'pyqint.basissets', 'pyqint.molecules', 'pyqint.blender'], include_package_data=True, classifiers=[ "Programming Language :: Python :: 3",