INF205 lab worksheet 2

Submit through Canvas by 21st February 2025.

The example use case that we will be working with is molecular modelling and simulation using the truncated-shifted Lennard-Jones potential model (see also here). This model is sometimes abbreviated as LJTS, for "Lennard-Jones truncated and shifted." In this worksheet, you will load positions and velocities of molecules as 3D coordinates from files and compute the energy of the system.

If you do this as group work (up to three group members are allowed), it is necessary to say for each problem on the worksheet what group member had the lead or main responsibility for that part, and every group member needs to be the lead on at least one of the problems 3, 4, or 5. Best submit a ZIP file with the code in a subdirectory and a PDF, odt, docx or similar document in the main directory. Do not include compiled files (object files *.o, or any executable binary files) in your submission; instead, preferably include a makefile so that your code can be compiled easily.

3. Molecule class

Implement a class called Molecule, where molecules are represented as point masses with three-dimensional position coordinates q[3] and velocities v[3]. Each molecule should have a unique integer identifier. Follow common coding conventions in C++ and reasonable design practice, and make the member variables (properties) of the class private.

4. File input

Load the data on the molecules in a system from files in XYZ format.

Write your code and include a makefile such that the binary is called readxyz and takes two command-line arguments: One for an XYZ file containing position coordinates, one for an XYZ file containing velocities. The molecule IDs should be given incrementally, starting from zero. You can ignore the first data item in each line (atom label such as 'H', 'C', or similar). If your code contains any manually allocated memory, take care of the deallocation as well, and also follow good practices in dealing with memory otherwise.

As a result, for example, it should be possible to execute:

./readxyz 5molecules-positions.xyz 5molecules-velocities.xyz

In that case, e.g., the second molecule should have the identifier 1, the position q = {2.0, 0.0, 0.0}, and the velocity v = {0.866025, 0.0, -0.866025}.

It should also be possible to pass only the first command-line argument, as follows:

./readxyz 5molecules-positions.xyz

In that case, the velocities should all be set to zero.

5. Energy calculation

Extend the code by functionality to compute and output to std::cout:

For the kinetic energy, we assume that each molecule has a mass of m = 1. It is then given by Ekin = Σi vi2/2, where i is an index running over all the N molecules in the system, 0≤i<N.

The potential energy is given as a sum over all pairs of molecules Epot = Σi Σj>i u(rij). Therein, rij is the distance between molecules i and j. The truncated-shifted Lennard-Jones potential u(rij) is given by:

There, uLJ is the Lennard-Jones potential, given by uLJ(r) = 4 · (1/r12 - 1/r6).

So the Lennard-Jones potential is being truncated at a cutoff distance of 2.5, and shifted by adding -uLJ(2.5) = +0.0163169. All pairs of molecules i and j with a distance below 2.5 contribute to the potential energy, while the others do not contribute; each pair is counted only once, i.e., if you include the interaction between i and j, do not include it again in reverse order, between j and i.

Follow an object oriented design approach and encapsulate the key parts of the functionality in classes and methods. For this purpose, introduce at least one more class in addition to Molecule. Also give some first thoughts to a numerically efficient implementation of the above. Explain them briefly in the PDF/docx (or other text file) that you submit together with your code as a ZIP archive. We will be paying more attention to computational efficency in the third worksheet.

As a result, for example, when running

./readxyz 5molecules-positions.xyz 5molecules-velocities.xyz

an output should be obtained such as follows:

E_kin = 3.75.
E_pot = -0.180826.
E_kin + E_pot = 3.56917.

The reason is that in the example files, each of the five molecules has the kinetic energy 0.75. Also, the adjacent molecule pairs 0-1, 1-2, 2-3, and 3-4 are each at a distance of exactly 2, while any other pairs are at a distance greater than the cutoff of 2.5. So the potential energy is Epot = 4u(2), where each u(2) = uLJ(2) - uLJ(2.5) = 4·(1/212 - 1/26) + 0.0163169 = -0.0452065.

Index