Thermopack

A package for the atomistic computation of free energy surfaces

Daniel Schwen and Alfredo Caro

The Gibbs Free Energy is a thermodynamic potential that contains every possible thermodynamic data for a given system at a given phase. Knowledge of the Free energy for example allows the construction of phase diagagrams (which show the most stable of many competing phases for given conditions).

Free energy is not immediately expressable as ensemble average. Mixed Hamiltonian with switching parameter $\lambda$ \[ \lambda \color{red}U + (1-\lambda) \color{green}W \] Free energy $F_{\color{green}W}$ of the harmonic oscilator system can be computed analytically \[ F_W=-kT\log\Omega_{\color{green}W}\quad,\quad\Omega_{\color{green}W}=\left(\frac T T_\text{Einstein}\right)^3 \] Evaluate the ensemble with respect to the pure Hamiltonians $\color{red}U$ (full interaction, CD-EAM) and $\color{green}W$ (harmonic oscillators) and compute \[ F_{\color{red}U}=F_{\color{green}W}+\underbrace{\int_0^1 \langle \color{red}U-\color{green}W\rangle d\lambda}_{\text{switching work}} \]

Documentation

Input files

The entire problem to be calculated is defined in a LAMMPS configuraton file fragment. Here all necessary physical data is specified, such as

Debye Temperature
expected melting point (approximate)
selected temperature for the Hamiltonian switching
Pressure at which to run
Simulation time (must be increased for mixtures of heavy and light elements)

This configuration file contains LAMMPS instructions as well as functional comment lines, which are interpreted by the runall.sh and analysis.sh scripts.

Example 1: Single run on a pure Ni crystal

# lattice lattice fcc 3.52 origin 0.25 0.25 0.25 region box block -8 8 -8 8 -8 8 create_box 2 box create_atoms 1 box # potential pair_style eam/alloy pair_coeff * * ../pot/Mishin-Ni-Al-2009.eam.alloy Ni Al neighbor 0.3 bin neigh_modify delay 3 # set Debye temperatures [in K] variable TDA equal 375 variable TDB equal 394 # set switching temperature [K] variable Tswitch equal 385 # set timestep [ps] # rule of thumb: sqrt(lightest mass/50)/1000 timestep 0.001 # set simulation step number multiplier (default is 1) # rule of thumb: sqrt(largest mass/lightest mass) variable Smul equal 1 # set approximate melting temperature [K] # (Overestimate rather than underestimate!) variable Tmelt equal 2000 # set upper end of the enthalpy temperature ramp [K] # (try to stay below the boiling point here) variable Thigh equal 3000 # set simulation pressure (kBar) # usually this should be 0.0 variable Spress equal 0.0 # the following magic comment tells thermopack not to iterate over concentrations, # but to run only the sample as is. #SINGLE_RUN

Example 2: Random solution FeCr alloy - G as a function of concentration

# lattice lattice bcc 2.87 origin 0.25 0.25 0.25 region box block -8 8 -8 8 -8 8 create_box 2 box create_atoms 1 box # potential pair_style eam/cd pair_coeff * * ../pot/FeCr.cdeam Fe Cr neighbor 0.3 bin neigh_modify delay 3 # set Debye temperatures [in K] variable TDA equal 426 variable TDB equal 395 # set switching temperature $T_0$ [K] variable Tswitch equal 410 # set timestep [ps] # rule of thumb: sqrt(lightest mass/50)/1000 timestep 0.001 # set approximate melting temperature [K] # (Overestimate rather than underestimate!) variable Tmelt equal 2100 # set upper end of the enthalpy temperature ramp [K] # (try to stay below the boiling point here) variable Thigh equal 4000 # set simulation step number multiplier (default is 1) # rule of thumb: sqrt(largest mass/lightest mass) variable Smul equal 1 # set simulation pressure (kBar) # usually this should be 0.0 variable Spress equal 0.0 # type of structure (replace Ni atoms) # the following incluse file will apply the concentration values to # generate a random solution sample based on the chrystal structure above include ../lib/in.random #CONCENTRATIONS (all numbers on the next line after the # are taken as concentrations to run the system at) # 1 0.985423 0.956268 0.927114 0.897959 0.868805 0.83965 0.810496 0.781341 0.752187 0.723032 0.693878 0.664723 0.635569 0.606414 0.577259 0.548105 0.51895 0.489796 0.460641 0.431487 0.402332 0.373178 0.344023 0.314869 0.285714 0.25656 0.227405 0.198251 0.169096 0.139942 0.110787 0.0816327 0.0524781 0.0451895 0.0422741 0.0379009 0.0349854 0.0306122 0.0276968 0.0233236 0.016035 0.00874636 0.0

Analysis

Thermopack includes scripts to process the data generated by the LAMMPS runs. The analyze scripts include a series of fitting procedures which are outlined below.

Performing an analysis

The data generated during the run phase is analyzed with the analyze.py Python script.

analyze.py takes the name of the run as an argument. It expects a configuration (in.config.name) and a run directory (run.name) in the current path. The generated data is saved to the run directory

Fitting the switching curves

The MD data $\langle U-W\rangle$ is fit using \[ f(\lambda)=a+b\lambda+c\lambda^2+\sqrt{d(1-\lambda)}+e\lambda^3+f\lambda^4 \] and integrated as \[ %\int_0^1 f(\lambda) d\lambda = \frac14e\lambda^4 + \frac15f\lambda^5 + \frac13c\lambda^3 + \frac12b\lambda^2 + ax - \frac23d(1.0-\lambda)^\frac32 \int_0^1 f(\lambda) d\lambda = a\lambda + \frac12b\lambda^2 + \frac13c\lambda^3 - \frac23d(1.0-\lambda)^\frac32 + \frac14e\lambda^4 + \frac15f\lambda^5 \] These integrals plotted as a function of concentration...

Fitting $G_0$

The free energy $G_0$ at the switching temperature $T_0$ is fit using a Redlich-Kister polynomial \[ %G_0(c) = \underbrace{cG_0(1)+(1-c)G_0(0)}_{\text{ideal solution}} + \underbrace{c(1-c)\left[\sum_{i=0}^5 f_i(1-2c)^i \right]+f_6(1-c)+f_7c}_{\text{Redlich-Kister polynomial}} G_0(c) = \underbrace{cG_0(1)+(1-c)G_0(0)}_{\text{ideal solution}} + \underbrace{c(1-c)\left[\sum_{i=0}^5 f_i(1-2c)^i \right]}_{\text{Redlich-Kister polynomial}} \] To obtain the free energy at all other temperatures we resort to the Gibbs-Duhem equation \[ \frac GT=\frac{G_0}{T_0} -\int_{T_0}^T \frac{H(\tau)}{\tau^2}d\tau \] Using MD we obtain enthalpy curves $H(T,c)$ in the desired $T$ range for all concentrations...

Fitting $H$

The enthalpy $H$ is fit using \[ H(T)=H_a+H_bT+H_cT^2+H_dT^3 \] Note that $H_{a,b,c,d}$ are functions of concentration...

Concentration dependent free energy surface

Free energy at switching temperature
Enthalpy at all temperatures and concentrations
Gibbs-Duhem equation \[ \frac GT=\frac{G_0}{T_0} -\int_{T_0}^T \frac{H(\tau)}{\tau^2}d\tau \]
Resulting fit conforms to CALPHAD form

$$ \begin{eqnarray} H_a(c) &=& c(1-c)\left[\sum_{i=0}^5 a_i(1-2c)^i \right]+a_6(1-c)+a_7c \\ H_{\kappa\in\{b,c,d\}}(c) &=& \sum_{i=0}^3 \kappa_i\cdot c^i \\ H_e(c) &=&-\frac{H_a(c)}{T_0} + H_c(c)T_0 + H_b(c)\log T_0 + H_d(c)\frac{T_0^2}2 \\ G_0(c) &=& c(1-c)\left[\sum_{i=0}^5 f_i(1-2c)^i \right]+f_6(1-c)+f_7c \\ \nonumber \color{red}{G(c,T)}&\color{red}=& \color{red}{G_0(c)\frac{T}{T_0} + H_a(c)-H_b(c)T\log T-H_c(c)T^2+H_d(c)\frac{T^3}2 + H_e(c)T + Tk_B\left(c\log c +(1-c)\log(1-c)\right)} \end{eqnarray} $$

Common tangent construction

The phase diagram is obtained through the common tangent construction by solving \[ \left.\frac{dG}{dc}\right|_{c_1} = \left.\frac{dG}{dc}\right|_{c_2} = \frac{G(c_1)-G(c_2)}{c_1-c_2} \] for every Temperature. (You might need to construct multiple free energy surfaces for different phases!)

The spinodal is found by solving \[ \left.\frac{d^2G}{dc^2}\right|_{c} = 0 \] for every temperature.

Modifications applied to LAMMPS

Modification to the stock LAMMPS, needed by and supplied along with Thermopack are listed below

pair hybrid/weighted, for mixing of multiple potentials (liquid)
fix mulforce, to set the contribution of EAM in the mixed Hamiltonian (solid)
restore_state command, to restore a configuration saved with fix store/state
variable tmass(n) function to obtain mass of atom type n
convenience function to calculate Debye frequencies and oscillator spring constants
fix vaf to compute velocity auto-correlation functions
and phonon spectra for verification of T_Debye
Allow for a spring constant $k=0$ in fix spring/self

Download

thermopack-24Jan15.tar.gz

Unpack the archive into the lammps src directory and install the package with make yes-THERMO. Please use an up-to-date LAMMPS version as the thermopack patches may not apply otherwise. Let me know if the thermopackage requires any updates. Contact info is contained in the README file in the archive.

Running the examples

After downloading and installing thermopack compile a new LAMMPS executable. To run the example problems please install the MANYBODY and USER-MISC packages.

make yes-thermo
make yes-manybody
make yes-user-misc
make mpi

Use the appropriate build configuration (mpi in this example). Finally set the LAMMPS environment variable to the full path of the newly built executable

export LAMMPS=~/Downloads/lammps-19Feb15.thermo/src/lmp_mpi

then in the THERMO/examples directory run the two examples by entering

cd THERMO/examples
. ./runexamples.sh

This will launch two runs that will take several hours to complete.

Citation

If you use thermopack in a publication, please cite our paper

Schwen, D., E. Martinez, and A. Caro. “On the Analytic Calculation of Critical Size for Alpha Prime Precipitation in FeCr.” Journal of Nuclear Materials 439, no. 1–3 (August 2013): 180–84. doi:10.1016/j.jnucmat.2013.03.057.