Example 03: Manual calculation of Cu melting temperature#
In this example, the melting temperature for Cu is calculated using a combination of free energy calculation and temperature sweep.
The EAM potential we will use is : Mishin, Y., M. J. Mehl, D. A. Papaconstantopoulos, A. F. Voter, and J. D. Kress. “Structural Stability and Lattice Defects in Copper: Ab Initio , Tight-Binding, and Embedded-Atom Calculations.” Physical Review B 63, no. 22 (May 21, 2001): 224106.
The calculation block gives the input conditions at which the calculation is carried out. First of all, the
ts, a temperature sweep over the temperature range given in the
temperature keyword will be done.
fcc lattice are chosen for both calculations, the former with reference phase
solid and the latter with the reference phase
liquid. The potential file is specified in
pair_coeff command in the
The calculation can be run by,
calphy -i input.yaml
Once submitted, it should give a message
Total number of 2 calculations found. It will also create a set of folders with the names
mode-lattice-referencephase-temperature-pressure. In this case, there will be
ts-fcc-liquid-1200-0. If there are any errors in the calculation, it will be recoreded in
ts-fcc-liquid-1200-0.sub.err. Once the calculation starts, a log file called
calphy.log will be created in the
aforementioned folders. The file gives some information about the preparation stage.
ts mode consists of two stages, the first step is the calculation of free energy at 1200 K, followed by a sweep until 1400 K. The results of the free energy calculation is recorded in
report.yaml file. The file is shown below:
average: density: 0.07557340045029688 vol_atom: 13.2345681628356 input: concentration: '1.0' element: Cu lattice: fcc pressure: 0.0 temperature: 1200 results: error: 0.0 free_energy: -4.058705625895811 pv: 0.0 reference_system: 1.9111881193348468 unit: eV/atom work: 4.742856409041955
In the file, the average and input quantities are recorded. The more interesting block is the
results block. Here the calculated free energy value in eV/atom is given in the
free_energy key. The free energy of the reference system is given in
reference_system and the work done in switching is under
error key gives the error in the calculation. In this case, its 0 as we ran only a single loop (see
The variation of the free energy within the temperature range is given in the
temperature_sweep.dat files instead of each of the folders. The file contains three columns, temperature, free energy and the error in free energy. The files are read in and plotted below.
import numpy as np import matplotlib.pyplot as plt
st, sfe, sferr = np.loadtxt("ts-fcc-solid-1200-0/temperature_sweep.dat", unpack=True) lt, lfe, lferr = np.loadtxt("ts-fcc-liquid-1200-0/temperature_sweep.dat", unpack=True)
plt.plot(st, sfe, color="#E53935", label="solid") plt.plot(lt, lfe, color="#0097A7", label="liquid") plt.xlabel("Temperature (K)", fontsize=12) plt.ylabel("F (ev/atom)", fontsize=12) plt.legend() plt.savefig("tm.png", dpi=300, bbox_inches="tight")
From the plot, at temperatures below 1300 K, solid is the more stable structure with lower free energy. Around 1340 K, liquid becomes the more stable structure. We can find the temperature at which the free energy of both phases are equal, which is the melting temperature.
args = np.argsort(np.abs(sfe-lfe)) print(st[args], "K")