Simulating the phase change of Argon

Using the HOOMD molucule dynamics software, I tried to reproduce the phase change of Argon. The model is a Lennard-Jones system with the constants defining the interaction from here: http://complex.elte.hu/~csabai/szamszim/simLec5.pdf .

First I tried to create a box of frozen Argon and warm it up, because trying to freeze a substance we always experience some overfreezing which makes it hard to detect the freezing point. In the other direction, this problem should not be that bad. But this way I could not reproduce an exact freezing point, the phase change was rather continuus. So in the end I'm showing system which were frozen with different freezing velocities.

I tried to analyze the phase change with constant pressure, so during the simulation I constantly downsize the box, to keep up with the changing temperature. But even with constant temperature the pressure has strong fluctuations, and it is not really possible to keep it very stable.

The exact temperature of phase change is not revealed in this simulation because cooling is too fast. I have tried to exactly simulate the boiling point too, but the strong fluctuation in the pressure, and the very-very slow convergence at a given temperature close to the boiling point made it unfeasible to exactly measure the real temperature of the phase transition.

Installing hoomd

conda config --add channels glotzer
conda install hoomd

Installing MDAnalysis

Load modules

In [1]:
import MDAnalysis
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
%matplotlib inline

Modify the quickstart scipt

In [2]:
def write_hoomd_script(n_steps):

    with open('nsteps'+str(n_steps)+'.hoomd','w') as f_h:
        f_h.write("""
#physical params
N=1000 #number of parts
L0=65 #box sizes
L_end=35

n_timesteps="""+str(n_steps)+""" #overall number of timesteps
dt=0.01 #integration dt
change_period=1e2 #changing properties at every period timestep
write_period=1e2 #writing at every period timestep


eps=1.0 # 1 is energy unit -> then T=E/kb= E * 119K
sigma=1.0
r_cut=3.0

T_0=90.0/119.8 #max temp in energy units
T_end=0.0/119.8 #min temp

#filenames
prefix='dump_dcd'

############################################
############################################
############################################

from hoomd_script import *
context.initialize()

# create random particles of name A
init.create_random(N=N, box=data.boxdim(L=L0), name='A')

# specify Lennard-Jones interactions between particle pairs
lj = pair.lj(r_cut=r_cut)
lj.pair_coeff.set('A', 'A', epsilon=eps, sigma=sigma)

# integrate mode
all = group.all();
integrate.mode_standard(dt=dt)
#recommended for rescale_temp
integrate.nve(group=all)
#integrate.nvt(group=all, T=1.2, tau=0.5)


# dump an xmle file
dump.xml(filename=prefix+'.xml',vis=True,velocity=True)

# dump a .dcd file for the trajectories
dump.dcd(filename=prefix+'.dcd', period=write_period)

#Cooling: goes from T_max at time step 0 to T_min at the end
update.rescale_temp(period=change_period,T=variant.linear_interp([(0, T_0), (n_timesteps, T_end)]))
update.box_resize(period=change_period, L = variant.linear_interp([(0,L0), (n_timesteps, L_end)]))


# analyze
analyze.log(filename='results"""+str(n_steps)+""".log', quantities=['pressure', 'temperature', 'kinetic_energy', 'potential_energy'],
            period=write_period, header_prefix='#')

# run 
run(n_timesteps)

""")

Create hoomd scripts for different cooling velocities

In [3]:
write_hoomd_script(int(1e4))
write_hoomd_script(int(1e5))
write_hoomd_script(int(1e6))

Run the simulations

In [4]:
%%bash
rm dump_dcd.*  hoomd.log
time hoomd nsteps10000.hoomd >> hoomd.log
rm dump_dcd.*
time hoomd nsteps100000.hoomd >> hoomd.log
rm dump_dcd.*
time hoomd nsteps1000000.hoomd >> hoomd.log

real	0m2.461s
user	0m2.384s
sys	0m0.308s

real	0m20.522s
user	0m20.048s
sys	0m0.696s

real	8m19.078s
user	8m12.864s
sys	0m5.784s

Load data

In [5]:
result4=pd.read_csv('results10000.log',sep='\t').iloc[1:,:] #skip first data
result5=pd.read_csv('results100000.log',sep='\t').iloc[1:,:] #skip first data
result6=pd.read_csv('results1000000.log',sep='\t').iloc[1:,:] #skip first data

Calculate some SI units

In [6]:
e_unit=1.65e-21
beta=119.8
d_unit=3.41e-10
p_unit=e_unit/d_unit**3
m_unit=1.660538921e-27

result4['temperature']*=beta
result4['pressure']*=p_unit
speed4=90/(1e4*0.01*np.sqrt(m_unit*d_unit**2/e_unit))

result5['temperature']*=beta
result5['pressure']*=p_unit
speed5=90/(1e5*0.01*np.sqrt(m_unit*d_unit**2/e_unit))

result6['temperature']*=beta
result6['pressure']*=p_unit
speed6=90/(1e6*0.01*np.sqrt(m_unit*d_unit**2/e_unit))

Plot the potential energy depending on the temperature

In the liquid state, the potential energy should be much higher than in the random gas state. With more points the cooling process becomes slower and the system has more time to start to form the liquid state. Unfortunately even the 1 million point simulation is actually very-very fast, and runs quite long. Running a simulation which condesates close to the theoretical boiling point looks hopeless.

In [7]:
fig,ax=plt.subplots()
fig.set_size_inches(12,9)

ax.plot(result4['temperature'],result4['potential_energy'],
        '.',mec='none',label='%e' % speed4 + r' $[K/s]$ (10K points)')
ax.plot(result5['temperature'],result5['potential_energy'],
        '.',mec='none',label='%e' % speed5 + r' $[K/s]$ (100K points)')
ax.plot(result6['temperature'],result6['potential_energy'],
        '.',mec='none',label='%e' % speed6 + r' $[K/s]$ (1M points)')

ax.axvline(83.81,c='m',lw=4,linestyle='dashed',label='theoretical melting point')
ax.axvline(87.302,c='orange',lw=4,linestyle='dashed',label='theoretical boiling point')

ax.legend(loc='best',fancybox=True,fontsize=16)
ax.set_ylabel(r'Potential energy [ $\epsilon$ ]',fontsize=16)
dump=ax.set_xlabel(r'$T [K]$',fontsize=16)

Check the pressure

  • It looks relatively constant, and close to atmospheric pressure
  • Note that int the condesated phase pressure starts to fluctuate very-very strongly.
In [8]:
fig,ax=plt.subplots()
fig.set_size_inches(12,9)
ax.plot(result6['temperature'],result6['pressure'],'.',mec='none',label='')
ax.axhline(1e5,lw=4,label='atmospheric pressure',c='g')
ax.legend(loc='best',fancybox=True,fontsize=16)
ax.set_xlabel(r'$T [K]$',fontsize=16)
ax.set_ylabel(r'$p [Pa]$',fontsize=16)
ax.set_yscale('log')
dump=ax.set_ylim(1e3,1e6)

Plot the particles

In [9]:
u = MDAnalysis.Universe('dump_dcd.xml','dump_dcd.dcd')
atoms=u.select_atoms('type A')

fig = plt.figure(figsize=(16,9))

u.trajectory[0] 
ax = fig.add_subplot(121, projection='3d')
plt.axis('off')
ax.set_title('gas phase')
dump=ax.plot(atoms.positions[:,0],atoms.positions[:,1],atoms.positions[:,2],'.')

u.trajectory[-1] 
ax = fig.add_subplot(122, projection='3d')
plt.axis('off')
ax.set_title('condensed phase')
dump=ax.plot(atoms.positions[:,0],atoms.positions[:,1],atoms.positions[:,2],'.')

Conclusion

I have managed to reconstruct the phase change of Argon, and the results look credible. Measuring the exact boiling point looks unfeasible because of the high pressure fluctuations, and the slow convergence around the boiling point. With much more particles, the pressure fluctuation could be remedied, and with much stronger hardware (i was running it on 1 core of an intel i7 ), the long simulations can also be carried out.