INF203 worksheet 3

Submit through Canvas by Monday, 16th June 2025, 17.30 CEST:

3.1 Configuration output in XYZ format

Implement functionality for file output of a configuration in XYZ format, which can be visualized with tools such as JMol. See configuration.xyz for an example of such a file from a simulation of the reference scenario.

Your code should create the following files:

The trajectory is a sequence of frames in XYZ format which are simply concatenated. (See trajectory.xyz example.) It can be used to generate video output. Do not include a frame after every MC step, but with a lower frequency, e.g., once every 200 steps.

3.2 Potential energy of a distorted box

Implement functionality by which the potential energy of the simulation box is calculated while all coordinates (including the size of the box) are subject to a distortion, i.e. all x, y, z coordinates and distances are scaled by a factor sx, sy, and sz, respectively. Determine the internal energy change ΔU = Epot(distorted) - Epot(undistorted).

Note that for this it is enough to rescale the intermolecular distance vectors when calculating the pair potential.

Also determine the surface area change ΔA = A(distorted) - A(undistorted), assuming that the interface is perpendicular to the y axis and that there are two interfaces in the system, so that A = 2LxLz.

3.3 Test area method

The test area method determines the surface tension γ from the expression ΔF = γΔA = -T ln ⟨ exp(-ΔU / T) ⟩. As usual, the angular brackets stand for the statistical average (over the sampled configurations).

Therefore, in order to implement this method, we need to sample the statistical average of the quantity exp(-ΔU / T), where ΔU is the potential energy change for a small distortion and T is the temperature.

Extend your code such that it determines the potential energy for the undistorted system once after every MC step, and additionally, also under two different distortion vectors that both conserve the volume (sx,1sy,1sz,1 = sx,2sy,2sz,2 = 1). Make it so that one of the distortions increases the surface area by a factor ζ, i.e., sx,1sz,1 = ζ, while the other distortion vector decreases the surface area by the same factor, sy,2 = ζ. That factor should be close to 1, for example ζ = 1.00001.

Pay attention to distinguish between an equilibration and a production stage in the simulation, and avoid including the initial equilibration phase - where configurations do not yet reflect the normal physical state of the system - in the sampling for the average value of exp(-ΔU / T), and any other quantities for which accumulated averages are taken, such as the potential energy of the system.

Do a separate sampling and determine separate results for each of the two distortion vectors: The one that increases the surface area and the one that decreases it.

3.4 Output of simulation results

Extend the console output from your code to include at least, for each of the two distortion vectors s1 and s2, the present value and the accumulated average of the sampled quantity exp(-ΔU / T). For each s, it should also include the value of the surface tension γ based on ⟨exp(-ΔU / T)⟩. It should also still continue the present value and the accumulated average of the potential energy Epot of the (undistorted) system.

Make sure that the console output is also stored in a file, either by redirecting the console output from the command line or by explicitly writing into a file from Python. The output should not be made after every MC step, but at some frequency, which can be higher than that for the XYZ output. (For example, leave this as you had it from no. 2.3 in worksheet 2.) In this way, your result files will not be very large and can be included in your submission.

3.5 Simulation run for the reference system

Try it out for the reference scenario with the vapour-liquid interface at T = 0.8. Run at least 10 000 steps for equilibration, i.e., values from these steps should not be counted for any quantities that you are computing, and at least 20 000 steps after that for production. (These are absolute minimum values - if it can be done, please do run your simulation much longer than this.)

Analyse your findings. Are they in agreement with Table 2 (p. 1515) in Vrabec et al. (2006), doi:10.1080/00268970600556774?

Note that it is perfectly OK if the values don't exactly agree, since our reference scenario is for a much smaller simulation volume than theirs (Vrabec et al. (2006) have N = 2872, whereas we have N = 162 …), and also the number of simulation steps will be comparably small, which limits the sampling. But if they are completely far off, it most likely points to a problem.

Index