test/hoomd_script/test_dihedral_table.py

   1 # -*- coding: iso-8859-1 -*-
   2 # Maintainer: joaander
   3
   4 from hoomd_script import *
   5 import unittest
   6 import os
   7 import math
   8
   9 # tests dihedral.table
  10 class dihedral_table_tests (unittest.TestCase):
  11     def setUp(self):
  12         print
  13         # create a polymer system and add a dihedral to it
  14         self.polymer1 = dict(bond_len=1.2, type=['A']*4, bond="linear", count=10);
  15         self.polymers = [self.polymer1]
  16         self.box = data.boxdim(L=35);
  17         self.separation=dict(A=0.35, B=0.35)
  18         self.sys = init.create_random_polymers(box=self.box, polymers=self.polymers, separation=self.separation);
  19
  20         for i in range(len(self.sys.particles)//4-1):
  21             self.sys.dihedrals.add('dihedralA', 4*i+0, 4*i+1, 4*i+2, 4*i+3);
  22
  23         sorter.set_params(grid=8)
  24
  25     # test to see that se can create a dihedral.table
  26     def test_create(self):
  27         dihedral.table(width=100);
  28
  29     # test setting the table
  30     def test_set_coeff(self):
  31         def har(theta, kappa, theta_0):
  32             V = 0.5 * kappa * (theta-theta_0)**2;
  33             T = -kappa*(theta-theta_0);
  34             return (V, T)
  35
  36         harmonic = dihedral.table(width=1000)
  37         harmonic.dihedral_coeff.set('dihedralA', func=har, coeff=dict(kappa=1,theta_0=.1))
  38         all = group.all();
  39         integrate.mode_standard(dt=0.005);
  40         integrate.nve(all);
  41         run(100);
  42
  43     # test coefficient not set checking
  44     def test_set_coeff_fail(self):
  45         harmonic = dihedral.table(width=123)
  46         all = group.all();
  47         integrate.mode_standard(dt=0.005);
  48         integrate.nve(all);
  49         self.assertRaises(RuntimeError, run, 100);
  50
  51     # compare against harmonic dihedral
  52     def test_harmonic_compare(self):
  53         harmonic_1 = dihedral.table(width=1000)
  54         harmonic_1.dihedral_coeff.set('dihedralA', func=lambda theta: (0.5*1*( 1 + math.cos(theta)), 0.5*1*math.sin(theta)),coeff=dict())
  55         harmonic_2 = dihedral.harmonic()
  56         harmonic_2.set_coeff('dihedralA', k=1.0, d=1,n=1)
  57         integrate.mode_standard(dt=0.005);
  58         all = group.all()
  59         integrate.nve(all)
  60         run(1)
  61         for i in range(len(self.sys.particles)):
  62             f_1 = harmonic_1.forces[i]
  63             f_2 = harmonic_2.forces[i]
  64             # we have to have a very rough tolerance (~10%) because
  65             # of 1) discretization of the potential and 2) different handling of precision issues in both potentials
  66             self.assertAlmostEqual(f_1.energy, f_2.energy,3)
  67             self.assertAlmostEqual(f_1.force[0], f_2.force[0],2)
  68             self.assertAlmostEqual(f_1.force[1], f_2.force[1],2)
  69             self.assertAlmostEqual(f_1.force[2], f_2.force[2],2)
  70
  71     def tearDown(self):
  72         del self.sys
  73         init.reset();
  74
  75
  76 if __name__ == '__main__':
  77     unittest.main(argv = ['test.py', '-v'])