share/hoomd/examples/bond_table

   1 hoomd_script::bond::table
   2 # run this script with "python -x filename" to skip the first line, or remove this header
   3
   4 # ---- bond_table.py ----
   5 from hoomd_script import *
   6 import math
   7
   8 # parameters
   9 phi_P = 0.0001
  10 n_poly = 1
  11 T = 1.2
  12
  13 # This polymer is A-B-A with the two bonds having two different types.
  14 polymer = dict(bond_len=1.2, type=['A'] + ['B'] + ['A'], bond=[(0,1, 'bond1'), (1,2, 'bond2')], count=n_poly)
  15
  16
  17 # perform some simple math to find the length of the box
  18 N = len(polymer['type']) * polymer['count'];
  19 L = math.pow(math.pi * N / (6.0 * phi_P), 1.0/3.0);
  20
  21 # generate the polymer system
  22 sys=init.create_random_polymers(box=hoomd.BoxDim(L), polymers=[polymer], separation=dict(A=0.35, B=0.35));
  23
  24 # Assign the particles positions for simplicity
  25 p = sys.particles[0]
  26 p.position = (0,0,0);
  27 p = sys.particles[1]
  28 p.position = (0,1,0);
  29 p = sys.particles[2]
  30 p.position = (0,2,0)
  31
  32 # force field setup
  33 def harmonic(r, rmin, rmax, kappa, r0):
  34     V = 0.5 * kappa * (r-r0)**2;
  35     F = -kappa*(r-r0);
  36     return (V, F)
  37
  38
  39 # Initializing the bond Table forces.
  40 btable = bond.table(width=1000)
  41 # For the first bond it reads from the file bond_table.dat.  This table specifies a harmonic bond of kappa =2, r0=0.084
  42 btable.set_from_file('bond1', 'bond_table.dat')
  43 #  For the second bond we use the harmonic function defined above
  44 btable.bond_coeff.set('bond2', func=harmonic, rmin=0.2, rmax=5.0, coeff=dict(kappa=330, r0=0.84))
  45
  46
  47 # Apply LJ forces to the particles
  48 lj = pair.lj(r_cut=3.0)
  49 lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, alpha=0.0)
  50 lj.pair_coeff.set('A', 'B', epsilon=1.0, sigma=1.0, alpha=0.0)
  51 lj.pair_coeff.set('B', 'B', epsilon=1.0, sigma=1.0, alpha=1.0)
  52
  53 # NVT integration
  54 all = group.all()
  55 integrate.mode_standard(dt=0.005)
  56 integrate.nvt(group=all, T=T, tau=0.5)
  57
  58 xml = dump.xml();
  59 xml.set_params(all=True);
  60 xml.set_params(image=False);
  61 xml.write('Init.xml')
  62
  63 dump.dcd(filename="trajectory.dcd", period=10, overwrite=True)
  64
  65 # warm up the system
  66 run(2000)
  67
  68 # The user should observe a polymer oscillating up and down on the screen due to the difference in bond types.