一個腳本兩步計算材料Raman譜（附數據處理和繪圖腳本）

????????在以往推送中已經介紹了相當多的計算材料Raman的方法，使用的軟件主要為Phonopy-Spectroscopy，相關軟件還有vasp，phonopy，phono3py等。

Phonopy-Spectroscopy計算材料紅外和Raman光譜

Phonopy-Spectroscopy 計算紅外和拉曼光譜

????也有一些教程介紹了無需額外軟件即可得到材料Raman信息的腳本vasp_raman.py，但腳本缺少自由度并且需要額外的數據處理。

????本篇推送即介紹這一腳本計算材料Raman信息的流程。

????腳本和案例可在作者處下載：

https://github.com/ZeliangSu/VaspRoutine/blob/main/raman-sc/

在本文最后也會附上0.6.0版本的腳本

準備工作

計算可得到無虛頻的聲子譜的結構文件

這里使用單晶硅作案例

 ? Si1.0 ? ? 5.4661639157319968 ? ?0.0000000000000000 ? ?0.0000000000000000 ? ? 0.0000000000000000 ? ?5.4661639157319968 ? ?0.0000000000000000 ? ? 0.0000000000000000 ? ?0.0000000000000000 ? ?5.4661639157319968 ? ? Si ? 8Direct ?0.8750000000000000 ?0.8750000000000000 ?0.8750000000000000 ?0.8750000000000000 ?0.3750000000000000 ?0.3750000000000000 ?0.3750000000000000 ?0.8750000000000000 ?0.3750000000000000 ?0.3750000000000000 ?0.3750000000000000 ?0.8750000000000000 ?0.1250000000000000 ?0.1250000000000000 ?0.1250000000000000 ?0.1250000000000000 ?0.6250000000000000 ?0.6250000000000000 ?0.6250000000000000 ?0.1250000000000000 ?0.6250000000000000 ?0.6250000000000000 ?0.6250000000000000 ?0.1250000000000000

這里已將聲子譜反折疊為原胞第一布里淵區內的聲子色散，并將縱坐標單位設置為（cm-1）

計算第一步：頻率信息

計算INCAR參考：

SYSTEM = Si_bulkISTART = 0 # From-scratch; job : 0-new 1-cont 2-samecut NWRITE = 3 Verbosity
! electronic relaxationENCUT = 300.0 # cut-off energyPREC = Accurate # precision : accurate/normal/low ISPIN = 1 # 1 - off, 2 - on (non spin-polarized calculation)ICHARG = 2 # > 10 for non-SC calculationIALGO = 38 # DAVidson, then RMM-DIISEDIFF = 1.0E-8 # defaultISMEAR = 0 # gaussianSIGMA = 0.05
! PAW'sLREAL = .FALSE. # default - Automatic choice of how projection is doneADDGRID = .TRUE.
! phononsIBRION = 5POTIM = 0.01
! parallelisationLPLANE = .FALSE.KPAR=8
! outputLWAVE = .FALSE. # WAVECAR fileLCHARG = .FALSE. # CHCAR fileLELF = .FALSE.LVTOT = .FALSE.

????特別注意參數：NWRITE =?3?Verbosity

計算完成后可自行處理聲子色散和相關數據。材料的頻率信息在OUTCAR中可自行查看。

計算第二步：宏觀介電張量

? ? 將第一步計算結果中的POSCAR復制為POSCAR.phon，將OUTCAR復制為OUTCAR.phon。

第二步INCAR參考：

SYSTEM = Si_bulkISTART = 0 # From-scratch; job : 0-new 1-cont 2-samecut NWRITE = 3 Verbosity
! electronic relaxationENCUT = 300.0 # cut-off energyPREC = Accurate # precision : accurate/normal/low ISPIN = 1 # 1 - off, 2 - on (non spin-polarized calculation)ICHARG = 2 # > 10 for non-SC calculationIALGO = 38 # DAVidson, then RMM-DIISEDIFF = 1.0E-8 # defaultISMEAR = 0 # gaussianSIGMA = 0.05
! PAW'sLREAL = .FALSE. # default - Automatic choice of how projection is doneADDGRID = .TRUE.
! phonons!! IBRION = 5!! POTIM = 0.01
LEPSILON=.TRUE.
! parallelisationLPLANE = .FALSE.KPAR=8
! outputLWAVE = .FALSE. # WAVECAR fileLCHARG = .FALSE. # CHCAR fileLELF = .FALSE.LVTOT = .FALSE.

此時需要調用腳本vasp_raman.py?進行第二步計算命令為：

export VASP_RAMAN_RUN='aprun -B /u/afonari/vasp.5.3.2/vasp.5.3/vasp &> job.out'export VASP_RAMAN_PARAMS='01_21_2_0.01'
python vasp_raman.py > vasp_raman.out

????第一行中?VASP_RAMAN_RUN為服務器等計算資源可使用的vasp的執行命令

注意：在集群使用隊列資源排隊進行計算的時候需要將上面三行命令都填寫進隊列系統腳本中，并替代原有的執行命令，具體請根據實際計算情況更改

?????第二行中VASP_RAMAN_PARAMS參數為Raman計算針對的頻率范圍和計算設置。

第一個數為起始頻率編號，01為起始頻率編號；

第二個數為截至頻率編號，最大值為POSCAR總原子數×3。

計算總任務數為第二個數減去第一個數再乘第三個數。

計算過程中腳本會讀取OUTCAR.phon中的頻率信息針對不同頻率的聲子所對應的原子振動模式對結構施加微擾。

同時會實時將新生成的OUTCAR另存最后統一數據處理。

計算結果保存在vasp_raman.dat

# mode ? ?freq(cm-1) ? ?alpha ? ?beta2 ? ?activity001 ? 517.72301 ?-0.0219904 ?520.9354925 ?3646.5702084002 ? 517.71996 ? 0.0177803 ?521.9308522 ?3653.5301918003 ? 517.71816 ?-0.0130798 ?521.3205054 ?3649.2512364004 ? 446.13270 ?-0.0022072 ? 0.0007120 ? 0.0052035005 ? 446.12649 ?-0.0026568 ? 0.0001072 ? 0.0010683006 ? 446.11229 ?-0.0091150 ? 0.0020607 ? 0.0181638007 ? 446.10094 ? 0.0056815 ? 0.0001899 ? 0.0027820008 ? 446.09466 ? 0.0108317 ? 0.0016816 ? 0.0170508009 ? 446.08791 ? 0.0082157 ? 0.0009541 ? 0.0097161010 ? 397.45270 ?-0.0004496 ? 0.0001839 ? 0.0012966011 ? 397.45083 ?-0.0050684 ? 0.0004061 ? 0.0039986012 ? 397.44702 ? 0.0011036 ? 0.0020508 ? 0.0144103013 ? 397.44649 ?-0.0029021 ? 0.0011642 ? 0.0085283014 ? 397.44295 ?-0.0007766 ? 0.0000371 ? 0.0002871015 ? 397.44138 ? 0.0047005 ? 0.0005162 ? 0.0046073016 ? 130.93534 ?-0.0010219 ? 0.0000026 ? 0.0000654017 ? 130.92851 ? 0.0026977 ? 0.0000319 ? 0.0005511018 ? 130.92701 ?-0.0001635 ? 0.0000030 ? 0.0000221019 ? 130.92623 ?-0.0006949 ? 0.0000040 ? 0.0000499020 ? 130.92383 ?-0.0003270 ? 0.0000030 ? 0.0000261021 ? 130.92351 ? 0.0000817 ? 0.0000035 ? 0.0000247

????可見Raman活性頻率為517cm-1，與Si實驗值520cm-1相當接近?(J.H. Parker,?et al., Phys Rev, 155, 712 (1967))。

????如果吹毛求疵，可以在計算頻率或聲子時便通過高精度結構優化和使用實驗值的晶格常數等方法將所得對應的頻率和實驗值矯正。

? ? 使用腳本將已經得到Raman信息處理（可自行選擇擬合方式Gaussian或Lorentzian）

 python broaden.py vasp_raman.dat

附：

vasp_raman.py腳本：

https://github.com/ZeliangSu/VaspRoutine/blob/main/raman-sc/vasp_raman.py

#!/usr/bin/env python## vasp_raman.py v. 0.6.0## Raman off-resonant activity calculator# using VASP as a back-end.## Contributors: Alexandr Fonari (Georgia Tech)# Shannon Stauffer (UT Austin)## URL: http://raman-sc.github.io## MIT license, 2013 - 2016#import reimport sysdef MAT_m_VEC(m, v): ? ?p = [ 0.0 for i in range(len(v)) ] ? ?for i in range(len(m)): ? ? ? ?assert len(v) == len(m[i]), 'Length of the matrix row is not equal to the length of the vector' ? ? ? ?p[i] = sum( [ m[i][j]*v[j] for j in range(len(v)) ] ) ? ?return pdef T(m): ? ?p = [[ m[i][j] for i in range(len( m[j] )) ] for j in range(len( m )) ] ? ?return pdef parse_poscar(poscar_fh): ? ?# modified subroutine from phonopy 1.8.3 (New BSD license) ? ?# ? ?poscar_fh.seek(0) # just in case ? ?lines = poscar_fh.readlines() ? ?# ? ?print(lines) ? ?scale = float(lines[1]) ? ?if scale < 0.0: ? ? ? ?print("[parse_poscar]: ERROR negative scale not implemented.") ? ? ? ?sys.exit(1) ? ?# ? ?b = [] ? ?for i in range(2, 5): ? ? ? ?b.append([float(x)*scale for x in lines[i].split()[:3]]) ? ?# ? ?vol = b[0][0]*b[1][1]*b[2][2] + b[1][0]*b[2][1]*b[0][2] + b[2][0]*b[0][1]*b[1][2] - \ ? ? ? ? ?b[0][2]*b[1][1]*b[2][0] - b[2][1]*b[1][2]*b[0][0] - b[2][2]*b[0][1]*b[1][0] ? ?# ? ?try: ? ? ? ?num_atoms = [int(x) for x in lines[5].split()] ? ? ? ?line_at = 6 ? ?except ValueError: ? ? ? ?symbols = [x for x in lines[5].split()] ? ? ? ?num_atoms = [int(x) for x in lines[6].split()] ? ? ? ?line_at = 7 ? ?nat = sum(num_atoms) ? ?# ? ?if lines[line_at][0].lower() == 's': ? ? ? ?line_at += 1 ? ?# ? ?if (lines[line_at][0].lower() == 'c' or lines[line_at][0].lower() == 'k'): ? ? ? ?is_scaled = False ? ?else: ? ? ? ?is_scaled = True ? ?# ? ?line_at += 1 ? ?# ? ?positions = [] ? ?for i in range(line_at, line_at + nat): ? ? ? ?pos = [float(x) for x in lines[i].split()[:3]] ? ? ? ?# ? ? ? ?if is_scaled: ? ? ? ? ? ?pos = MAT_m_VEC(T(b), pos) ? ? ? ?# ? ? ? ?positions.append(pos) ? ?# ? ?poscar_header = ''.join(lines[1:line_at-1]) # will add title and 'Cartesian' later ? ?return nat, vol, b, positions, poscar_headerdef parse_env_params(params): ? ?tmp = params.strip().split('_') ? ?if len(tmp) != 4: ? ? ? ?print("[parse_env_params]: ERROR there should be exactly four parameters") ? ? ? ?sys.exit(1) ? ?# ? ?[first, last, nderiv, step_size] = [int(tmp[0]), int(tmp[1]), int(tmp[2]), float(tmp[3])] ? ?# ? ?return first, last, nderiv, step_size#### subs for the output from VTST toolsdef parse_freqdat(freqdat_fh, nat): ? ?freqdat_fh.seek(0) # just in case ? ?# ? ?eigvals = [ 0.0 for i in range(nat*3) ] ? ?# ? ?for i in range(nat*3): # all frequencies should be supplied, regardless of requested to calculate ? ? ? ?tmp = freqdat_fh.readline().split() ? ? ? ?eigvals[i] = float(tmp[0]) ? ?# ? ?return eigvals#def parse_modesdat(modesdat_fh, nat): ? ?from math import sqrt ? ?modesdat_fh.seek(0) # just in case ? ?# ? ?eigvecs = [ 0.0 for i in range(nat*3) ] ? ?norms = ? [ 0.0 for i in range(nat*3) ] ? ?# ? ?for i in range(nat*3): # all frequencies should be supplied, regardless of requested to calculate ? ? ? ?eigvec = [] ? ? ? ?for j in range(nat): ? ? ? ? ? ?tmp = modesdat_fh.readline().split() ? ? ? ? ? ?eigvec.append([ float(tmp[x]) for x in range(3) ]) ? ? ? ?# ? ? ? ?modesdat_fh.readline().split() # empty line ? ? ? ?eigvecs[i] = eigvec ? ? ? ?norms[i] = sqrt( sum( [abs(x)**2 for sublist in eigvec for x in sublist] ) ) ? ?# ? ?return eigvecs, norms#### end subs for VTST#def get_modes_from_OUTCAR(outcar_fh, nat): ? ?from math import sqrt ? ?eigvals = [ 0.0 for i in range(nat*3) ] ? ?eigvecs = [ 0.0 for i in range(nat*3) ] ? ?norms ? = [ 0.0 for i in range(nat*3) ] ? ?# ? ?outcar_fh.seek(0) # just in case ? ?while True: ? ? ? ?line = outcar_fh.readline() ? ? ? ?if not line: ? ? ? ? ? ?break ? ? ? ?# ? ? ? ?if "Eigenvectors after division by SQRT(mass)" in line: ? ? ? ? ? ?outcar_fh.readline() # empty line ? ? ? ? ? ?outcar_fh.readline() # Eigenvectors and eigenvalues of the dynamical matrix ? ? ? ? ? ?outcar_fh.readline() # ---------------------------------------------------- ? ? ? ? ? ?outcar_fh.readline() # empty line ? ? ? ? ? ?# ? ? ? ? ? ?for i in range(nat*3): # all frequencies should be supplied, regardless of those requested to calculate ? ? ? ? ? ? ? ?outcar_fh.readline() # empty line ? ? ? ? ? ? ? ?p = re.search(r'^\s*(\d+).+?([\.\d]+) cm-1', outcar_fh.readline()) ? ? ? ? ? ? ? ?eigvals[i] = float(p.group(2)) ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ?outcar_fh.readline() # X ? ? ? ? Y ? ? ? ? Z ? ? ? ? ? dx ? ? ? ? ?dy ? ? ? ? ?dz ? ? ? ? ? ? ? ?eigvec = [] ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ?for j in range(nat): ? ? ? ? ? ? ? ? ? ?tmp = outcar_fh.readline().split() ? ? ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ? ? ?eigvec.append([ float(tmp[x]) for x in range(3,6) ]) ? ? ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ?eigvecs[i] = eigvec ? ? ? ? ? ? ? ?norms[i] = sqrt( sum( [abs(x)**2 for sublist in eigvec for x in sublist] ) ) ? ? ? ? ? ?# ? ? ? ? ? ?return eigvals, eigvecs, norms ? ? ? ?# ? ?print("[get_modes_from_OUTCAR]: ERROR Couldn't find 'Eigenvectors after division by SQRT(mass)' in OUTCAR. Use 'NWRITE=3' in INCAR. Exiting...") ? ?sys.exit(1)#def get_epsilon_from_OUTCAR(outcar_fh): ? ?epsilon = [] ? ?# ? ?outcar_fh.seek(0) # just in case ? ?while True: ? ? ? ?line = outcar_fh.readline() ? ? ? ?if not line: ? ? ? ? ? ?break ? ? ? ?# ? ? ? ?if "MACROSCOPIC STATIC DIELECTRIC TENSOR" in line: ? ? ? ? ? ?try: ? ? ? ? ? ? ? ?outcar_fh.readline() ? ? ? ? ? ? ? ?epsilon.append([float(x) for x in outcar_fh.readline().split()]) ? ? ? ? ? ? ? ?epsilon.append([float(x) for x in outcar_fh.readline().split()]) ? ? ? ? ? ? ? ?epsilon.append([float(x) for x in outcar_fh.readline().split()]) ? ? ? ? ? ?except: ? ? ? ? ? ? ? ?from lxml import etree as ET ? ? ? ? ? ? ? ?doc = ET.parse('vasprun.xml') ? ? ? ? ? ? ? ?epsilon = [[float(x) for x in c.text.split()] for c in doc.xpath("/modeling/calculation/varray")[3].getchildren()] ? ? ? ? ? ?return epsilon ? ?# ? ?raise RuntimeError("[get_epsilon_from_OUTCAR]: ERROR Couldn't find dielectric tensor in OUTCAR") ? ?return 1#if __name__ == '__main__': ? ?from math import pi ? ?from shutil import move ? ?import os ? ?import datetime ? ?import time ? ?#import argparse ? ?import optparse ? ?# ? ?print("") ? ?print(" ? ?Raman off-resonant activity calculator,") ? ?print(" ? ?using VASP as a back-end.") ? ?print("") ? ?print(" ? ?Contributors: Alexandr Fonari ?(Georgia Tech)") ? ?print(" ? ? ? ? ? ? ? ? ?Shannon Stauffer (UT Austin)") ? ?print(" ? ?MIT License, 2013") ? ?print(" ? ?URL: http://raman-sc.github.io") ? ?print(" ? ?Started at: "+datetime.datetime.now().strftime("%Y-%m-%d %H:%M")) ? ?print("") ? ?# ? ?description ?= "Before run, set environment variables:\n" ? ?description += " ? ?VASP_RAMAN_RUN='mpirun vasp'\n" ? ?description += " ? ?VASP_RAMAN_PARAMS='[first-mode]_[last-mode]_[nderiv]_[step-size]'\n\n" ? ?description += "bash one-liner is:\n" ? ?description += "VASP_RAMAN_RUN='mpirun vasp' VASP_RAMAN_PARAMS='1_2_2_0.01' python vasp_raman.py" ? ?# ? ?parser = optparse.OptionParser(description=description) ? ?parser.add_option('-g', '--gen', help='Generate POSCAR only', action='store_true') ? ?parser.add_option('-u', '--use_poscar', help='Use provided POSCAR in the folder, USE WITH CAUTION!!', action='store_true') ? ?(options, args) = parser.parse_args() ? ?#args = vars(parser.parse_args()) ? ?args = vars(options) ? ?# ? ?VASP_RAMAN_RUN = os.environ.get('VASP_RAMAN_RUN') ? ?if VASP_RAMAN_RUN == None: ? ? ? ?print("[__main__]: ERROR Set environment variable 'VASP_RAMAN_RUN'") ? ? ? ?print("") ? ? ? ?parser.print_help() ? ? ? ?sys.exit(1) ? ?print("[__main__]: VASP_RAMAN_RUN='"+VASP_RAMAN_RUN+"'") ? ?# ? ?VASP_RAMAN_PARAMS = os.environ.get('VASP_RAMAN_PARAMS') ? ?if VASP_RAMAN_PARAMS == None: ? ? ? ?print("[__main__]: ERROR Set environment variable 'VASP_RAMAN_PARAMS'") ? ? ? ?print("") ? ? ? ?parser.print_help() ? ? ? ?sys.exit(1) ? ?print("[__main__]: VASP_RAMAN_PARAMS='"+VASP_RAMAN_PARAMS+"'") ? ?# ? ?first, last, nderiv, step_size = parse_env_params(VASP_RAMAN_PARAMS) ? ?assert first >= 1, ? ?'[__main__]: First mode should be equal or larger than 1' ? ?assert last >= first, '[__main__]: Last mode should be equal or larger than first mode' ? ?if args['gen']: assert last == first, "[__main__]: '-gen' mode -> only generation for the one mode makes sense" ? ?assert nderiv == 2, ? '[__main__]: At this time, nderiv = 2 is the only supported' ? ?disps = [-1, 1] ? ? ?# hardcoded for ? ?coeffs = [-0.5, 0.5] # three point stencil (nderiv=2) ? ?# ? ?try: ? ? ? ?poscar_fh = open('POSCAR.phon', 'r') ? ?except IOError: ? ? ? ?print("[__main__]: ERROR Couldn't open input file POSCAR.phon, exiting...\n") ? ? ? ?sys.exit(1) ? ?# ? ?# nat, vol, b, poscar_header = parse_poscar_header(poscar_fh) ? ?nat, vol, b, pos, poscar_header = parse_poscar(poscar_fh) ? ?print(pos) ? ?#print poscar_header ? ?#sys.exit(0) ? ?# ? ?# either use modes from vtst tools or VASP ? ?if os.path.isfile('freq.dat') and os.path.isfile('modes_sqrt_amu.dat'): ? ? ? ?try: ? ? ? ? ? ?freqdat_fh = open('freq.dat', 'r') ? ? ? ?except IOError: ? ? ? ? ? ?print("[__main__]: ERROR Couldn't open freq.dat, exiting...\n") ? ? ? ? ? ?sys.exit(1) ? ? ? ?# ? ? ? ?eigvals = parse_freqdat(freqdat_fh, nat) ? ? ? ?freqdat_fh.close() ? ? ? ?# ? ? ? ?try:  ? ? ? ? ? ?modes_fh = open('modes_sqrt_amu.dat' , 'r') ? ? ? ?except IOError: ? ? ? ? ? ?print("[__main__]: ERROR Couldn't open modes_sqrt_amu.dat, exiting...\n") ? ? ? ? ? ?sys.exit(1) ? ? ? ?# ? ? ? ?eigvecs, norms = parse_modesdat(modes_fh, nat) ? ? ? ?modes_fh.close() ? ?# ? ?elif os.path.isfile('OUTCAR.phon'): ? ? ? ?try: ? ? ? ? ? ?outcar_fh = open('OUTCAR.phon', 'r') ? ? ? ?except IOError: ? ? ? ? ? ?print("[__main__]: ERROR Couldn't open OUTCAR.phon, exiting...\n") ? ? ? ? ? ?sys.exit(1) ? ? ? ?# ? ? ? ?eigvals, eigvecs, norms = get_modes_from_OUTCAR(outcar_fh, nat) ? ? ? ?outcar_fh.close() ? ?# ? ?else: ? ? ? ?print("[__main__]: Neither OUTCAR.phon nor freq.dat/modes_sqrt_amu.dat were found, nothing to do, exiting...") ? ? ? ?sys.exit(1) ? ?# ? ?output_fh = open('vasp_raman.dat', 'w') ? ?output_fh.write("# mode ? ?freq(cm-1) ? ?alpha ? ?beta2 ? ?activity\n") ? ?for i in range(first-1, last): ? ? ? ?eigval = eigvals[i] ? ? ? ?eigvec = eigvecs[i] ? ? ? ?norm = norms[i] ? ? ? ?# ? ? ? ?print("") ? ? ? ?print("[__main__]: Mode #%i: frequency %10.7f cm-1; norm: %10.7f" % ( i+1, eigval, norm )) ? ? ? ?# ? ? ? ?ra = [[0.0 for x in range(3)] for y in range(3)] ? ? ? ?for j in range(len(disps)): ? ? ? ? ? ?disp_filename = 'OUTCAR.%04d.%+d.out' % (i+1, disps[j]) ? ? ? ? ? ?# ? ? ? ? ? ?try: ? ? ? ? ? ? ? ?outcar_fh = open(disp_filename, 'r') ? ? ? ? ? ? ? ?print("[__main__]: File "+disp_filename+" exists, parsing...") ? ? ? ? ? ?except IOError: ? ? ? ? ? ? ? ?if args['use_poscar'] != True: ? ? ? ? ? ? ? ? ? ?print("[__main__]: File "+disp_filename+" not found, preparing displaced POSCAR") ? ? ? ? ? ? ? ? ? ?poscar_fh = open('POSCAR', 'w') ? ? ? ? ? ? ? ? ? ?poscar_fh.write("%s %4.1e \n" % (disp_filename, step_size)) ? ? ? ? ? ? ? ? ? ?poscar_fh.write(poscar_header) ? ? ? ? ? ? ? ? ? ?poscar_fh.write("Cartesian\n") ? ? ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ? ? ?for k in range(nat): ? ? ? ? ? ? ? ? ? ? ? ?pos_disp = [ pos[k][l] + eigvec[k][l]*step_size*disps[j]/norm for l in range(3)] ? ? ? ? ? ? ? ? ? ? ? ?poscar_fh.write( '%15.10f %15.10f %15.10f\n' % (pos_disp[0], pos_disp[1], pos_disp[2]) ) ? ? ? ? ? ? ? ? ? ? ? ?#print '%10.6f %10.6f %10.6f %10.6f %10.6f %10.6f' % (pos[k][0], pos[k][1], pos[k][2], dis[k][0], dis[k][1], dis[k][2]) ? ? ? ? ? ? ? ? ? ?poscar_fh.close() ? ? ? ? ? ? ? ?else: ? ? ? ? ? ? ? ? ? ?print("[__main__]: Using provided POSCAR") ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ?if args['gen']: # only generate POSCARs ? ? ? ? ? ? ? ? ? ?poscar_fn = 'POSCAR.%+d.out' % disps[j] ? ? ? ? ? ? ? ? ? ?move('POSCAR', poscar_fn) ? ? ? ? ? ? ? ? ? ?print("[__main__]: '-gen' mode -> "+poscar_fn+" with displaced atoms have been generated") ? ? ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ? ? ?if j+1 == len(disps): # last iteration for the current displacements list ? ? ? ? ? ? ? ? ? ? ? ?print("[__main__]: '-gen' mode -> POSCAR files with displaced atoms have been generated, exiting now") ? ? ? ? ? ? ? ? ? ? ? ?sys.exit(0) ? ? ? ? ? ? ? ?else: # run VASP here ? ? ? ? ? ? ? ? ? ?print("[__main__]: Running VASP...") ? ? ? ? ? ? ? ? ? ?os.system(VASP_RAMAN_RUN) ? ? ? ? ? ? ? ? ? ?try: ? ? ? ? ? ? ? ? ? ? ? ?move('OUTCAR', disp_filename) ? ? ? ? ? ? ? ? ? ?except IOError: ? ? ? ? ? ? ? ? ? ? ? ?print("[__main__]: ERROR Couldn't find OUTCAR file, exiting...") ? ? ? ? ? ? ? ? ? ? ? ?sys.exit(1) ? ? ? ? ? ? ? ? ? ?# ? ? ? ? ? ? ? ? ? ?outcar_fh = open(disp_filename, 'r') ? ? ? ? ? ?# ? ? ? ? ? ?try: ? ? ? ? ? ? ? ?eps = get_epsilon_from_OUTCAR(outcar_fh) ? ? ? ? ? ? ? ?outcar_fh.close() ? ? ? ? ? ?except Exception as err: ? ? ? ? ? ? ? ?print(err) ? ? ? ? ? ? ? ?print("[__main__]: Moving "+disp_filename+" back to 'OUTCAR' and exiting...") ? ? ? ? ? ? ? ?move(disp_filename, 'OUTCAR') ? ? ? ? ? ? ? ?sys.exit(1) ? ? ? ? ? ?# ? ? ? ? ? ?for m in range(3): ? ? ? ? ? ? ? ?for n in range(3): ? ? ? ? ? ? ? ? ? ?ra[m][n] ? += eps[m][n] * coeffs[j]/step_size * norm * vol/(4.0*pi) ? ? ? ? ? ?#units: A^2/amu^1/2 = ? ? ? ? dimless ? * 1/A ? ? ? ? * 1/amu^1/2 ?* A^3 ? ? ? ?# ? ? ? ?alpha = (ra[0][0] + ra[1][1] + ra[2][2])/3.0 ? ? ? ?beta2 = ( (ra[0][0] - ra[1][1])**2 + (ra[0][0] - ra[2][2])**2 + (ra[1][1] - ra[2][2])**2 + 6.0 * (ra[0][1]**2 + ra[0][2]**2 + ra[1][2]**2) )/2.0 ? ? ? ?print("") ? ? ? ?print("! %4i ?freq: %10.5f ?alpha: %10.7f ?beta2: %10.7f ?activity: %10.7f " % (i+1, eigval, alpha, beta2, 45.0*alpha**2 + 7.0*beta2)) ? ? ? ?output_fh.write("%03i ?%10.5f ?%10.7f ?%10.7f ?%10.7f\n" % (i+1, eigval, alpha, beta2, 45.0*alpha**2 + 7.0*beta2)) ? ? ? ?output_fh.flush() ? ?# ? ?output_fh.close()

繪圖腳本：

https://github.com/ZeliangSu/VaspRoutine/blob/main/raman-sc/broaden.py（有修改）

#!/usr/bin/env python#def to_plot(hw,ab,gam=0.001,type='lorentzian'): ? ?import numpy as np ? ?# ? ?fmin = min(hw) ? ?fmax = max(hw) ? ?erange = np.arange(fmin-40*gam,fmax+40*gam,gam/10)#np.arange(fmin-40*gam,fmax+40*gam,gam/10) ? ?spectrum = 0.0*erange ? ?for i in range(len(hw)): ? ? ? ?if type=='Gaussian': ? ? ? ? ? ?spectrum += (2*np.pi)**(-.5)/gam*np.exp(np.clip(-1.0*(hw[i]-erange)**2/(2*gam**2),-300,300)) ? ? ? ?elif type=='Lorentzian': ? ? ? ? ? ?spectrum += ab[i]*1/np.pi*gam/((hw[i]-erange)**2+gam**2) ? ?# ? ?return erange, spectrum#if __name__ == '__main__': ? ?import numpy as np ? ?import sys ? ?# ? ?hw=np.genfromtxt(sys.argv[1], dtype=float) ? ?cm1 = hw[:,1] ? ?int1 = hw[:,4] ? ?int1 /= np.max(np.abs(int1),axis=0) ? ?Es1,Spectrum1 = to_plot(cm1, int1, 15.0, 'Lorentzian')# ? ?filename = "new-broaden.dat" ? ?print( "Writing", filename) ? ?f = open(filename,'w') ? ?f.write('# freq/cm-1 ?Intensity \n')
 ? ?for i in range(len(Es1)): ? ? ? ?f.write('%.5e ? %.5e\n' % (Es1[i],Spectrum1[i])) ? ?f.close() ? ?import matplotlib.pyplot as plt ? ? ?import numpy as np ? ?
 ? ?data = np.loadtxt("new-broaden.dat", unpack=True) ? ?
plt.plot(data[0], data[1]) ? ? plt.show()