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