Changed directories names (src->pizza, test_dit->tests), removed old example dir.

# compute pairwise energy between 2 groups of atoms

# Syntax: group_energy.py data.file dump.file1 dump.file2 ...

# Author: Paul Crozier (Sandia)

# return distance sq between 2 atoms with PBC

def distance(box,x1,y1,z1,x2,y2,z2):

  delx = x2 - x1

  dely = y2 - y1

  delz = z2 - z1

  xprd = box[3] - box[0]

  yprd = box[4] - box[1]

  zprd = box[5] - box[2]

  if abs(delx) > 0.5*xprd:

    if delx < 0.0:

      delx += xprd

    else:

      delx -= xprd

  if abs(dely) > 0.5*yprd:

    if dely < 0.0:

      dely += yprd

    else:

      dely -= yprd

  if abs(delz) > 0.5*zprd:

    if delz < 0.0:

      delz += zprd

    else:

      delz -= zprd

  distsq = delx*delx + dely*dely + delz*delz

  return distsq

# main script

if len(argv) < 3:

  raise StandardError,"group_energy.py data.file dump.file1 dump.file2 ..."

dt = data(argv[1])				# data file

q = dt.get("Atoms",4)

files = ' '.join(argv[2:])		        # dump files

d = dump(files,0)

d.map(1,"id",2,"type",3,"x",4,"y",5,"z")

p = pair("lj/charmm/coul/charmm")

p.coeff(dt)

cut1 = 8.0                                      # potential cutoffs

cut2 = 10.0

cut3 = 8.0

cut4 = 10.0

p.init(cut1,cut2,cut3,cut4)

maxcut = cut1

if cut2 > maxcut: maxcut = cut2

if cut3 > maxcut: maxcut = cut3

if cut4 > maxcut: maxcut = cut4

maxcut_sq = maxcut*maxcut

while 1:

  time = d.next()

  if time < 0: break

  d.unscale(time)

  box = (d.snaps[0].xlo,d.snaps[0].ylo,d.snaps[0].zlo,

         d.snaps[0].xhi,d.snaps[0].yhi,d.snaps[0].zhi)

  d.aselect.all(time)

  d.aselect.test("$id >= 14306 and $id <= 14516",time)          # 1st group

  id1,type1,x1,y1,z1 = d.vecs(time,"id","type","x","y","z")

  d.aselect.all(time)                                           # 2nd group

  d.aselect.test("$id >= 1 and $id <= 7243 or $id >= 7274 and $id <= 14283",

                 time)

  id2,type2,x2,y2,z2 = d.vecs(time,"id","type","x","y","z")

  id1 = map(int,id1)

  id2 = map(int,id2)

  type1 = map(int,type1)

  type2 = map(int,type2)

  n1 = len(type1)

  n2 = len(type2)

  for i in xrange(n1): 

    id1[i] -= 1

    type1[i] -= 1

  for i in xrange(n2): 

    id2[i] -= 1

    type2[i] -= 1

  e_coul_sum = 0.0

  e_vdwl_sum = 0.0

  for i in xrange(n1):

    typei = type1[i]

    qi = q[id1[i]]

    for j in xrange(n2):

      rsq = distance(box,x1[i],y1[i],z1[i],x2[j],y2[j],z2[j])

      if rsq < maxcut_sq:

        eng_coul,eng_vdwl = p.single(rsq,typei,type2[j],qi,q[id2[j]])

        e_coul_sum += eng_coul     

        e_vdwl_sum += eng_vdwl    

  print "eng_coul = %g at timestep %d" % (e_coul_sum,time)

  print "eng_vdwl = %g at timestep %d" % (e_vdwl_sum,time)

  d.tselect.none()

  d.tselect.one(time)

  d.delete()

# movie of bucky-ball data

d = dump("dump.bucky")

d.set("$center = ((int($id)-1)/61+1)*61")

d.owrap("center")

r = raster(d)

r.file = "bucky"

r.acol([1,2,3],["blue","red","green"])

r.arad([1,2,3],[0.5,0.5,2.5])

r.box(1)

# spin

d.tselect.test("$t == 0")

r.pan(60,150,1,60,30,1)

r.all(0,25,0)

# slice

r.pan()

r.rotate(60,30)

r.select = "$x < -17 + (1-%g) * 34 and $type > 1"

r.all(0,25,25)

# timestep animation

d.tselect.all()

d.aselect.test("$type > 1")

r.rotate(60,30)

r.select = ""

r.all(50)

# movie of flow around obstacle

d = dump("dump.flow")

d.map(1,"id",2,"type",3,"x",4,"y",5,"z",6,"vx",7,"vy")

d.set("$ke = sqrt($vx*$vx + $vy*$vy)")

d.spread("vx",100,"color")

d.atype = "color"

r = raster(d)

r.acol(range(100),["red","red","red","red","yellow","green","blue","purple","purple","purple","purple"])

r.arad(range(100),0.5)

r.rotate(0,-90)

r.zoom(1.5)

r.file = "flow"

r.all()

# movie of melting LJ solid

a = dump("dump.melt")

a.tselect.test("$t == 0")

a.scale()

a.set("$ix = int($x * 4)")

a.set("$iy = int($y * 4)")

a.set("$iz = int($z * 4)")

a.set("$type = ($ix + $iy + $iz) % 2 + 1")

a.unscale()

a.tselect.all()

a.clone(0,"type")

r = raster(a)

r.acol([1,2],["red","green"])

r.arad([1,2],0.5)

r.file = "melt"

r.pan(130,25,1,60,135,0.6)

r.all()

# movie of self-assembling micelles

d = dump("dump.micelle")

s = svg(d)

s.acol([1,2,3,4],["blue","red","cyan","yellow"])

s.arad(range(4),0.5)

s.zoom(1.5)

s.file = "micelle"

s.all()