Entropy.py

go back to Main Page, Computational Resources, Scripts, Scripts for VASP or VASP

  File:Get mag.tgz 

This script will do the following:

Reads OUTCAR and calculate the entropy contributions of vibrations at desired temperature;

Note 1:

To get the vibration information, we need set:

IBRION = 5 
POTIM  = 0.02 (or smaller 0.01) 
NSW = 1 
NFREE = 2 

Note 2:

For gas and some surface adsorbates, we will get some imaginary frequencies. In this script, you can choose to conclude these imaginary frequencies or not.

To use this script:

1) make sure that OUTCAR is in the path when you run this script.

2) Use the command below:

   python entropy.py Temperature 

3) Choose to consider imaginary frequencies or not.

Enjoy!

---

# -*- coding: utf-8 -*-
#!/usr/bin/env python 
# Version-4.2: Changed the Reminder from Gas into Imaginary frequencies.
# -*- coding: utf-8 -*-
# Writen By Qiang Li from ICIQ in 23-09-2017
# Ref: Atkin's Physical chemistry, 10th Edition, Statistical thermodynamics: the vibrational Contribution, Page: 642
# To use it: python entropy-v4.py 273.15  (273.15 is the Temperature!!!)
# Then choose the calculation type by follow the suggestion: Gas calculation or others 

import sys
import math 
from scipy import constants as con

script, Tem = sys.argv
h_p = con.h   # 6.62606957E-34 # J*s Plank Constant
k_b = con.k   # 1.38064852E-23 # m²*kg*s⁻²*K⁻¹ Boltzman Constant 
R_gas = con.R # 8.3144598      # J*mol⁻¹*K⁻¹ Gas Constant 
l_s = con.c   # 299792458      # light speed m * s ⁻¹

Tem = float(Tem)  # Temperature 
beta = 1/(k_b * Tem)

def get_pf(miu): # get partition function
    x_i = h_p * float(miu) * l_s * beta
    pf_l  = x_i / (math.exp(x_i) -1)   # Lleft part of the partition fuction
    pf_r  = math.log(1-math.exp(-x_i)) 
    pf    = pf_l - pf_r 
    return pf 

look_up1 = 'cm-1'
look_up2 = 'f/i'
miu_list = []

status = raw_input('\nDo you want to consider the imaginary frequencies? 1: Yes. 2 No.: Please type 1 or 2\n>>>>\t')

def miu(status):
    with open('OUTCAR') as data_in:
        if status in ['1', 1, 'Yes', 'Y', 'y', 'yes']:  
            for i in data_in.readlines():
                if look_up1 in i.rstrip() and look_up2 not in i.rstrip():
                    miu_list.append(float(i.split()[7]))
                elif look_up2  in i.rstrip():
                    miu_list.append(float(i.split()[6]))
        elif status in  ['2', 2, 'No', 'N', 'n', 'no']: 
            for i in data_in.readlines():
                if look_up1 in i.rstrip() and look_up2 not in i.rstrip():
                    miu_list.append(float(i.split()[7]))                
        else: 
            print "\nPlease Run it again and Type: 1 or 2 !!!\n"
miu(status)
miu_list= list(set(miu_list))

sum_pf =  R_gas * sum(get_pf(i*100) for i in miu_list) # Calculate the Entropy S 

TS= Tem * sum_pf  # entropy contribution: T * S 

print 'S in J * K⁻¹ * mol⁻¹:\t%s' %(sum_pf)
print 'S in eV * K⁻¹ * mol⁻¹:\t%s' %(sum_pf/1000/96.4869)
print 'T*S in eV * mol⁻¹:\t%s' %(TS/1000/96.4869)