NaCl-crystal.py

From GeoMod

To run the code, once you have VPython installed, you can either.

1. Download this file: NaCl-crystal.py
2. or, copy and paste the following code to a file (named "NaCl-crystal.py" for example) using your favorite text editor.

from visual import *

class Na:
    '''sodium'''
    def __init__(self, pos=vector(0,0,0), scale=1.0):
        self.pos = pos
        self.scale = scale
        self.atomic_number = 11
        
        self.atom = sphere()
        self.atom.pos = pos
        self.atom.radius = 1.02  #Angstroms ionic radius from http://environmentalchemistry.com/yogi/periodic/ionicradius.html
        self.atom.color = (0.9,0.9,0.9)
        self.label = label(pos=pos,text="Na", height =40, opacity = 0, box=0)

    def move(self, newpos):
        self.atom.pos = newpos
        self.label.pos = newpos

class Cl:
    '''chlorine'''
    def __init__(self, pos=vector(0,0,0), scale=1.0):
        self.scale = scale
        self.atomic_number = 17
        
        self.atom = sphere()
        self.atom.pos = pos
        self.atom.radius = 1.81  #Angstroms ionic radius from http://environmentalchemistry.com/yogi/periodic/ionicradius.html
        self.atom.color = (0.7,1.0,0.7)
        self.atom.pos = pos

        self.label = label(pos=pos,text="Cl", height =40, opacity = 0, box=0)

    def move(self, newpos):
        self.atom.pos = newpos
        self.label.pos = newpos


class NaCl:
    def __init__(self, pos = vector(0,0,0), scale=1.0):
        self.pos = pos
        self.scale = scale

##        self.Na = Na()
##        self.Cl = Cl()
##        
##        self.Na.atom.pos.x = (self.Na.atom.radius + self.Cl.atom.radius)/2.0
##        self.Cl.atom.pos.x = -(self.Na.atom.radius + self.Cl.atom.radius)/2.0

        self.atoms = {}
        self.atoms['Na'] = Na()
        self.atoms['Cl'] = Cl()
        
        #arrange atoms
        self.r = (self.atoms['Na'].atom.radius + self.atoms['Cl'].atom.radius) #/2.0
        self.atoms['Na'].move( self.pos + vector(self.r , 0, 0))
        self.atoms['Cl'].move( self.pos)

        self.atoms['Na'].label.visible = 0
        self.atoms['Cl'].label.visible = 0
        
    def rotate(self, angle, axis = vector(0,0,1)):
        for element,atom in self.atoms.iteritems():
            lpos = atom.atom.pos - self.pos
            lpos = rotate(lpos, angle=angle, axis=axis)
            atom.move( self.pos + lpos ) 
        
    def move(self, newpos):
        for element,atom in self.atoms.iteritems():
            lpos = atom.atom.pos - self.pos
            atom.move( lpos + newpos ) 
        

scene.background = (1,1,1)
scene.width=704 #+4+4
scene.height=576 #+24+4
scene.x = 0
scene.y = 0
scene.title = "Halite"

# arrange halite crystal by layers
s1 = NaCl()
s2 = NaCl(pos=vector(s1.r,0,s1.r))
s2.rotate(pi)
s3 = NaCl(pos=vector(s1.r,0,-s1.r))
s3.rotate(pi)
s4 = NaCl(pos=vector(s1.r*2.0,0,0))
s4.rotate(angle=pi/2.0, axis=vector(0,1,0))
s5 = NaCl(pos=vector(s1.r*2.0,0,0))
s5.rotate(angle=-pi/2.0, axis=vector(0,1,0))

s1u = NaCl(vector(s1.r,s1.r,0))
s1u.rotate(pi)
s2u = NaCl(vector(0,s1.r,s1.r))
#s2u.rotate(pi)
s3u = NaCl(vector(0,s1.r,-s1.r))
#s3u.rotate(pi)
s4u = NaCl(vector(s1.r*2.0,s1.r,-s1.r))
s4u.rotate(angle=-pi/2.0, axis=vector(0,1,0))
s5u = NaCl(vector(s1.r*2.0,s1.r,s1.r))
s5u.rotate(angle=pi/2.0, axis=vector(0,1,0))

s1d = NaCl(vector(s1.r,-s1.r,0))
s1d.rotate(pi)
s2d = NaCl(vector(0,-s1.r,s1.r))
#s2u.rotate(pi)
s3d = NaCl(vector(0,-s1.r,-s1.r))
#s3u.rotate(pi)
s4d = NaCl(vector(s1.r*2.0,-s1.r,-s1.r))
s4d.rotate(angle=-pi/2.0, axis=vector(0,1,0))
s5d = NaCl(vector(s1.r*2.0,-s1.r,s1.r))
s5d.rotate(angle=pi/2.0, axis=vector(0,1,0))

scene.center = (s1.r,0,0)

while 1:
    rate(50)
    scene.forward = rotate(scene.forward, angle=pi/360.0, axis=(0,1,0))