Source code for psi4.driver.p4util.inpsight

#
# @BEGIN LICENSE
#
# Psi4: an open-source quantum chemistry software package
#
# Copyright (c) 2007-2023 The Psi4 Developers.
#
# The copyrights for code used from other parties are included in
# the corresponding files.
#
# This file is part of Psi4.
#
# Psi4 is free software; you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, version 3.
#
# Psi4 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with Psi4; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# @END LICENSE
#

__all__ = ["InPsight"]

import math
import os
from datetime import date


# yapf: disable
[docs] class InPsight: """POV-Ray visualization.""" # POV-Ray defines defines = {} defines['Shadows'] = 'false' defines['Background_Color'] = '<0.6,0.6,0.6>' defines['Output_File_Type'] = 'N' defines['Output_Alpha'] = 'true' defines['Light_Color'] = '<1,1,1>' defines['Filename'] = 'inpsight' defines['Filepath'] = os.getcwd() defines['Antialias'] = 'true' defines['Antialias_Threshold'] = '0.1' # Molecule geometry atoms = [] # (Z,x,y,z,R,r,g,b,t) in bohr bonds = [] # (x1,y1,z1,R1,x2,y2,z2,R2,r,g,b,t) # Molecular geometry defines colors = [] radii = [] radial_scale = 0.25 bond_width = 0.2 # bohr bohr_per_ang = 1.8897161646320724 bonding_alpha = 0.65 # Used to select/reject bonds via sum of vDW radii # View defines (high-level) azimuth = 0.0 elevation = 0.0 zoom = 0.5 height = 900 width = 1200 # Camera positions (low-level) location = [1.0,0.0,0.0] up = [0.0,0.75,0.0] right = [1.0,0.0,0.0] sky = [0.0,-1.0,0.0] look_at = [0.0,0.0,0.0] light = [1.0,0.0,0.0] light_color = [0.6,0.6,0.6] # Standard Jmol colors, 256-based colors.append([0,0,0]) colors.append([255,255,255]) colors.append([217,255,255]) colors.append([204,128,255]) colors.append([194,255,0]) colors.append([255,181,181]) colors.append([144,144,144]) colors.append([48,80,248]) colors.append([255,13,13]) colors.append([144,224,80]) colors.append([179,227,245]) colors.append([171,92,242]) colors.append([138,255,0]) colors.append([191,166,166]) colors.append([240,200,160]) colors.append([255,128,0]) colors.append([255,255,48]) colors.append([31,240,31]) colors.append([128,209,227]) colors.append([143,64,212]) colors.append([61,255,0]) colors.append([230,230,230]) colors.append([191,194,199]) colors.append([166,166,171]) colors.append([138,153,199]) colors.append([156,122,199]) colors.append([224,102,51]) colors.append([240,144,160]) colors.append([80,208,80]) colors.append([200,128,51]) colors.append([125,128,176]) colors.append([194,143,143]) colors.append([102,143,143]) colors.append([189,128,227]) colors.append([255,161,0]) colors.append([166,41,41]) colors.append([92,184,209]) colors.append([112,46,176]) colors.append([0,255,0]) colors.append([148,255,255]) colors.append([148,224,224]) colors.append([115,194,201]) colors.append([84,181,181]) colors.append([59,158,158]) colors.append([36,143,143]) colors.append([10,125,140]) colors.append([0,105,133]) colors.append([192,192,192]) colors.append([255,217,143]) colors.append([166,117,115]) colors.append([102,128,128]) colors.append([158,99,181]) colors.append([212,122,0]) colors.append([148,0,148]) colors.append([66,158,176]) colors.append([87,23,143]) colors.append([0,201,0]) colors.append([112,212,255]) colors.append([255,255,199]) colors.append([217,255,199]) colors.append([199,255,199]) colors.append([163,255,199]) colors.append([143,255,199]) colors.append([97,255,199]) colors.append([69,255,199]) colors.append([48,255,199]) colors.append([31,255,199]) colors.append([0,255,156]) colors.append([0,230,117]) colors.append([0,212,82]) colors.append([0,191,56]) colors.append([0,171,36]) colors.append([77,194,255]) colors.append([77,166,255]) colors.append([33,148,214]) colors.append([38,125,171]) colors.append([38,102,150]) colors.append([23,84,135]) colors.append([208,208,224]) colors.append([255,209,35]) colors.append([184,184,208]) colors.append([166,84,77]) colors.append([87,89,97]) colors.append([158,79,181]) colors.append([171,92,0]) colors.append([117,79,69]) colors.append([66,130,150]) colors.append([66,0,102]) colors.append([0,125,0]) colors.append([112,171,250]) colors.append([0,186,255]) colors.append([0,161,255]) colors.append([0,143,255]) colors.append([0,128,255]) colors.append([0,107,255]) colors.append([84,92,242]) colors.append([120,92,227]) colors.append([138,79,227]) colors.append([161,54,212]) colors.append([179,31,212]) colors.append([179,31,186]) colors.append([179,13,166]) colors.append([189,13,135]) colors.append([199,0,102]) colors.append([204,0,89]) colors.append([209,0,79]) colors.append([217,0,69]) colors.append([224,0,56]) colors.append([230,0,46]) colors.append([235,0,38]) # Approximate vDW radii in angstrom radii.append(2.0) radii.append(1.001) radii.append(1.012) radii.append(0.825) radii.append(1.408) radii.append(1.485) radii.append(1.452) radii.append(1.397) radii.append(1.342) radii.append(1.287) radii.append(1.243) radii.append(1.144) radii.append(1.364) radii.append(1.639) radii.append(1.716) radii.append(1.705) radii.append(1.683) radii.append(1.639) radii.append(1.595) radii.append(1.485) radii.append(1.474) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.562) radii.append(1.650) radii.append(1.727) radii.append(1.760) radii.append(1.771) radii.append(1.749) radii.append(1.727) radii.append(1.628) radii.append(1.606) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.639) radii.append(1.672) radii.append(1.804) radii.append(1.881) radii.append(1.892) radii.append(1.892) radii.append(1.881) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) radii.append(2.0) def __init__(self,molecule): self.molecule = molecule self.molecule.update_geometry() self.update_geometry()
[docs] def update_geometry(self): # Atoms natom = self.molecule.natom() self.atoms = [] for k in range(0,natom): x = self.molecule.x(k) y = self.molecule.y(k) z = self.molecule.z(k) Z = self.molecule.Z(k) atom = Z, x, y, z, self.radial_scale * self.bohr_per_ang * self.radii[Z], self.colors[Z][0] / 256.0, \ self.colors[Z][1] / 256.0, self.colors[Z][2] / 256.0, 0.0 self.atoms.append(atom) # Bonds self.bonds = [] for k in range(1,natom): for l in range (0, k): Z1 = self.atoms[k][0] Z2 = self.atoms[l][0] R1 = self.bohr_per_ang*self.radii[Z1] R2 = self.bohr_per_ang*self.radii[Z2] x1 = self.atoms[k][1] y1 = self.atoms[k][2] z1 = self.atoms[k][3] x2 = self.atoms[l][1] y2 = self.atoms[l][2] z2 = self.atoms[l][3] r1 = self.atoms[k][5] g1 = self.atoms[k][6] b1 = self.atoms[k][7] t1 = self.atoms[k][8] r2 = self.atoms[l][5] g2 = self.atoms[l][6] b2 = self.atoms[l][7] t2 = self.atoms[l][8] R = math.sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)) if (R < self.bonding_alpha*(R1 + R2)): omega = R2 / (R1 + R2) xc = omega * (x1 - x2) + x2 yc = omega * (y1 - y2) + y2 zc = omega * (z1 - z2) + z2 bond1 = x1,y1,z1,self.bond_width, xc,yc,zc,self.bond_width,r1,g1,b1,t1 bond2 = x2,y2,z2,self.bond_width, xc,yc,zc,self.bond_width,r2,g2,b2,t2 self.bonds.append(bond1) self.bonds.append(bond2)
[docs] def set_define(self, key, value): self.defines[key] = value
[docs] def set_color(self, Z, color): self.colors[Z] = color
[docs] def set_radius(self, Z, radius): self.radii[Z] = radius
[docs] def position_camera(self): xc = self.molecule.center_of_mass() self.look_at = [xc[0], xc[1], xc[2]] Rmax = 0.0 natom = self.molecule.natom() for k in range(0,natom): x = [self.molecule.x(k), self.molecule.y(k), self.molecule.z(k)] R = math.sqrt((x[0] - xc[0])*(x[0] - xc[0]) + (x[1] - xc[1])*(x[1] - xc[1]) + (x[2] - xc[2])*(x[2] - xc[2])) if R > Rmax: Rmax = R Rmax = Rmax / self.zoom if (self.width < self.height): self.right = [Rmax, 0.0, 0.0] self.up = [0.0, self.right[0]*self.height/self.width, 0.0] else: self.up = [0.0, Rmax, 0.0] self.right = [self.up[1]*self.width/self.height, 0.0, 0.0] phi = math.pi*(-self.azimuth)/180.0 theta = math.pi*(90.0 - self.elevation)/180.0 delta = [Rmax*math.cos(phi)*math.sin(theta), Rmax*math.sin(phi)*math.sin(theta), Rmax*math.cos(theta)] self.location = [xc[0] + delta[0], xc[1] + delta[1], xc[2] + delta[2]] phi = math.pi*(-(self.azimuth + 30.0))/180.0 theta = math.pi*(90.0 - (self.elevation + 30.0))/180.0 delta = [Rmax*math.cos(phi)*math.sin(theta), Rmax*math.sin(phi)*math.sin(theta), Rmax*math.cos(theta)] self.light = [xc[0] + delta[0], xc[1] + delta[1], xc[2] + delta[2]]
[docs] def set_view(self,azimuth, elevation, zoom = 0.7): self.azimuth = azimuth self.elevation = elevation self.zoom = zoom self.position_camera()
[docs] def set_size(self, width,height): self.width = width self.height = height
[docs] def set_camera(self, location, sky, up, right, look_at, light, light_color): self.location = location self.sky = sky self.up = up self.right = right self.look_at = look_at self.light = light self.light_color = light_color
[docs] def save_molecule(self, filename): if (filename != ''): self.defines['Filename'] = filename ini_filename = self.defines['Filepath'] + '/' + self.defines['Filename'] + '.pov.ini' pov_filename = self.defines['Filepath'] + '/' + self.defines['Filename'] + '.pov' png_filename = self.defines['Filepath'] + '/' + self.defines['Filename'] + '.png' pov_file = self.defines['Filename'] + '.pov' png_file = self.defines['Filename'] + '.png' # Write the pov.ini file fh = open(ini_filename,'w') fh.write('; InPsight: visualization in Psi4\n') fh.write('; by Rob Parrish\n') fh.write('; .pov.ini file\n') fh.write('; Created %s\n' % str(date.today())) fh.write('\n') fh.write('Input_File_Name=%s\n' % pov_file) fh.write('Output_to_File=true\n') fh.write('Output_File_Type=%s\n' % self.defines['Output_File_Type']) fh.write('Output_File_Name=%s\n' % png_file) fh.write('Height=%s\n' % str(self.height)) fh.write('Width=%s\n' % str(self.width)) fh.write('Output_Alpha=%s\n' % self.defines['Output_Alpha']) fh.write('Antialias=%s\n' % self.defines['Antialias']) fh.write('Antialias_Threshold=%s\n' % self.defines['Antialias_Threshold']) fh.write('Display=true\n') fh.write('Warning_Level=5\n') fh.write('Verbose=false\n') fh.close() # Write the pov file fh = open(pov_filename, 'w') fh.write('// InPsight: visualization in Psi4\n') fh.write('// by Rob Parrish\n') fh.write('// .pov file (adopted from Jmol)\n') fh.write('// Created %s\n' % str(date.today())) fh.write('#declare Width = %s;\n' % str(self.width)) fh.write('#declare Height = %s;\n' % str(self.height)) fh.write('#declare Shadows = %s; \n' % self.defines['Shadows']) fh.write('\n') fh.write('camera{\n') fh.write(' orthographic\n') fh.write(' location < %s, %s, %s>\n' % (str(self.location[0]),str(self.location[1]),str(self.location[2]) )) fh.write(' sky < %s, %s, %s>\n' % (str(self.sky[0]), str(self.sky[1]), str(self.sky[2]) )) fh.write(' up < %s, %s, %s>\n' % (str(self.up[0]), str(self.up[1]), str(self.up[2]) )) fh.write(' right < %s, %s, %s>\n' % (str(self.right[0]), str(self.right[1]), str(self.right[2]) )) fh.write(' look_at < %s, %s, %s>\n' % (str(self.look_at[0]), str(self.look_at[1]), str(self.look_at[2]) )) fh.write('}\n') fh.write('\n') fh.write('background { color rgb %s }\n' % self.defines['Background_Color']) fh.write('light_source { <%s,%s,%s> rgb <%s,%s,%s> }\n' % (str(self.light[0]),str(self.light[1]),str(self.light[2]), str(self.light_color[0]),str(self.light_color[1]),str(self.light_color[2]))) fh.write('\n') fh.write('// ***********************************************\n') fh.write('// macros for atom/bond shapes\n') fh.write('// ***********************************************\n') fh.write('#macro check_shadow()\n') fh.write(' #if (!Shadows)\n') fh.write(' no_shadow \n') fh.write(' #end\n') fh.write('#end\n') fh.write('\n') fh.write('#macro translucentFinish(T)\n') fh.write(' #local shineFactor = T;\n') fh.write(' #if (T <= 0.25)\n') fh.write(' #declare shineFactor = (1.0-4*T);\n') fh.write(' #end\n') fh.write(' #if (T > 0.25)\n') fh.write(' #declare shineFactor = 0;\n') fh.write(' #end\n') fh.write(' finish {\n') fh.write(' ambient 0.45\n') fh.write(' diffuse 0.84\n') fh.write(' specular 0.22\n') fh.write(' roughness .00001\n') fh.write(' metallic shineFactor\n') fh.write(' phong 0.9*shineFactor\n') fh.write(' phong_size 120*shineFactor\n') fh.write('}#end\n') fh.write('\n') fh.write('#macro a(X,Y,Z,RADIUS,R,G,B,T)\n') fh.write(' sphere{<X,Y,Z>,RADIUS\n') fh.write(' pigment{rgbt<R,G,B,T>}\n') fh.write(' translucentFinish(T)\n') fh.write(' check_shadow()}\n') fh.write('#end\n') fh.write('\n') fh.write('#macro b(X1,Y1,Z1,RADIUS1,X2,Y2,Z2,RADIUS2,R,G,B,T)\n') fh.write(' cone{<X1,Y1,Z1>,RADIUS1,<X2,Y2,Z2>,RADIUS2\n') fh.write(' pigment{rgbt<R,G,B,T>}\n') fh.write(' translucentFinish(T)\n') fh.write(' check_shadow()}\n') fh.write('#end \n') for bond in self.bonds: fh.write('b(%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s)\n' % (str(bond[0]),str(bond[1]),str(bond[2]),str(bond[3]), str(bond[4]),str(bond[5]),str(bond[6]),str(bond[7]), str(bond[8]),str(bond[9]),str(bond[10]),str(bond[11]))) for atom in self.atoms: fh.write('a(%s,%s,%s,%s,%s,%s,%s,%s)\n' % (str(atom[1]),str(atom[2]),str(atom[3]),str(atom[4]), str(atom[5]),str(atom[6]),str(atom[7]),str(atom[8]))) fh.close()
[docs] def save_density(self,filename='rho',overlap = 2.0,n = [40,40,40],caxis = [0.0,1.0]): if (filename != ''): self.defines['Filename'] = filename # grid = GridProp() # GridProp seems to have been retired grid = None grid.set_n(n[0],n[1],n[2]) grid.set_caxis(caxis[0],caxis[1]) grid.set_filename(self.defines['Filename']) grid.add('RHO') grid.compute() df3_file = filename + '.RHO.df3' l = [grid.get_l(0),grid.get_l(1),grid.get_l(2)] o = [grid.get_o(0),grid.get_o(1),grid.get_o(2)] ini_filename = self.defines['Filepath'] + '/' + self.defines['Filename'] + '.pov.ini' pov_filename = self.defines['Filepath'] + '/' + self.defines['Filename'] + '.pov' png_filename = self.defines['Filepath'] + '/' + self.defines['Filename'] + '.png' pov_file = self.defines['Filename'] + '.pov' png_file = self.defines['Filename'] + '.png' # Write the pov.ini file fh = open(ini_filename,'w') fh.write('; InPsight: visualization in Psi4\n') fh.write('; by Rob Parrish\n') fh.write('; .pov.ini file\n') fh.write('; Created %s\n' % str(date.today())) fh.write('\n') fh.write('Input_File_Name=%s\n' % pov_file) fh.write('Output_to_File=true\n') fh.write('Output_File_Type=%s\n' % self.defines['Output_File_Type']) fh.write('Output_File_Name=%s\n' % png_file) fh.write('Height=%s\n' % str(self.height)) fh.write('Width=%s\n' % str(self.width)) fh.write('Output_Alpha=%s\n' % self.defines['Output_Alpha']) fh.write('Antialias=%s\n' % self.defines['Antialias']) fh.write('Antialias_Threshold=%s\n' % self.defines['Antialias_Threshold']) fh.write('Display=true\n') fh.write('Warning_Level=5\n') fh.write('Verbose=false\n') fh.close() # Write the pov file fh = open(pov_filename, 'w') fh.write('// InPsight: visualization in Psi4\n') fh.write('// by Rob Parrish\n') fh.write('// .pov file (adopted from Jmol)\n') fh.write('// Created %s\n' % str(date.today())) fh.write('#declare Shadows = %s; \n' % self.defines['Shadows']) fh.write('\n') fh.write('camera{\n') fh.write(' orthographic\n') fh.write(' location < %s, %s, %s>\n' % (str(self.location[0]),str(self.location[1]),str(self.location[2]) )) fh.write(' sky < %s, %s, %s>\n' % (str(self.sky[0]), str(self.sky[1]), str(self.sky[2]) )) fh.write(' up < %s, %s, %s>\n' % (str(self.up[0]), str(self.up[1]), str(self.up[2]) )) fh.write(' right < %s, %s, %s>\n' % (str(self.right[0]), str(self.right[1]), str(self.right[2]) )) fh.write(' look_at < %s, %s, %s>\n' % (str(self.look_at[0]), str(self.look_at[1]), str(self.look_at[2]) )) fh.write('}\n') fh.write('\n') fh.write('background { color rgb %s }\n' % self.defines['Background_Color']) fh.write('light_source { <%s,%s,%s> rgb <%s,%s,%s> }\n' % (str(self.light[0]),str(self.light[1]),str(self.light[2]), str(self.light_color[0]),str(self.light_color[1]),str(self.light_color[2]))) fh.write('\n') fh.write('// ***********************************************\n') fh.write('// macros for atom/bond shapes\n') fh.write('// ***********************************************\n') fh.write('#macro check_shadow()\n') fh.write(' #if (!Shadows)\n') fh.write(' no_shadow \n') fh.write(' #end\n') fh.write('#end\n') fh.write('\n') fh.write('#macro translucentFinish(T)\n') fh.write(' #local shineFactor = T;\n') fh.write(' #if (T <= 0.25)\n') fh.write(' #declare shineFactor = (1.0-4*T);\n') fh.write(' #end\n') fh.write(' #if (T > 0.25)\n') fh.write(' #declare shineFactor = 0;\n') fh.write(' #end\n') fh.write(' finish {\n') fh.write(' ambient 0.45\n') fh.write(' diffuse 0.84\n') fh.write(' specular 0.22\n') fh.write(' roughness .00001\n') fh.write(' metallic shineFactor\n') fh.write(' phong 0.9*shineFactor\n') fh.write(' phong_size 120*shineFactor\n') fh.write('}#end\n') fh.write('\n') fh.write('#macro a(X,Y,Z,RADIUS,R,G,B,T)\n') fh.write(' sphere{<X,Y,Z>,RADIUS\n') fh.write(' pigment{rgbt<R,G,B,T>}\n') fh.write(' translucentFinish(T)\n') fh.write(' check_shadow()}\n') fh.write('#end\n') fh.write('\n') fh.write('#macro b(X1,Y1,Z1,RADIUS1,X2,Y2,Z2,RADIUS2,R,G,B,T)\n') fh.write(' cone{<X1,Y1,Z1>,RADIUS1,<X2,Y2,Z2>,RADIUS2\n') fh.write(' pigment{rgbt<R,G,B,T>}\n') fh.write(' translucentFinish(T)\n') fh.write(' check_shadow()}\n') fh.write('#end \n') for bond in self.bonds: fh.write('b(%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s)\n' % (str(bond[0]),str(bond[1]),str(bond[2]),str(bond[3]), str(bond[4]),str(bond[5]),str(bond[6]),str(bond[7]), str(bond[8]),str(bond[9]),str(bond[10]),str(bond[11]))) for atom in self.atoms: fh.write('a(%s,%s,%s,%s,%s,%s,%s,%s)\n' % (str(atom[1]),str(atom[2]),str(atom[3]),str(atom[4]), str(atom[5]),str(atom[6]),str(atom[7]),str(atom[8]))) fh.close()
# yapf: enable