share/hoomd/examples/dihedral_table

   1 hoomd_script::dihedral::table
   2 # run this script with "python -x filename" to skip the first line, or remove this header
   3
   4 # ---- angle_table.py ----
   5 from hoomd_script import *
   6 import math
   7
   8 # generate the polymer system
   9 sys=init.read_xml('angle_table.xml')
  10
  11 bonds = bond.harmonic()
  12 bonds.bond_coeff.set('polymer', k=330, r0=1.0)
  13
  14 # force field setup
  15 def harmonic(theta, kappa, theta_0):
  16     V = 0.5 * kappa * (theta-theta_0)**2;
  17     T = -kappa*(theta-theta_0);
  18     return (V, T)
  19
  20 # initialize the angle table force
  21 dtable = dihedral.table(width=1000)
  22 dtable.dihedral_coeff.set('dihedral', func=harmonic, coeff=dict(kappa=330, theta_0=math.pi/2))
  23
  24 # NVT integration
  25 all = group.all()
  26 integrate.mode_standard(dt=0.005)
  27 integrate.bdnvt(group=all, T=1.0)
  28
  29 dump.mol2('dihedral_table.mol2');
  30
  31 dump.dcd(filename="trajectory.dcd", period=10, overwrite=True)
  32
  33 # warm up the system
  34 run(2000)
  35
  36 # A visualization of trajectory.dcd should show a polymer wiggling while maintaining roughly the r0 and theta_0 values
  37 # for bonds and angles