# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html # Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories # # Copyright (2005) Sandia Corporation. Under the terms of Contract # DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains # certain rights in this software. This software is distributed under # the GNU General Public License. # for python3 compatibility from __future__ import print_function # gl tool oneline = "3d interactive visualization via OpenGL" docstr = """ g = gl(d) create OpenGL display for data in d d = atom snapshot object (dump, data) g.bg("black") set background color (def = "black") g.size(N) set image size to NxN g.size(N,M) set image size to NxM g.rotate(60,135) view from z theta and azimuthal phi (def = 60,30) g.shift(x,y) translate by x,y pixels in view window (def = 0,0) g.zoom(0.5) scale image by factor (def = 1) g.box(0/1/2) 0/1/2 = none/variable/fixed box g.box(0/1/2,"green") set box color g.box(0/1/2,"red",4) set box edge thickness g.file = "image" file prefix for created images (def = "image") g.show(N) show image of snapshot at timestep N g.all() make images of all selected snapshots g.all(P) images of all, start file label at P g.all(N,M,P) make M images of snapshot N, start label at P g.pan(60,135,1.0,40,135,1.5) pan during all() operation g.pan() no pan during all() (default) args = z theta, azimuthal phi, zoom factor at beginning and end values at each step are interpolated between beginning and end values g.select = "$x > %g*3.0" string to pass to d.aselect.test() during all() g.select = "" no extra aselect (default) %g varies from 0.0 to 1.0 from beginning to end of all() g.acol(2,"green") set atom colors by atom type (1-N) g.acol([2,4],["red","blue"]) 1st arg = one type or list of types g.acol(0,"blue") 2nd arg = one color or list of colors g.acol(range(20),["red","blue"]) if list lengths unequal, interpolate g.acol(range(10),"loop") assign colors in loop, randomly ordered if 1st arg is 0, set all types to 2nd arg if list of types has a 0 (e.g. range(10)), +1 is added to each value interpolate means colors blend smoothly from one value to the next g.arad([1,2],[0.5,0.3]) set atom radii, same rules as acol() g.bcol() set bond color, same args as acol() g.brad() set bond thickness, same args as arad() g.tcol() set triangle color, same args as acol() g.tfill() set triangle fill, 0 fill, 1 line, 2 both g.lcol() set line color, same args as acol() g.lrad() set line thickness, same args as arad() g.adef() set atom/bond/tri/line properties to default g.bdef() default = "loop" for colors, 0.45 for radii g.tdef() default = 0.25 for bond/line thickness g.ldef() default = 0 fill by default 100 types are assigned if atom/bond/tri/line has type > # defined properties, is an error from vizinfo import colors access color list print(colors) list defined color names and RGB values colors["nickname"] = [R,G,B] set new RGB values from 0 to 255 140 pre-defined colors: red, green, blue, purple, yellow, black, white, etc Settings specific to gl tool: g.q(10) set quality of image (def = 5) g.axis(0/1) turn xyz axes off/on g.ortho(0/1) perspective (0) vs orthographic (1) view g.clip('xlo',0.25) clip in xyz from lo/hi at box fraction (0-1) g.reload() force all data to be reloaded g.cache = 0/1 turn off/on GL cache lists (def = on) theta,phi,x,y,scale,up = g.gview() grab all current view parameters g.sview(theta,phi,x,y,scale,up) set all view parameters data reload is necessary if dump selection is used to change the data cache lists usually improve graphics performance gview returns values to use in other commands: theta,phi are args to rotate() x,y are args to shift() scale is arg to zoom() up is a 3-vector arg to sview() """ # History # 9/05, Steve Plimpton (SNL): original version # ToDo list # when do aselect with select str while looping N times on same timestep # would not let you grow # of atoms selected # Variables # ztheta = vertical angle from z-azis of viewpoint # azphi = azimuthal angle of viewpoint # xshift,yshift = xy translation of scene (in pixels) # distance = size of simulation box (largest dim) # eye = viewpoint distance from center of scene # file = filename prefix to use for images produced # boxflag = 0/1/2 for drawing simulation box: none/variable/fixed # bxcol = color of box # bxthick = thickness of box lines # bgcol = color of background # vizinfo = scene attributes # center[3] = center point of simulation box # view[3] = direction towards eye in simulation box (unit vector) # up[3] = screen up direction in simulation box (unit vector) # right[3] = screen right direction in simulation box (unit vector) # Imports and external programs from math import sin,cos,sqrt,pi,acos from OpenGL.Tk import * from OpenGL.GLUT import * import Image from vizinfo import vizinfo # Class definition class gl: # -------------------------------------------------------------------- def __init__(self,data): self.data = data self.root = None self.xpixels = 512 self.ypixels = 512 self.ztheta = 60 self.azphi = 30 self.scale = 1.0 self.xshift = self.yshift = 0 self.file = "image" self.boxflag = 0 self.bxcol = [1,1,0] self.bxthick = 0.3 self.bgcol = [0,0,0] self.labels = [] self.panflag = 0 self.select = "" self.axisflag = 0 self.orthoflag = 1 self.nslices = 5 self.nstacks = 5 self.nsides = 10 self.theta_amplify = 2 self.shiny = 2 self.clipflag = 0 self.clipxlo = self.clipylo = self.clipzlo = 0.0 self.clipxhi = self.clipyhi = self.clipzhi = 1.0 self.nclist = 0 self.calllist = [0] # indexed by 1-Ntype, so start with 0 index self.cache = 1 self.cachelist = 0 self.boxdraw = [] self.atomdraw = [] self.bonddraw = [] self.tridraw = [] self.linedraw = [] self.ready = 0 self.create_window() self.vizinfo = vizinfo() self.adef() self.bdef() self.tdef() self.ldef() self.center = 3*[0] self.view = 3*[0] self.up = 3*[0] self.right = 3*[0] self.viewupright() # -------------------------------------------------------------------- def bg(self,color): from vizinfo import colors self.bgcol = [colors[color][0]/255.0,colors[color][1]/255.0, colors[color][2]/255.0] self.w.tkRedraw() # -------------------------------------------------------------------- def size(self,xnew,ynew=None): self.xpixels = xnew if not ynew: self.ypixels = self.xpixels else: self.ypixels = ynew self.create_window() # -------------------------------------------------------------------- def axis(self,value): self.axisflag = value self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def create_window(self): if self.root: self.root.destroy() from __main__ import tkroot self.root = Toplevel(tkroot) self.root.title('Pizza.py gl tool') self.w = MyOpengl(self.root,width=self.xpixels,height=self.ypixels, double=1,depth=1) self.w.pack(expand=YES) # self.w.pack(expand=YES,fill=BOTH) glViewport(0,0,self.xpixels,self.ypixels) glEnable(GL_LIGHTING); glEnable(GL_LIGHT0); glEnable(GL_DEPTH_TEST); glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,GL_TRUE); glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) self.rtrack = self.xpixels if self.ypixels > self.xpixels: self.rtrack = self.ypixels self.w.redraw = self.redraw self.w.parent = self self.w.tkRedraw() tkroot.update_idletasks() # force window to appear # -------------------------------------------------------------------- def clip(self,which,value): if which == "xlo": self.clipxlo = value if value > self.clipxhi: self.clipxlo = self.clipxhi elif which == "xhi": self.clipxhi = value if value < self.clipxlo: self.clipxhi = self.clipxlo elif which == "ylo": self.clipylo = value if value > self.clipyhi: self.clipylo = self.clipyhi elif which == "yhi": self.clipyhi = value if value < self.clipylo: self.clipyhi = self.clipylo elif which == "zlo": self.clipzlo = value if value > self.clipzhi: self.clipzlo = self.clipzhi elif which == "zhi": self.clipzhi = value if value < self.clipzlo: self.clipzhi = self.clipzlo oldflag = self.clipflag if self.clipxlo > 0 or self.clipylo > 0 or self.clipzlo > 0 or \ self.clipxhi < 1 or self.clipyhi < 1 or self.clipzhi < 1: self.clipflag = 1 else: self.clipflag = 0 if oldflag == 0 and self.clipflag == 0: return self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def q(self,value): self.nslices = value self.nstacks = value self.make_atom_calllist() self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def ortho(self,value): self.orthoflag = value self.w.tkRedraw() # -------------------------------------------------------------------- # set unit vectors for view,up,right from ztheta,azphi # assume +z in scene should be up on screen (unless looking down z-axis) # right = up x view def viewupright(self): self.view[0] = cos(pi*self.azphi/180) * sin(pi*self.ztheta/180) self.view[1] = sin(pi*self.azphi/180) * sin(pi*self.ztheta/180) self.view[2] = cos(pi*self.ztheta/180) if self.ztheta == 0.0: self.up[0] = cos(pi*self.azphi/180) self.up[1] = -sin(pi*self.azphi/180) self.up[2] = 0.0 elif self.ztheta == 180.0: self.up[0] = cos(pi*self.azphi/180) self.up[1] = sin(pi*self.azphi/180) self.up[2] = 0.0 else: dot = self.view[2] # dot = (0,0,1) . view self.up[0] = -dot*self.view[0] # up projected onto v = dot * v self.up[1] = -dot*self.view[1] # up perp to v = up - dot * v self.up[2] = 1.0 - dot*self.view[2] self.up = vecnorm(self.up) self.right = veccross(self.up,self.view) # -------------------------------------------------------------------- # reset ztheta,azphi and thus view,up.right # called as function from Pizza.py def rotate(self,ztheta,azphi): self.ztheta = ztheta self.azphi = azphi self.viewupright() self.setview() self.w.tkRedraw() # -------------------------------------------------------------------- # return all view params to reproduce current display via sview() def gview(self): return self.ztheta,self.azphi,self.xshift,self.yshift,self.scale,self.up # -------------------------------------------------------------------- # set current view, called by user with full set of view params # up is not settable via any other call, all other params are def sview(self,ztheta,azphi,xshift,yshift,scale,up): self.ztheta = ztheta self.azphi = azphi self.xshift = xshift self.yshift = yshift self.scale = scale self.up[0] = up[0] self.up[1] = up[1] self.up[2] = up[2] self.up = vecnorm(self.up) self.view[0] = cos(pi*self.azphi/180) * sin(pi*self.ztheta/180) self.view[1] = sin(pi*self.azphi/180) * sin(pi*self.ztheta/180) self.view[2] = cos(pi*self.ztheta/180) self.right = veccross(self.up,self.view) self.setview() self.w.tkRedraw() # -------------------------------------------------------------------- # rotation triggered by mouse trackball # project old,new onto unit trackball surf # rotate view,up around axis of rotation = old x new # right = up x view # reset ztheta,azphi from view def mouse_rotate(self,xnew,ynew,xold,yold): # change y pixels to measure from bottom of window instead of top yold = self.ypixels - yold ynew = self.ypixels - ynew # vold = unit vector to (xold,yold) projected onto trackball # vnew = unit vector to (xnew,ynew) projected onto trackball # return (no rotation) if either projection point is outside rtrack vold = [0,0,0] vold[0] = xold - (0.5*self.xpixels + self.xshift) vold[1] = yold - (0.5*self.ypixels + self.yshift) vold[2] = self.rtrack*self.rtrack - vold[0]*vold[0] - vold[1]*vold[1] if vold[2] < 0: return vold[2] = sqrt(vold[2]) vold = vecnorm(vold) vnew = [0,0,0] vnew[0] = xnew - (0.5*self.xpixels + self.xshift) vnew[1] = ynew - (0.5*self.ypixels + self.yshift) vnew[2] = self.rtrack*self.rtrack - vnew[0]*vnew[0] - vnew[1]*vnew[1] if vnew[2] < 0: return vnew[2] = sqrt(vnew[2]) vnew = vecnorm(vnew) # rot = trackball rotation axis in screen ref frame = vold x vnew # theta = angle of rotation = sin(theta) for small theta # axis = rotation axis in body ref frame described by right,up,view rot = veccross(vold,vnew) theta = sqrt(rot[0]*rot[0] + rot[1]*rot[1] + rot[2]*rot[2]) theta *= self.theta_amplify axis = [0,0,0] axis[0] = rot[0]*self.right[0] + rot[1]*self.up[0] + rot[2]*self.view[0] axis[1] = rot[0]*self.right[1] + rot[1]*self.up[1] + rot[2]*self.view[1] axis[2] = rot[0]*self.right[2] + rot[1]*self.up[2] + rot[2]*self.view[2] axis = vecnorm(axis) # view is changed by (axis x view) scaled by theta # up is changed by (axis x up) scaled by theta # force up to be perp to view via up_perp = up - (up . view) view # right = up x view delta = veccross(axis,self.view) self.view[0] -= theta*delta[0] self.view[1] -= theta*delta[1] self.view[2] -= theta*delta[2] self.view = vecnorm(self.view) delta = veccross(axis,self.up) self.up[0] -= theta*delta[0] self.up[1] -= theta*delta[1] self.up[2] -= theta*delta[2] dot = vecdot(self.up,self.view) self.up[0] -= dot*self.view[0] self.up[1] -= dot*self.view[1] self.up[2] -= dot*self.view[2] self.up = vecnorm(self.up) self.right = veccross(self.up,self.view) # convert new view to ztheta,azphi self.ztheta = acos(self.view[2])/pi * 180.0 if (self.ztheta == 0.0): self.azphi = 0.0 else: self.azphi = acos(self.view[0]/sin(pi*self.ztheta/180.0))/pi * 180.0 if self.view[1] < 0: self.azphi = 360.0 - self.azphi self.setview() self.w.tkRedraw() # -------------------------------------------------------------------- def shift(self,x,y): self.xshift = x; self.yshift = y; self.setview() self.w.tkRedraw() # -------------------------------------------------------------------- def zoom(self,scale): self.scale = scale self.setview() self.w.tkRedraw() # -------------------------------------------------------------------- # set view params needed by redraw # input: center = center of box # distance = size of scene (longest box length) # scale = zoom factor (1.0 = no zoom) # xshift,yshift = translation factor in pixels # view = unit vector from center to viewpoint # up = unit vector in up direction in scene # right = unit vector in right direction in scene # output: eye = distance to view scene from # xto,yto,zto = point to look to # xfrom,yfrom,zfrom = point to look from def setview(self): if not self.ready: return # no distance since no scene yet self.eye = 3 * self.distance / self.scale xfactor = 0.5*self.eye*self.xshift/self.xpixels yfactor = 0.5*self.eye*self.yshift/self.ypixels self.xto = self.center[0] - xfactor*self.right[0] - yfactor*self.up[0] self.yto = self.center[1] - xfactor*self.right[1] - yfactor*self.up[1] self.zto = self.center[2] - xfactor*self.right[2] - yfactor*self.up[2] self.xfrom = self.xto + self.eye*self.view[0] self.yfrom = self.yto + self.eye*self.view[1] self.zfrom = self.zto + self.eye*self.view[2] # -------------------------------------------------------------------- # box attributes, also used for triangle lines def box(self,*args): self.boxflag = args[0] if len(args) > 1: from vizinfo import colors self.bxcol = [colors[args[1]][0]/255.0,colors[args[1]][1]/255.0, colors[args[1]][2]/255.0] if len(args) > 2: self.bxthick = args[2] self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- # grab all selected snapshots from data object # add GL-specific info to each bond def reload(self): print("Loading data into gl tool ...") data = self.data self.timeframes = [] self.boxframes = [] self.atomframes = [] self.bondframes = [] self.triframes = [] self.lineframes = [] box = [] if self.boxflag == 2: box = data.maxbox() flag = 0 while 1: which,time,flag = data.iterator(flag) if flag == -1: break time,boxone,atoms,bonds,tris,lines = data.viz(which) if self.boxflag < 2: box = boxone if bonds: self.bonds_augment(bonds) self.timeframes.append(time) self.boxframes.append(box) self.atomframes.append(atoms) self.bondframes.append(bonds) self.triframes.append(tris) self.lineframes.append(lines) print(time,end='') sys.stdout.flush() print() self.nframes = len(self.timeframes) self.distance = compute_distance(self.boxframes[0]) self.center = compute_center(self.boxframes[0]) self.ready = 1 self.setview() # -------------------------------------------------------------------- def nolabel(self): self.cachelist = -self.cachelist self.labels = [] # -------------------------------------------------------------------- # show a single snapshot # distance from snapshot box or max box for all selected steps def show(self,ntime): data = self.data which = data.findtime(ntime) time,box,atoms,bonds,tris,lines = data.viz(which) if self.boxflag == 2: box = data.maxbox() self.distance = compute_distance(box) self.center = compute_center(box) if bonds: self.bonds_augment(bonds) self.boxdraw = box self.atomdraw = atoms self.bonddraw = bonds self.tridraw = tris self.linedraw = lines self.ready = 1 self.setview() self.cachelist = -self.cachelist self.w.tkRedraw() self.save() # -------------------------------------------------------------------- def pan(self,*list): if len(list) == 0: self.panflag = 0 else: self.panflag = 1 self.ztheta_start = list[0] self.azphi_start = list[1] self.scale_start = list[2] self.ztheta_stop = list[3] self.azphi_stop = list[4] self.scale_stop = list[5] # -------------------------------------------------------------------- def all(self,*list): data = self.data if len(list) == 0: nstart = 0 ncount = data.nselect elif len(list) == 1: nstart = list[0] ncount = data.nselect else: ntime = list[0] nstart = list[2] ncount = list[1] if self.boxflag == 2: box = data.maxbox() # loop over all selected steps # distance from 1st snapshot box or max box for all selected steps # recompute box center on 1st step or if panning if len(list) <= 1: n = nstart i = flag = 0 while 1: which,time,flag = data.iterator(flag) if flag == -1: break fraction = float(i) / (ncount-1) if self.select != "": newstr = self.select % fraction data.aselect.test(newstr,time) time,boxone,atoms,bonds,tris,lines = data.viz(which) if self.boxflag < 2: box = boxone if n == nstart: self.distance = compute_distance(box) if n < 10: file = self.file + "000" + str(n) elif n < 100: file = self.file + "00" + str(n) elif n < 1000: file = self.file + "0" + str(n) else: file = self.file + str(n) if self.panflag: self.ztheta = self.ztheta_start + \ fraction*(self.ztheta_stop - self.ztheta_start) self.azphi = self.azphi_start + \ fraction*(self.azphi_stop - self.azphi_start) self.scale = self.scale_start + \ fraction*(self.scale_stop - self.scale_start) self.viewupright() if n == nstart or self.panflag: self.center = compute_center(box) if bonds: self.bonds_augment(bonds) self.boxdraw = box self.atomdraw = atoms self.bonddraw = bonds self.tridraw = tris self.linedraw = lines self.ready = 1 self.setview() self.cachelist = -self.cachelist self.w.tkRedraw() self.save(file) print(time,end='') sys.stdout.flush() i += 1 n += 1 # loop ncount times on same step # distance from 1st snapshot box or max box for all selected steps # recompute box center on 1st step or if panning else: which = data.findtime(ntime) n = nstart for i in range(ncount): fraction = float(i) / (ncount-1) if self.select != "": newstr = self.select % fraction data.aselect.test(newstr,ntime) time,boxone,atoms,bonds,tris,lines = data.viz(which) if self.boxflag < 2: box = boxone if n == nstart: self.distance = compute_distance(box) if n < 10: file = self.file + "000" + str(n) elif n < 100: file = self.file + "00" + str(n) elif n < 1000: file = self.file + "0" + str(n) else: file = self.file + str(n) if self.panflag: self.ztheta = self.ztheta_start + \ fraction*(self.ztheta_stop - self.ztheta_start) self.azphi = self.azphi_start + \ fraction*(self.azphi_stop - self.azphi_start) self.scale = self.scale_start + \ fraction*(self.scale_stop - self.scale_start) self.viewupright() if n == nstart or self.panflag: self.center = compute_center(box) if bonds: self.bonds_augment(bonds) self.boxdraw = box self.atomdraw = atoms self.bonddraw = bonds self.tridraw = tris self.linedraw = lines self.ready = 1 self.setview() self.cachelist = -self.cachelist self.w.tkRedraw() self.save(file) print(n,end='') sys.stdout.flush() n += 1 print("\n%d images" % ncount) # -------------------------------------------------------------------- def display(self,index): self.boxdraw = self.boxframes[index] self.atomdraw = self.atomframes[index] self.bonddraw = self.bondframes[index] self.tridraw = self.triframes[index] self.linedraw = self.lineframes[index] self.ready = 1 self.cachelist = -self.cachelist self.w.tkRedraw() return (self.timeframes[index],len(self.atomdraw)) # -------------------------------------------------------------------- # draw the GL scene def redraw(self,o): # clear window to background color glClearColor(self.bgcol[0],self.bgcol[1],self.bgcol[2],0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) # not ready if no scene yet if not self.ready: return # set view from eye, distance, 3 lookat vectors (from,to,up) glMatrixMode(GL_PROJECTION) glLoadIdentity() if self.orthoflag: glOrtho(-0.25*self.eye,0.25*self.eye,-0.25*self.eye,0.25*self.eye, self.eye-2*self.distance,self.eye+2*self.distance) else: gluPerspective(30.0,1.0,0.01,10000.0) glMatrixMode(GL_MODELVIEW) glLoadIdentity() gluLookAt(self.xfrom,self.yfrom,self.zfrom,self.xto,self.yto,self.zto, self.up[0],self.up[1],self.up[2]) # draw scene from display list if caching allowed and list hasn't changed # else redraw and store as new display list if caching allowed if self.cache and self.cachelist > 0: glCallList(self.cachelist); else: if self.cache: if self.cachelist < 0: glDeleteLists(-self.cachelist,1) self.cachelist = glGenLists(1) glNewList(self.cachelist,GL_COMPILE_AND_EXECUTE) # draw box, clip-box, xyz axes, lines glDisable(GL_LIGHTING) if self.boxflag: self.draw_box(0) if self.clipflag: self.draw_box(1) if self.axisflag: self.draw_axes() ncolor = self.vizinfo.nlcolor for line in self.linedraw: itype = int(line[1]) if itype > ncolor: raise StandardError("line type too big") red,green,blue = self.vizinfo.lcolor[itype] glColor3f(red,green,blue) thick = self.vizinfo.lrad[itype] glLineWidth(thick) glBegin(GL_LINES) glVertex3f(line[2],line[3],line[4]) glVertex3f(line[5],line[6],line[7]) glEnd() glEnable(GL_LIGHTING) # draw non-clipped scene = atoms, bonds, triangles # draw atoms as collection of points # cannot put PointSize inside glBegin # so probably need to group atoms by type for best performance # or just allow one radius # need to scale radius appropriately with box size # or could leave it at absolute value # use POINT_SMOOTH to enable anti-aliasing and round points # multiple timesteps via vcr::play() is still not fast # caching makes it fast for single frame, but multiple frames is slow # need to enable clipping # if not self.clipflag: # glDisable(GL_LIGHTING) # glEnable(GL_POINT_SMOOTH) # glPointSize(self.vizinfo.arad[int(self.atomdraw[0][1])]) # glBegin(GL_POINTS) # for atom in self.atomdraw: # red,green,blue = self.vizinfo.acolor[int(atom[1])] # glColor(red,green,blue) # glVertex3d(atom[2],atom[3],atom[4]) # glEnd() # glEnable(GL_LIGHTING) if not self.clipflag: for atom in self.atomdraw: glTranslatef(atom[2],atom[3],atom[4]); glCallList(self.calllist[int(atom[1])]); glTranslatef(-atom[2],-atom[3],-atom[4]); if self.bonddraw: bound = 0.25 * self.distance ncolor = self.vizinfo.nbcolor for bond in self.bonddraw: if bond[10] > bound: continue itype = int(bond[1]) if itype > ncolor: raise StandardError("bond type too big") red,green,blue = self.vizinfo.bcolor[itype] rad = self.vizinfo.brad[itype] glPushMatrix() glTranslatef(bond[2],bond[3],bond[4]) glRotatef(bond[11],bond[12],bond[13],0.0) glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); obj = gluNewQuadric() gluCylinder(obj,rad,rad,bond[10],self.nsides,self.nsides) glPopMatrix() if self.tridraw: fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] if fillflag != 1: if fillflag: glEnable(GL_POLYGON_OFFSET_FILL) glPolygonOffset(1.0,1.0) glBegin(GL_TRIANGLES) ncolor = self.vizinfo.ntcolor for tri in self.tridraw: itype = int(tri[1]) if itype > ncolor: raise Exception("tri type too big") red,green,blue = self.vizinfo.tcolor[itype] glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); glNormal3f(tri[11],tri[12],tri[13]) glVertex3f(tri[2],tri[3],tri[4]) glVertex3f(tri[5],tri[6],tri[7]) glVertex3f(tri[8],tri[9],tri[10]) glEnd() if fillflag: glDisable(GL_POLYGON_OFFSET_FILL) if fillflag: glDisable(GL_LIGHTING) glPolygonMode(GL_FRONT_AND_BACK,GL_LINE) glLineWidth(self.bxthick) glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) glBegin(GL_TRIANGLES) for tri in self.tridraw: glVertex3f(tri[2],tri[3],tri[4]) glVertex3f(tri[5],tri[6],tri[7]) glVertex3f(tri[8],tri[9],tri[10]) glEnd() glEnable(GL_LIGHTING) glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) # draw clipped scene = atoms, bonds, triangles else: box = self.boxdraw xlo = box[0] + self.clipxlo*(box[3] - box[0]) xhi = box[0] + self.clipxhi*(box[3] - box[0]) ylo = box[1] + self.clipylo*(box[4] - box[1]) yhi = box[1] + self.clipyhi*(box[4] - box[1]) zlo = box[2] + self.clipzlo*(box[5] - box[2]) zhi = box[2] + self.clipzhi*(box[5] - box[2]) for atom in self.atomdraw: x,y,z = atom[2],atom[3],atom[4] if x >= xlo and x <= xhi and y >= ylo and y <= yhi and \ z >= zlo and z <= zhi: glTranslatef(x,y,z); glCallList(self.calllist[int(atom[1])]); glTranslatef(-x,-y,-z); if self.bonddraw: bound = 0.25 * self.distance ncolor = self.vizinfo.nbcolor for bond in self.bonddraw: xmin = min2(bond[2],bond[5]) xmax = max2(bond[2],bond[5]) ymin = min2(bond[3],bond[6]) ymax = max2(bond[3],bond[6]) zmin = min2(bond[4],bond[7]) zmax = max2(bond[4],bond[7]) if xmin >= xlo and xmax <= xhi and \ ymin >= ylo and ymax <= yhi and zmin >= zlo and zmax <= zhi: if bond[10] > bound: continue itype = int(bond[1]) if itype > ncolor: raise Exception("bond type too big") red,green,blue = self.vizinfo.bcolor[itype] rad = self.vizinfo.brad[itype] glPushMatrix() glTranslatef(bond[2],bond[3],bond[4]) glRotatef(bond[11],bond[12],bond[13],0.0) glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); obj = gluNewQuadric() gluCylinder(obj,rad,rad,bond[10],self.nsides,self.nsides) glPopMatrix() if self.tridraw: fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] if fillflag != 1: if fillflag: glEnable(GL_POLYGON_OFFSET_FILL) glPolygonOffset(1.0,1.0) glBegin(GL_TRIANGLES) ncolor = self.vizinfo.ntcolor for tri in self.tridraw: xmin = min3(tri[2],tri[5],tri[8]) xmax = max3(tri[2],tri[5],tri[8]) ymin = min3(tri[3],tri[6],tri[9]) ymax = max3(tri[3],tri[6],tri[9]) zmin = min3(tri[4],tri[7],tri[10]) zmax = max3(tri[4],tri[7],tri[10]) if xmin >= xlo and xmax <= xhi and \ ymin >= ylo and ymax <= yhi and \ zmin >= zlo and zmax <= zhi: itype = int(tri[1]) if itype > ncolor: raise Exception("tri type too big") red,green,blue = self.vizinfo.tcolor[itype] glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION, [red,green,blue,1.0]); glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); glNormal3f(tri[11],tri[12],tri[13]) glVertex3f(tri[2],tri[3],tri[4]) glVertex3f(tri[5],tri[6],tri[7]) glVertex3f(tri[8],tri[9],tri[10]) glEnd() if fillflag: glDisable(GL_POLYGON_OFFSET_FILL) if fillflag: glDisable(GL_LIGHTING) glPolygonMode(GL_FRONT_AND_BACK,GL_LINE) glLineWidth(self.bxthick) glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) glBegin(GL_TRIANGLES) for tri in self.tridraw: xmin = min3(tri[2],tri[5],tri[8]) xmax = max3(tri[2],tri[5],tri[8]) ymin = min3(tri[3],tri[6],tri[9]) ymax = max3(tri[3],tri[6],tri[9]) zmin = min3(tri[4],tri[7],tri[10]) zmax = max3(tri[4],tri[7],tri[10]) if xmin >= xlo and xmax <= xhi and \ ymin >= ylo and ymax <= yhi and \ zmin >= zlo and zmax <= zhi: glVertex3f(tri[2],tri[3],tri[4]) glVertex3f(tri[5],tri[6],tri[7]) glVertex3f(tri[8],tri[9],tri[10]) glEnd() glEnable(GL_LIGHTING) glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) if self.cache: glEndList() glFlush() # -------------------------------------------------------------------- # make new call list for each atom type # called when atom color/rad/quality is changed def make_atom_calllist(self): # extend calllist array if necessary if self.vizinfo.nacolor > self.nclist: for i in range(self.vizinfo.nacolor-self.nclist): self.calllist.append(0) self.nclist = self.vizinfo.nacolor # create new calllist for each atom type for itype in xrange(1,self.vizinfo.nacolor+1): if self.calllist[itype]: glDeleteLists(self.calllist[itype],1) ilist = glGenLists(1) self.calllist[itype] = ilist glNewList(ilist,GL_COMPILE) red,green,blue = self.vizinfo.acolor[itype] rad = self.vizinfo.arad[itype] glColor3f(red,green,blue); # glPointSize(10.0*rad) # glBegin(GL_POINTS) # glVertex3f(0.0,0.0,0.0) # glEnd() glMaterialfv(GL_FRONT,GL_EMISSION,[red,green,blue,1.0]); glMaterialf(GL_FRONT,GL_SHININESS,self.shiny); glutSolidSphere(rad,self.nslices,self.nstacks) glEndList() # -------------------------------------------------------------------- # augment bond info returned by viz() with info needed for GL draw # info = length, theta, -dy, dx for bond orientation def bonds_augment(self,bonds): for bond in bonds: dx = bond[5] - bond[2] dy = bond[6] - bond[3] dz = bond[7] - bond[4] length = sqrt(dx*dx + dy*dy + dz*dz) dx /= length dy /= length dz /= length theta = acos(dz)*180.0/pi bond += [length,theta,-dy,dx] # -------------------------------------------------------------------- def draw_box(self,flag): xlo,ylo,zlo,xhi,yhi,zhi = self.boxdraw if flag: tmp = xlo + self.clipxlo*(xhi - xlo) xhi = xlo + self.clipxhi*(xhi - xlo) xlo = tmp tmp = ylo + self.clipylo*(yhi - ylo) yhi = ylo + self.clipyhi*(yhi - ylo) ylo = tmp tmp = zlo + self.clipzlo*(zhi - zlo) zhi = zlo + self.clipzhi*(zhi - zlo) zlo = tmp glLineWidth(self.bxthick) glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) glBegin(GL_LINE_LOOP) glVertex3f(xlo,ylo,zlo) glVertex3f(xhi,ylo,zlo) glVertex3f(xhi,yhi,zlo) glVertex3f(xlo,yhi,zlo) glEnd() glBegin(GL_LINE_LOOP) glVertex3f(xlo,ylo,zhi) glVertex3f(xhi,ylo,zhi) glVertex3f(xhi,yhi,zhi) glVertex3f(xlo,yhi,zhi) glEnd() glBegin(GL_LINES) glVertex3f(xlo,ylo,zlo) glVertex3f(xlo,ylo,zhi) glVertex3f(xhi,ylo,zlo) glVertex3f(xhi,ylo,zhi) glVertex3f(xhi,yhi,zlo) glVertex3f(xhi,yhi,zhi) glVertex3f(xlo,yhi,zlo) glVertex3f(xlo,yhi,zhi) glEnd() # -------------------------------------------------------------------- def draw_axes(self): xlo,ylo,zlo,xhi,yhi,zhi = self.boxdraw delta = xhi-xlo if yhi-ylo > delta: delta = yhi-ylo if zhi-zlo > delta: delta = zhi-zlo delta *= 0.1 glLineWidth(self.bxthick) glBegin(GL_LINES) glColor3f(1,0,0) glVertex3f(xlo-delta,ylo-delta,zlo-delta) glVertex3f(xhi-delta,ylo-delta,zlo-delta) glColor3f(0,1,0) glVertex3f(xlo-delta,ylo-delta,zlo-delta) glVertex3f(xlo-delta,yhi-delta,zlo-delta) glColor3f(0,0,1) glVertex3f(xlo-delta,ylo-delta,zlo-delta) glVertex3f(xlo-delta,ylo-delta,zhi-delta) glEnd() # -------------------------------------------------------------------- def save(self,file=None): self.w.update() # force image on screen to be current before saving it pstring = glReadPixels(0,0,self.xpixels,self.ypixels, GL_RGBA,GL_UNSIGNED_BYTE) snapshot = Image.fromstring("RGBA",(self.xpixels,self.ypixels),pstring) snapshot = snapshot.transpose(Image.FLIP_TOP_BOTTOM) if not file: file = self.file snapshot.save(file + ".png") # -------------------------------------------------------------------- def adef(self): self.vizinfo.setcolors("atom",range(100),"loop") self.vizinfo.setradii("atom",range(100),0.45) self.make_atom_calllist() self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def bdef(self): self.vizinfo.setcolors("bond",range(100),"loop") self.vizinfo.setradii("bond",range(100),0.25) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def tdef(self): self.vizinfo.setcolors("tri",range(100),"loop") self.vizinfo.setfills("tri",range(100),0) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def ldef(self): self.vizinfo.setcolors("line",range(100),"loop") self.vizinfo.setradii("line",range(100),0.25) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def acol(self,atypes,colors): self.vizinfo.setcolors("atom",atypes,colors) self.make_atom_calllist() self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def arad(self,atypes,radii): self.vizinfo.setradii("atom",atypes,radii) self.make_atom_calllist() self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def bcol(self,btypes,colors): self.vizinfo.setcolors("bond",btypes,colors) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def brad(self,btypes,radii): self.vizinfo.setradii("bond",btypes,radii) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def tcol(self,ttypes,colors): self.vizinfo.setcolors("tri",ttypes,colors) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def tfill(self,ttypes,flags): self.vizinfo.setfills("tri",ttypes,flags) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def lcol(self,ltypes,colors): self.vizinfo.setcolors("line",ltypes,colors) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def lrad(self,ltypes,radii): self.vizinfo.setradii("line",ltypes,radii) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- # derived class from Togl's Opengl # overwrite redraw, translate, rotate, scale methods # latter 3 are mouse-motion methods class MyOpengl(Opengl): def __init__(self, master, cnf={}, **kw): args = (self,master,cnf) Opengl.__init__(*args,**kw) Opengl.autospin_allowed = 0 # redraw Opengl scene # call parent redraw() method def tkRedraw(self,*dummy): if not self.initialised: return self.tk.call(self._w,'makecurrent') self.redraw(self) self.tk.call(self._w,'swapbuffers') # left button translate # access parent xshift/yshift and call parent trans() method def tkTranslate(self,event): dx = event.x - self.xmouse dy = event.y - self.ymouse x = self.parent.xshift + dx y = self.parent.yshift - dy self.parent.shift(x,y) self.tkRedraw() self.tkRecordMouse(event) # middle button trackball # call parent mouse_rotate() method def tkRotate(self,event): self.parent.mouse_rotate(event.x,event.y,self.xmouse,self.ymouse) self.tkRedraw() self.tkRecordMouse(event) # right button zoom # access parent scale and call parent zoom() method def tkScale(self,event): scale = 1 - 0.01 * (event.y - self.ymouse) if scale < 0.001: scale = 0.001 elif scale > 1000: scale = 1000 scale *= self.parent.scale self.parent.zoom(scale) self.tkRedraw() self.tkRecordMouse(event) # -------------------------------------------------------------------- # draw a line segment def segment(p1,p2): glVertex3f(p1[0],p1[1],p1[2]) glVertex3f(p2[0],p2[1],p2[2]) # -------------------------------------------------------------------- # normalize a 3-vector to unit length def vecnorm(v): length = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]) return [v[0]/length,v[1]/length,v[2]/length] # -------------------------------------------------------------------- # dot product of two 3-vectors def vecdot(v1,v2): return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] # -------------------------------------------------------------------- # cross product of two 3-vectors def veccross(v1,v2): v = [0,0,0] v[0] = v1[1]*v2[2] - v1[2]*v2[1] v[1] = v1[2]*v2[0] - v1[0]*v2[2] v[2] = v1[0]*v2[1] - v1[1]*v2[0] return v # -------------------------------------------------------------------- # return characteristic distance of simulation domain = max dimension def compute_distance(box): distance = box[3]-box[0] if box[4]-box[1] > distance: distance = box[4]-box[1] if box[5]-box[2] > distance: distance = box[5]-box[2] return distance # -------------------------------------------------------------------- # return center of box as 3 vector def compute_center(box): c = [0,0,0] c[0] = 0.5 * (box[0] + box[3]) c[1] = 0.5 * (box[1] + box[4]) c[2] = 0.5 * (box[2] + box[5]) return c # -------------------------------------------------------------------- # return min of 2 values def min2(a,b): if b < a: a = b return a # -------------------------------------------------------------------- # return max of 2 values def max2(a,b): if b > a: a = b return a # -------------------------------------------------------------------- # return min of 3 values def min3(a,b,c): if b < a: a = b if c < a: a = c return a # -------------------------------------------------------------------- # return max of 3 values def max3(a,b,c): if b > a: a = b if c > a: a = c return a