import os import sys import math import threading from tkinter import * try: import tklib.tkimport as imp except Exception as e: print() print("######################################################################") print("########### ERROR ERROR ERROR ERROR ERROR ERROR #####################") print("######################################################################") print(f"# Failed to import [tklib.tkimport] module ({e}).") print(f"# Add [tkProg]{os.sep}tklib{os.sep}python to PYTHONPATH variable") print(f"# Current PYTHONPATH:", sys.path) print("######################################################################") input("Press ENTER to terminate>>") exit() np = imp.import_lib("numpy", stop_by_error = False) pmg = imp.import_lib("pymatgen", stop_by_error = False) ogl = imp.import_lib("OpenGL", stop_by_error = False) pil = imp.import_lib("PIL", stop_by_error = False) imp.messages(stop_by_error = True) from numpy import sqrt, exp, log, log10 from OpenGL.GL import * from OpenGL.GLUT import * from OpenGL.GLU import * from PIL import Image #from pymatgen.core.structure import Structure #from pymatgen.core import Composition #from pymatgen.symmetry.analyzer import SpacegroupAnalyzer #from pymatgen.io.cif import CifWriter, CifParser from tklib.tkprogvars import VESTADBDir from tklib.tkutils import pint, pfloat #from tklib.tkutils import print_data, pint, pfloat #from tklib.tkvariousdata import tkVariousData from tklib.tkapplication import tkApplication from tklib.tkparams import tkParams from tklib.tkcrystal.tkcif import tkCIF from tklib.tkcrystal.tkcrystal import tkCrystal from tklib.tkgraphic.tkopengl import draw_sphere, draw_cylinder, draw_arrow def save_image(filename, width, height): glPixelStorei(GL_PACK_ALIGNMENT, 1) data = glReadPixels(0, 0, width, height, GL_RGB, GL_UNSIGNED_BYTE) image = Image.frombytes("RGB", (width, height), data) image = image.transpose(Image.FLIP_TOP_BOTTOM) image.save(filename) def read_cif(infile, app = None, print_level = 1): if print_level: print("") print(f"Read [{infile}]") if not os.path.isfile(infile): app.terminate(f"Error in read_cif(): Can not read [{infile}]") cif = tkCIF() cifdata = cif.ReadCIF(infile, find_valid_structure = 1) cif.Close() if not cifdata: if print_level: terminate(f"Error in read_cif(): Could not get cifdata from ciffile [{ciffile}]") return None cry = cifdata.GetCrystal() # cry.PrintInf() return cry def read_VESTA_elements(path = None, print_level = 1): if path is None: path = os.path.join(VESTADBDir, 'elements.ini') if print_level: print(f"\nread_VESTA_elements(): Read atom data from [{path}]") if path is None or not os.path.exists(path): return None atom_dict = {} fp = open(path, "r") for line in fp: a = line.split() atom_name = a[1] _inf = { "Z": pfloat(a[0]), "atom_name": atom_name, "ratom": pfloat(a[2]), "rvdw" : pfloat(a[3]), "rion" : pfloat(a[4]), "color": [pfloat(a[5]), pfloat(a[6]), pfloat(a[7])], } atom_dict[atom_name] = _inf fp.close() return atom_dict def draw_text(text, x, y, z, sx, sy, sz, R = 0.5, G = 0.5, B = 0.5): glColor3f(R, G, B) glRasterPos3f(x, y, z) glPushMatrix() # glTranslatef(0, 0, 0) glScalef(sx, sy, sz) for ch in text: glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18, ord(ch)) glPopMatrix() def initialize(app): app.cfg = tkParams() app.cfg.figsize = [8, 6] arg_config_file = app.replace_path(None, template = ["{dirname}", "{filebody}_arg_config.xlsx"]) print(f"Read config file [{arg_config_file}]") ret = app.read_arg_config_from_file(arg_config_file) if not arg_config_file: app.terminate("\nError in initialize(): Can not read config file [{arg_config_file}]") return app.cfg class tkCrystalView(): def __init__(self, app = None, cry = None, center = [0, 0, 0], r_mode = 'rion', max_distance = 2.7, draw_range = [2, 2, 2], bond_r = 0.1, k_rion = 0.3, color_bg = [1.0, 0.0, 1.0, 1.0], alpha_cell = 0.5, alpha_atom = 0.8, alpha_bond = 0.5, alpha_vector = 0.5, nslices = 32, nstacks = 32, VESTA_elements_db = None, title = b"Crystal Structure", window_w = 800, window_h = 600, debug = False): self.app = app self.width = window_w self.height = window_h self.title = title self.debug = debug self.print_level = 1 if debug else 0 # 状態変数 self.model_x = 0.0 self.model_y = 0.0 self.model_z = -10.0 self.mouse_down = False self.mouse_x = 0.0 self.mouse_y = 0.0 self.angle_x = 0.0 self.angle_y = 0.0 self.angle_z = 0.0 self.scale = 1.0 # 描画設定 self.color_bg = color_bg # 結晶構造 self.center = center self.cry = cry self.r_mode = r_mode self.k_rion = k_rion self.max_distance = max_distance self.bond_r = bond_r self.draw_range = draw_range # self.nxrange = (-1, 2) # self.nyrange = (-1, 2) # self.nzrange = (-1, 2) self.nslices = nslices self.nstacks = nstacks self.VESTA_elements_db = VESTA_elements_db self.atom_dict = read_VESTA_elements(VESTA_elements_db, print_level = self.print_level) self.alpha_cell = alpha_cell self.alpha_atom = alpha_atom self.alpha_bond = alpha_bond self.alpha_vector = alpha_vector self.init_canvas(color_bg = self.color_bg, title = self.title) def init_canvas(self, color_bg, title): glutInit(sys.argv) glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH) glutInitWindowSize(self.width, self.height) self.hWindow = glutCreateWindow(title) print("window handle: ", self.hWindow) glClearColor(*color_bg) glEnable(GL_BLEND) glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA) glEnable(GL_DEPTH_TEST) # glShadeModel(GL_FLAT) glShadeModel(GL_SMOOTH) # ライトの設定 glEnable(GL_LIGHTING) # ライティングの有効化 glEnable(GL_LIGHT0) # ライト0の有効化 light_diffuse = [1.0, 1.0, 1.0, 1.0] # 拡散光の色 light_specular = [1.0, 1.0, 1.0, 1.0] # 鏡面反射光の色 glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse) glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular) # マテリアルの設定 mat_ambient = [0.2, 0.2, 0.2, 1.0] # 環境光反射成分 mat_diffuse = [1.0, 1.0, 1.0, 1.0] # 拡散反射成分 mat_specular = [1.0, 1.0, 1.0, 1.0] # 鏡面反射成分 mat_shininess = [50.0] # 輝き glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient) glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse) glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular) glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess) self.init_projection(self.width, self.height) # glEnable(GL_DEPTH_TEST) def camera_distance(self, W, D, H, fov_y): diagonal = math.sqrt(W**2 + H**2 + D**2) fov_rad = math.radians(fov_y) distance = diagonal / (2.0 * math.tan(fov_rad / 2.0)) return distance def init_projection(self, width, height): fov_y = 45.0 glMatrixMode(GL_PROJECTION) glLoadIdentity() glViewport(0, 0, width, height) gluPerspective(fov_y, width / height, 0.1, 100.0) # gluOrtho2D(0, width, 0, height) camera_distance = self.camera_distance(*(self.center * 2 + 3.0), fov_y) glMatrixMode(GL_MODELVIEW) glLoadIdentity() glTranslatef(-self.center[0], -self.center[1], -self.center[2]) gluLookAt(*np.array([-0.5, 2.0, 0.2]) * camera_distance, # カメラの位置 0.0, 0.0, 0.0, # 注視点の位置 0.0, 0.0, 1.0) # 上方向のベクトル glRotatef(180, 0, 0, 1) # glRotatef(90, 1, 0, 0) light_position = [3.0, -camera_distance, 5.0, 0.0] # ライトの位置 glLightfv(GL_LIGHT0, GL_POSITION, light_position) # callback for mouse/keyboard def reshape(self, w, h): self.width = w self.height = h self.init_projection(w, h) def mouse(self, button, state, x, y): if button == GLUT_LEFT_BUTTON: if state == GLUT_DOWN: self.mouse_down = True self.mouse_x = x self.mouse_y = y elif state == GLUT_UP: self.mouse_down = False def motion(self, x, y): if self.mouse_down: self.angle_x -= (y - self.mouse_y) self.angle_z += (x - self.mouse_x) self.mouse_x = x self.mouse_y = y glutPostRedisplay() if self.debug: print(f"mouse: ({self.mouse_x}, {self.mouse_y})") print(f"angle: ({self.angle_x}, {self.angle_y}, {self.angle_z})") def keyboard(self, key, x, y): if key == b'u' or key == b'U': self.scale += 0.1 elif key == b'd' or key == b'D': self.scale -= 0.1 if self.scale < 0.1: self.scale = 0.1 elif key == b'h' or key == b'H': self.model_x -= 0.1 elif key == b'l' or key == b'L': self.model_x += 0.1 elif key == b'j' or key == b'J': self.model_z -= 0.1 elif key == b'k' or key == b'K': self.model_z += 0.1 elif key == b'<' or key == b',': self.model_y -= 2.0 elif key == b'>' or key == b'.': self.model_y += 2.0 glutPostRedisplay() if self.debug: print(f"scale: {self.scale}") def draw_cell(self, aij, origin, radius = 0.1, R = 0.7, G = 0.7, B = 0.7, alpha = 0.5): nx, ny, nz = self.draw_range nx = pint(nx + 0.9999) ny = pint(ny + 0.9999) nz = pint(nz + 0.9999) o = origin a000 = o a100 = o + aij[0] a101 = o + aij[0] + aij[2] a110 = o + aij[0] + aij[1] a111 = o + aij[0] + aij[1] + aij[2] a010 = o + aij[1] a011 = o + aij[1] + aij[2] a001 = o + aij[2] kwargs = { "radius": radius, "R": R, "G": G, "B": B, "alpha": alpha, "slices": self.nslices, "stacks": self.nstacks } for _nx in range(0, nx): for _ny in range(0, ny): for _nz in range(0, nz): _a000 = a000 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] _a100 = a100 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] _a101 = a101 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] _a110 = a110 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] _a111 = a111 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] _a010 = a010 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] _a011 = a011 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] _a001 = a001 + _nx * aij[0] + _ny * aij[1] + _nz * aij[2] draw_cylinder(*_a000, *_a100, **kwargs) draw_cylinder(*_a000, *_a010, **kwargs) draw_cylinder(*_a000, *_a001, **kwargs) if True: # if _nx == nx - 1 and _ny == ny - 1: draw_cylinder(*_a110, *_a111, **kwargs) # elif _nx == nx - 1 and _nz == nz - 1: draw_cylinder(*_a101, *_a111, **kwargs) # elif _ny == ny - 1 and _nz == nz - 1: draw_cylinder(*_a011, *_a111, **kwargs) # elif _nx == nx - 1: draw_cylinder(*_a100, *_a110, **kwargs) draw_cylinder(*_a100, *_a101, **kwargs) # elif _ny == ny - 1: draw_cylinder(*_a010, *_a110, **kwargs) draw_cylinder(*_a010, *_a011, **kwargs) # elif _nz == nz - 1: draw_cylinder(*_a001, *_a101, **kwargs) draw_cylinder(*_a001, *_a011, **kwargs) def draw_vectors(self, sites_drawn, length = 0.6, radius = 0.1, arrow_h = 0.3, arrow_r = 0.3, R = 0.0, G = 0.5, B = 0.5, alpha = 0.5): for i, site_drawn in enumerate(sites_drawn): i, site, nx, ny, nz, x, y, z, xc, yc, zc = site_drawn atom_name = site.atom_name_only() # print(f"#{i}: {atom_name} ({x}, {y}, {z}) ({xc}, {yc}, {zc})") atom_inf = self.atom_dict[atom_name] # print("atom_inf=", atom_inf) color = atom_inf["color"] kwargs = { "radius": radius, "arrow_h": arrow_h, "arrow_r": arrow_r, "R": color[0], "G": color[1], "B": color[2], "alpha": alpha, "slices": self.nslices, "stacks": self.nstacks } draw_arrow(xc, yc, zc, xc + length, yc + length, zc, **kwargs) def draw_atoms(self, sites_drawn, alpha = 1.0): cry = self.cry kwargs = { "alpha": alpha, "slices": self.nslices, "stacks": self.nstacks } for i, site_drawn in enumerate(sites_drawn): i, site, nx, ny, nz, x, y, z, xc, yc, zc = site_drawn atom_name = site.atom_name_only() # print(f"#{i}: {atom_name} ({x}, {y}, {z}) ({xc}, {yc}, {zc})") atom_inf = self.atom_dict[atom_name] # print("atom_inf=", atom_inf) rion = atom_inf[self.r_mode] * self.k_rion color = atom_inf["color"] draw_sphere(xc, yc, zc, radius = rion, R = color[0], G = color[1], B = color[2], **kwargs) return sites_drawn def get_drawn_sites(self): eps = 0.05 cry = self.cry max_distance = self.app.cfg.get("max_distance", 10.0) nx, ny, nz = self.draw_range inx = pint(nx + 0.9999) iny = pint(ny + 0.9999) inz = pint(nz + 0.9999) sites = cry.expanded_atom_site_list() sites_drawn = [] for i, site in enumerate(sites): atom_name = site.atom_name_only() x, y, z = site.position(is_reduce01 = 1) for _nx in range(-1, inx + 1): for _ny in range(-1, iny + 1): for _nz in range(-1, inz + 1): _x = x + _nx _y = y + _ny _z = z + _nz # print(f"#{i}: {atom_name} ({x}, {y}, {z}) ({_x}, {_y}, {_z})") if _x < -eps or nx + eps < _x: continue if _y < -eps or ny + eps < _y: continue if _z < -eps or nz + eps < _z: continue # print(f" to be drawn") xc, yc, zc = cry.fractional2cartesian(_x, _y, _z) sites_drawn.append([i, site, _nx, _ny, _nz, _x, _y, _z, xc, yc, zc]) return sites_drawn def draw_bonds(self, sites_drawn, max_distance = 2.75, bond_r = 0.5, R = 0.5, G = 0.5, B = 0.5, alpha = 0.5): cry = self.cry kwargs = { "radius": bond_r, "R": R, "G": G, "B": B, "alpha": alpha, "slices": self.nslices, "stacks": self.nstacks } # print() for i, site0 in enumerate(sites_drawn): i0, _site0, nx0, ny0, nz0, x0, y0, z0, xc0, yc0, zc0 = site0 for j, site1 in enumerate(sites_drawn): i1, _site1, nx1, ny1, nz1, x1, y1, z1, xc1, yc1, zc1 = site1 d2 = (xc1 - xc0)**2 + (yc1 - yc0)**2 + (zc1 - zc0)**2 d = sqrt(d2) if d > max_distance: continue # print(f"dis:", i, j, name_i, name_j, [xc0, yc0, zc0], [xc1, yc1, zc1], d) draw_cylinder(xc0, yc0, zc0, xc1, yc1, zc1, **kwargs) def draw_axes(self, X0, aij, sx, sy, sz, R_arrow, G_arrow, B_arrow, R_font, G_font, B_font): kwargs = { "radius": 0.15, "arrow_h": 0.5, "arrow_r": 0.3, "R": R_arrow, "G": G_arrow, "B": B_arrow, "alpha": 1.0, "slices": self.nslices, "stacks": self.nstacks } X0 = np.array(X0) Xa = X0 + aij[0] Xb = X0 + aij[1] Xc = X0 + aij[2] scale = [sx, sy, sz] fontcolor = [R_font, G_font, B_font] draw_arrow(*X0, *Xa, **kwargs) draw_arrow(*X0, *Xb, **kwargs) draw_arrow(*X0, *Xc, **kwargs) draw_text("a", *Xa, *scale, *fontcolor) draw_text("b", *Xb, *scale, *fontcolor) draw_text("c", *Xc, *scale, *fontcolor) def draw_crystal(self): cry = self.cry aij = cry.lattice_vectors() sites = cry.atom_site_list() ''' print("aij=", aij) for i, site in enumerate(sites): atom_name = site.atom_name_only() x, y, z = site.position(is_reduce01 = 1) print(f"#{i}: {atom_name} ({x}, {y}, {z}") ''' glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) self.init_projection(self.width, self.height) glMatrixMode(GL_MODELVIEW); glTranslatef(self.model_x, self.model_y, self.model_z) glScalef(self.scale, self.scale, self.scale) glRotatef(self.angle_x, 1.0, 0.0, 0.0) glRotatef(self.angle_y, 0.0, 1.0, 0.0) glRotatef(self.angle_z, 0.0, 0.0, 1.0) # glScalefで大きさの変わった法線ベクトルを正規化 glEnable(GL_NORMALIZE) X0 = [-3, -3, -3] scale = [2, 2, 2] arrowcolor = [0.5, 0.5, 0.5] fontcolor = [1, 0, 0] self.draw_axes(X0, aij, *scale, *arrowcolor, *fontcolor) self.draw_cell(aij, origin = [0, 0, 0], radius = 0.02, R = 0.7, G = 0.7, B = 0.7, alpha = self.alpha_cell) #球内のvectorを見せるため、vectorを先に描く sites_drawn = self.get_drawn_sites() # self.draw_vectors(sites_drawn, length = 1.5, radius = 0.15, arrow_h = 0.6, arrow_r = 0.3, # alpha = self.alpha_vector) self.draw_bonds(sites_drawn, max_distance = self.max_distance, bond_r = self.bond_r, R = 0.5, G = 0.5, B = 0.5, alpha = self.alpha_bond) self.draw_atoms(sites_drawn, alpha = self.alpha_atom) glutSwapBuffers() def menu_callback(self, value): print("value=", value) if value == 1: print("Option 1 selected") elif value == 2: print("Option 2 selected") elif value == 3: print("Option 3 selected") elif value == 4: glutLeaveMainLoop() return 0 def create_menu(self): self.menu = glutCreateMenu(self.menu_callback) glutAddMenuEntry("Option 1", 1) glutAddMenuEntry("Option 2", 2) glutAddMenuEntry("Option 3", 3) glutAddMenuEntry("Exit", 4) glutAttachMenu(GLUT_RIGHT_BUTTON) def draw_opengl(app, cfg): cry = read_cif(cfg.infile, app = app, print_level = 1) cry.PrintInf() nslices = cfg.get("nslices", 10) nstacks = cfg.get("nstacks", 10) alpha_cell = cfg.get("alpha_cell", 0.5) alpha_atom = cfg.get("alpha_atom", 0.5) alpha_bond = cfg.get("alpha_bond", 0.5) alpha_vector = cfg.get("alpha_vector", 0.8) r_mode = cfg.get("r_mode", "rion") k_rion = cfg.get("k_rion", 0.3) bond_r = cfg.get("bond_r", 0.1) max_distance = cfg.get("max_distance", 2.7) draw_range = [pfloat(s) for s in cfg.get("draw_range", "2,2,2").split(',')] latt = cry.lattice_parameters() aij = cry.lattice_vectors() print() print("Draw configuration:") print(f"draw_range: ", draw_range) print(f"max_distance: {max_distance}") print(f"Input lattice parameters: ", latt) center = (draw_range[0] * aij[0] + draw_range[1] * aij[1] + draw_range[2] * aij[2]) / 2.0 center += -3.0 print(f"Center of cells: ", center) cv = tkCrystalView(cry = cry, app = app, center = center, max_distance = max_distance, draw_range = draw_range, r_mode = r_mode, k_rion = k_rion, bond_r = bond_r, alpha_cell = alpha_cell, alpha_atom = alpha_atom, alpha_bond = alpha_bond, alpha_vector = alpha_vector, color_bg = [1.0, 1.0, 1.0, 1.0], nslices = nslices, nstacks = nstacks) cv.create_menu() glutDisplayFunc(cv.draw_crystal) glutReshapeFunc(cv.reshape) glutMouseFunc(cv.mouse) glutMotionFunc(cv.motion) glutKeyboardFunc(cv.keyboard) glutMainLoop() def key_pressed(key, x, y): global animation_running if key == b' ': # スペースキーを押すとアニメーションの開始/停止を切り替える animation_running = not animation_running def start_animation(): global animation_running animation_running = True def stop_animation(): global animation_running animation_running = False def draw(app, cfg): root = Tk() root.title("Control Panel") start_button = Button(root, text="Start", command=start_animation) start_button.pack(side=LEFT) stop_button = Button(root, text="Stop", command=stop_animation) stop_button.pack(side=RIGHT) glut_thread = threading.Thread(target = lambda: draw_opengl(app, cfg)) glut_thread.start() root.mainloop() def main(): app = tkApplication() cfg = initialize(app) # cfg.base_dir = app.check_arg('--base_dir', defval = '.', vartype = 'str') #cfg.formulaからログファイル名を作り、console出力をredirectする cfg.logfile = app.replace_path(None, template = ["{dirname}", "{filebody}-out.txt"], ext_dict = {}) app.redirect(heading = f"Open logfile [{cfg.logfile}]", targets = ["stdout", cfg.logfile], mode = 'w') #起動時引数で与えられたパラメータをcfgに設定 print() print("Update parameters from command line arguments:") app.update_vars(cfg, apply_default = True) app.cfg.print_parameters() if cfg.mode == 'draw': draw(app, cfg) else: app.terminate(message = "\n", usage = app.usage, post_message = f"Error in main: Invalide mode [{cfg.mode}]", pause = True, ) app.terminate(pause = True) if __name__ == "__main__": main()