Lennard Jones reference configuration

In this example, an instantaneous configuration is compared against published values for the potential energy.

First, the potential energy between two particles is compared to the analytical value. Next, the energy of an entire configuration is compared against the NIST SRSW reference calculations.

import unittest
import feasst as fst

class TestMonteCarloLJRef(unittest.TestCase):
    """The FEASST implementation of the LJ potential is tested against known cases."""
    def test_two_particle(self):
        """Test the LJ potential against analytical calculation of two particles"""
        system = fst.System()
        system.add(fst.MakeConfiguration(fst.args({"cubic_box_length": "8",
                                                   "particle_type0": fst.install_dir() + "/forcefield/lj.fstprt"})))
        system.add(fst.MakePotential(fst.MakeLennardJones()))
        system.add(fst.MakePotential(fst.MakeLongRangeCorrections()))
        displacement = 1.2345
        system.get_configuration().add_particle_of_type(0)
        system.get_configuration().add_particle_of_type(0)
        system.get_configuration().update_positions(fst.Double2DVector([[0., 0., 0.], [displacement, 0., 0.]]))

        # compute the energy of the system
        system.energy()

        # compute the expected analytical LJ and LRC energies
        enlj = 4*(displacement**(-12) - displacement**(-6))
        rcut = system.configuration().model_params().select("cutoff").value(0)
        enlrc = (8./3.)*fst.PI*system.configuration().num_particles()**2/ \
            system.configuration().domain().volume()*((1./3.)*rcut**(-9) - rcut**(-3))

        # Compare the analytical results with the FEASST computed energies.
        # The energies of the individual potentials (e.g., LJ and LRC) are stored as profiles with
        # indices based on the order that the potentials were initialized.
        # Thus, profile index 0 refers to LJ while 1 refers to LRC.
        # In addition, the last computed value of the energy of all potentials is also stored.
        self.assertAlmostEqual(enlj, system.stored_energy_profile()[0], 15)
        self.assertAlmostEqual(enlrc, system.stored_energy_profile()[1], 15)
        self.assertAlmostEqual(enlj + enlrc, system.stored_energy(), 15)

    def test_srsw_ref_config(self):
        """Test the LJ potential against a configuration of 30 particles.
        In particular, the 4th configuration of the LJ SRSW reference:
        https://www.nist.gov/mml/csd/chemical-informatics-research-group/lennard-jones-fluid-reference-calculations
        """
        system = fst.System()
        system.add(fst.MakeConfiguration(fst.args({"cubic_box_length": "8",
                                                   "particle_type0": fst.install_dir() + "/forcefield/lj.fstprt",
                                                   "xyz_file": fst.install_dir() + "/plugin/configuration/test/data/lj_sample_config_periodic4.xyz"
                                                  })))
        system.add(fst.MakePotential(fst.MakeLennardJones()))
        system.add(fst.MakePotential(fst.MakeLongRangeCorrections()))
        self.assertEqual(30, system.configuration().num_particles())
        system.energy()
        enlj = -16.790321304625856
        enlrc = -0.5451660014945704
        self.assertAlmostEqual(enlj, system.stored_energy_profile()[0], 15)
        self.assertAlmostEqual(enlrc, system.stored_energy_profile()[1], 15)
        self.assertAlmostEqual(enlj + enlrc, system.energy(), 15)

    def test_srsw_ref_config_triclinic(self):
        """Test the LJ potential against a configuration of 300 particles in a trinclinic cell.
        In particular, the 3th configuration of the triclinic LJ SRSW reference:
        https://www.nist.gov/mml/csd/chemical-informatics-group/lennard-jones-fluid-reference-calculations-non-cuboid-cell
        """
        system = fst.System()
        system.add(fst.MakeConfiguration(fst.args({
            "side_length0": "10.0",
            "side_length1": "9.84807753012208",
            "side_length2": "9.64974312607518",
            "xy": "1.7364817766693041",
            "xz": "2.5881904510252074",
            "yz": "0.42863479791864567",
            "particle_type0": fst.install_dir() + "/forcefield/lj.fstprt",
            "xyz_file": fst.install_dir() + "/plugin/configuration/test/data/lj_triclinic_sample_config_periodic3.xyz"})))
        system.add(fst.MakePotential(fst.MakeLennardJones()))
        system.add(fst.MakePotential(fst.MakeLongRangeCorrections()))
        self.assertEqual(300, system.configuration().num_particles())
        system.energy()
        enlj = -505.78567945268367
        enlrc = -29.37186430697248
        self.assertAlmostEqual(enlj, system.stored_energy_profile()[0], 15)
        self.assertAlmostEqual(enlrc, system.stored_energy_profile()[1], 15)
        self.assertAlmostEqual(enlj + enlrc, system.energy(), 15)

Run the tests.

unittest.main(argv=[''], verbosity=2, exit=False)

test_srsw_ref_config (__main__.TestMonteCarloLJRef)
Test the LJ potential against a configuration of 30 particles. ... ok
test_srsw_ref_config_triclinic (__main__.TestMonteCarloLJRef)
Test the LJ potential against a configuration of 300 particles in a trinclinic cell. ... ok
test_two_particle (__main__.TestMonteCarloLJRef)
Test the LJ potential against analytical calculation of two particles ... ok

----------------------------------------------------------------------
Ran 3 tests in 0.024s

OK

<unittest.main.TestProgram at 0x7fb84452dcd0>

Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!