share/hoomd/examples/init_create_empty

   1 hoomd_script::init::create_empty hoomd_script::pair::lj hoomd_script::integrate::mode_standard hoomd_script::integrate::nve
   2 # Due to deficiencies in doxygen, the commands used in this example are listed explicitly here
   3 # run this script with "python -x filename" to skip the first line, or remove this header
   4
   5 # ---- init_create_empty.py ----
   6 from __future__ import division
   7 from hoomd_script import *
   8 import random
   9
  10 # parameters of the simple cubic crystal
  11 a = 1.2
  12 m = 10
  13
  14 system = init.create_empty(N=m*m*m, box=data.boxdim(L=m*a), particle_types=['A', 'B'])
  15
  16 # initialize a simple cubic array of particles
  17 lo = - m*a / 2.0;
  18 for p in system.particles:
  19     (i, j, k) = (p.tag % m, p.tag//m % m, p.tag//m**2  % m)
  20     p.position = (lo + i*a + a/2, lo + j*a + a/2, lo + k*a + a/2)
  21
  22 # make the top half type B
  23 top = group.cuboid(name="top", ymin=0)
  24 for p in top:
  25     p.type = 'B'
  26
  27 # initialize the velocities to be a thermal distribution
  28 random.seed(1234);
  29 T = 1.0
  30 px = py = pz = 0.0;
  31 for p in system.particles:
  32     mass = p.mass;
  33     vx = random.gauss(0, T / mass)
  34     vy = random.gauss(0, T / mass)
  35     vz = random.gauss(0, T / mass)
  36
  37     p.velocity = (vx, vy, vz)
  38
  39     # sum the total system momentum
  40     px += mass*vx;
  41     py += mass*vy;
  42     pz += mass*vz;
  43
  44 # compute average momentum
  45 px /= m*m*m;
  46 py /= m*m*m;
  47 pz /= m*m*m;
  48
  49 # subtract that average momentum from each particle
  50 for p in system.particles:
  51     mass = p.mass;
  52     v = p.velocity;
  53     p.velocity = (v[0] - px/mass, v[1] - py/mass, v[2] - pz/mass);
  54
  55 # simple lennard jones potential
  56 lj = pair.lj(r_cut=3.0)
  57 lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0)
  58 lj.pair_coeff.set('A', 'B', epsilon=1.0, sigma=1.0)
  59 lj.pair_coeff.set('B', 'B', epsilon=1.0, sigma=1.0)
  60
  61 # integrate forward in the nve ensemble
  62 all = group.all()
  63 integrate.mode_standard(dt=0.005)
  64 integrate.nve(group=all)
  65
  66 # run for 20k steps
  67 run(20000)
  68