From 6c8e25569d6c757affec5270a485709e761715a8 Mon Sep 17 00:00:00 2001 From: Tao Cheng Date: Wed, 4 Apr 2018 17:34:25 -0700 Subject: [PATCH] debug --- lib/mytype.py | 2 + tools/dftb/get_input.py | 11 ++-- tools/ion/getcheng.py | 2 +- tools/ion/tocheng.py | 2 +- tools/vasp/freq/vib-1.py | 137 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 148 insertions(+), 6 deletions(-) create mode 100755 tools/vasp/freq/vib-1.py diff --git a/lib/mytype.py b/lib/mytype.py index 759a39e..2e65d2d 100644 --- a/lib/mytype.py +++ b/lib/mytype.py @@ -209,6 +209,8 @@ def assignEleTypes(self,): counter = 0 for i in self.atoms: b = i.name.strip() + if "_" in b: + b = b.split("_")[0] match = pattern.match(b) if match: i.element= match.group(1) diff --git a/tools/dftb/get_input.py b/tools/dftb/get_input.py index ac01b5c..d1065f6 100644 --- a/tools/dftb/get_input.py +++ b/tools/dftb/get_input.py @@ -3,15 +3,16 @@ # http://www.dftb.org/parameters/download/3ob/3ob-3-1-cc/ HD = { "C": -0.1492, "H": -0.1857, "O":-0.1575, "N":-0.1535, "Cu":-0.2000, "S":-0.11, "Na":-0.0454, "F":-0.1623, - "Ca":-0.0340, "Cl": -0.0697, "Mg":-0.02, - } + "Ca":-0.0340, "Cl": -0.0697, "Mg":-0.02, "P": -0.14, + "K": -0.0339, } # MaxAngularMomentum MM = { "C": "p", "H": "s", "O":"p", "N":"p", "Cu":"d", "S":"d", - "Na": "p", "F":"p", "Ca":"p", "Cl": "p", "Mg":"p",} + "Na": "p", "F":"p", "Ca":"p", "Cl": "p", "Mg":"p", "P": "p", + "K": "p",} # CovalentRadius # http://periodictable.com/Properties/A/CovalentRadius.html CR = {"C":0.77, "H":0.37, "O": 0.73, "N":0.75, "Cu":1.38, "S":1.02, - "Na":0.71, "F":0.71, "Ca":1.76, "Cl":0.99,} + "Na":0.71, "F":0.71, "Ca":1.76, "Cl":0.99, "P":1.07, "K":2.03,} # HybridPolarisations HP = { "C": [1.382, 1.382, 1.382, 1.064, 1.064, 1.064, 3.8, 3.8, 3.8, 3.8, 3.8, 3.8, 2.50], @@ -22,7 +23,9 @@ "S": [1.030, 1.030, 1.090, 1.090, 1.090, 1.090, 3.8, 3.8, 3.8, 3.8, 3.8, 3.8, 2.82], "Na": [0.386, 0.386, 0.000, 0.000, 0.000, 0.000, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 0.80], # fake values +"K": [0.386, 0.386, 0.000, 0.000, 0.000, 0.000, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 0.80], "F": [1.030, 1.030, 1.090, 1.090, 1.090, 1.090, 3.8, 3.8, 3.8, 3.8, 3.8, 3.8, 2.82], +"P": [1.030, 1.030, 1.090, 1.090, 1.090, 1.090, 3.8, 3.8, 3.8, 3.8, 3.8, 3.8, 2.82], "Ca": [1.030, 1.030, 1.090, 1.090, 1.090, 1.090, 3.8, 3.8, 3.8, 3.8, 3.8, 3.8, 2.82], "Cl": [1.030, 1.030, 1.090, 1.090, 1.090, 1.090, 3.8, 3.8, 3.8, 3.8, 3.8, 3.8, 2.82], } diff --git a/tools/ion/getcheng.py b/tools/ion/getcheng.py index ca8b4bd..fe1cf71 100755 --- a/tools/ion/getcheng.py +++ b/tools/ion/getcheng.py @@ -9,7 +9,7 @@ if len(sys.argv) > 1: for i in range(1, len(sys.argv)): files += sys.argv[i] + " " - cmd += "tao@131.215.26.220:" + cmd += "tao@131.215.26.245:" cmd += files cmd += " ." print cmd diff --git a/tools/ion/tocheng.py b/tools/ion/tocheng.py index 4d9086f..23a2b1b 100755 --- a/tools/ion/tocheng.py +++ b/tools/ion/tocheng.py @@ -11,7 +11,7 @@ def main(): args = parser.parse_args() files = args.fname folder = args.o - ip = "131.215.26.220" + ip = "131.215.26.219" usr = "tao" cmd = "scp -r " cmd += " ".join(files) diff --git a/tools/vasp/freq/vib-1.py b/tools/vasp/freq/vib-1.py new file mode 100755 index 0000000..66e6632 --- /dev/null +++ b/tools/vasp/freq/vib-1.py @@ -0,0 +1,137 @@ +#!/usr/bin/env python + +############################################################################### +# Python script to calculate free energy contribution (ZPE included) +# of vibrational modes +# Written by Hai Xiao +############################################################################### + +import sys +from math import exp,log +#import os + +if len(sys.argv) == 4: + T = float(sys.argv[1]) + outcar = open(sys.argv[2]) + state = int(sys.argv[3]) + + # Get vibrational eigenvalues from OUTCAR + vibmatch = "f =" + ivibmatch = "f/i=" + Evib = [] + Nivib = 0 + Eivib = [] + for line in outcar: + if vibmatch in line: + # Read in vibrational mode in unit of meV (hv, NOT 1/2*hv) + Evib.append(float(line.split()[9])) + if ivibmatch in line: + Eivib.append(float(line.split()[8])) + Nivib = Nivib + 1 + outcar.close() + ###for i in Evib: + ### print "%10.6f"%(i) + + # Calculate characteristic vibrational temperatures + # conversion factor from http://physics.nist.gov/cgi-bin/cuu/Value?evk + # using the source: 2014 CODATA recommended values + # T = hv/k + meV2K = 11.6045221 + Tvib = [] + for i in Evib: + Tvib.append(i*meV2K) + ###for i in Tvib: + ### print "%12.6f"%(i) + + # Set the threshold for low frequency modes, which might lead to unphysically large entropy contribution + # Currently use the peak frequency 60 cm-1 of ice Ih spectra, corresponding to + # the acoustic translational mode of the six-member rings + # Planck constant in eV*s from http://physics.nist.gov/cgi-bin/cuu/Value?hev + # using the source: 2014 CODATA recommended values + h = float("4.135667662E-15") + # speed of light in vacuum in unit of cm/s (exact) + c = float("2.99792458E10") + # conversion factor for wavenumber 1/lambda (cm-1): h*c + cm2eV = h*c + ###print "%10.6E"%(cm2eV) + threshold = 60.0*cm2eV*1000.0*meV2K + ###print "%12.6f"%(threshold) + for i in range(len(Tvib)): + if Tvib[i] < threshold: + print "Warning: very small frequency found, which might lead to unphysically large entropy contribution" + print "very small frequency %3i: %12.6f K"%(i+1,Tvib[i]) + Tvib[i] = threshold + print " treated as real: %12.6f K"%(Tvib[i]) + ###for i in Tvib: + ### print "%12.6f"%(i) + + # Important check of small imaginary frequencies + if state == 0: + if Nivib > 0: + print "Warning: imaginary frequencies found for supposed local minimum !!!" + print "Warning: all imaginary frequencies treated as threshold ones for now" + for i in range(Nivib): + print "imaginary frequency %i: %12.6f K"%(i+1,Eivib[i]*meV2K) + Tvib.append(threshold) + print " treated as real: %12.6f K"%(threshold) + elif state == 1: + if Nivib == 0: + print "Error: no imaginary frequency found for supposed transition state !!!" + sys.exit(0) + elif Nivib > 1: + print "Warning: more than 1 imaginary frequency found for supposed transition state !!!" + print "Warning: all imaginary frequencies, except the lowest, treated as threshold ones for now" + for i in range(Nivib-1): + print "imaginary frequency %i: %12.6f K"%(i+1,Eivib[i]*meV2K) + Tvib.append(threshold) + print " treated as real: %12.6f K"%(threshold) + print "imaginary frequency %i: %12.6f K"%(Nivib,Eivib[Nivib-1]*meV2K) + print "treated as the only 1 imaginary frequency for supposed transition state" + else: + print "Error: only accept state 0 or 1 for now" + sys.exit(0) + + # Free energy contribution, using formulas from + # http://www.gaussian.com/g_whitepap/thermo.htm + # R = N_A * k + # use k here, in order to get final energies in the unit of eV + # Boltzmann constant from http://physics.nist.gov/cgi-bin/cuu/Value?tkev + # using the source: 2014 CODATA recommended values + k = float("8.6173303E-5") + # ZPE + H_zpe = 0.0 + for i in Tvib: + H_zpe = H_zpe + k * i * 0.5 + ###print "%10.6f"%(H_zpe) + ###H_zpe_test = 0.0 + ###for i in Evib: + ### H_zpe_test = H_zpe_test + 0.5 * i / 1000.0 + ###print "%10.6f"%(H_zpe_test) + # enthalpy and entropy contributions + H_vib = 0.0 + S_vib = 0.0 + if T == 0.0: + H_vib = 0.0 + S_vib = 0.0 + else: + for i in Tvib: + q = i/T + H_vib = H_vib + k * i / (exp(q) - 1.0) + S_vib = S_vib + k * (q / (exp(q) - 1.0) - log(1.0 - exp(-q))) + ###print "%13.6E"%(H_vib) + ###print "%13.6E"%(S_vib) + G_vib = H_vib - T * S_vib + total_vib = H_zpe + G_vib + print "ZPE contribution is %10.6f eV"%(H_zpe) + print "Enthalpy contribution H at %.1f K is %10.6f eV"%(T,H_vib) + print "Entropy S at %.1f K is %10.6f eV/K"%(T,S_vib) + print "H - TS at %.1f K is %10.6f eV"%(T,G_vib) + print "Total contribution at %.1f K with ZPE is %10.6f eV"%(T,total_vib) +else: + # Print usage + print "" + print "Usage:" + print "vib2G.py [T(K)] [OUTCAR] [0 or 1 for LM or TS]" + print "" + sys.exit(0) +