Raman+VASP计算的Python脚本运行出现问题

muxijiao

最近在尝试用vasp计算Raman光谱，在https://github.com/raman-sc/VASP#installation 找到个python程序，按照说明运行时出现个问题 ，不知道是什么原因？求Python大神们解答，不过运行他自带的例子就没有问题。
python代码如下

1. #!/usr/bin/env python
   
2. #
   
3. # vasp_raman.py v. 0.6.0
   
4. #
   
5. # Raman off-resonant activity calculator
   
6. # using VASP as a back-end.
   
7. #
   
8. # Contributors: Alexandr Fonari (Georgia Tech)
   
9. # Shannon Stauffer (UT Austin)
   
10. #
    
11. # URL: http://raman-sc.github.io
    
12. #
    
13. # MIT license, 2013 - 2016
    
14. #
    
15. 
    
16. 
    
17. import re
    
18. import sys
    
19. 
    
20. 
    
21. def MAT_m_VEC(m, v):
    
22.     p = [ 0.0 for i in range(len(v)) ]
    
23.     for i in range(len(m)):
    
24.         assert len(v) == len(m[i]), 'Length of the matrix row is not equal to the length of the vector'
    
25.         p[i] = sum( [ m[i][j]*v[j] for j in range(len(v)) ] )
    
26.     return p
    
27. 
    
28. 
    
29. def T(m):
    
30.     p = [[ m[i][j] for i in range(len( m[j] )) ] for j in range(len( m )) ]
    
31.     return p
    
32. 
    
33. 
    
34. def parse_poscar(poscar_fh):
    
35.     # modified subroutine from phonopy 1.8.3 (New BSD license)
    
36.     #
    
37.     poscar_fh.seek(0) # just in case
    
38.     lines = poscar_fh.readlines()
    
39.     #
    
40.     scale = float(lines[1])
    
41.     if scale < 0.0:
    
42.         print "[parse_poscar]: ERROR negative scale not implemented."
    
43.         sys.exit(1)
    
44.     #
    
45.     b = []
    
46.     for i in range(2, 5):
    
47.         b.append([float(x)*scale for x in lines[i].split()[:3]])
    
48.     #
    
49.     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] - \
    
50.           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]
    
51.     #
    
52.     try:
    
53.         num_atoms = [int(x) for x in lines[5].split()]
    
54.         line_at = 6
    
55.     except ValueError:
    
56.         symbols = [x for x in lines[5].split()]
    
57.         num_atoms = [int(x) for x in lines[6].split()]
    
58.         line_at = 7
    
59.     nat = sum(num_atoms)
    
60.     #
    
61.     if lines[line_at][0].lower() == 's':
    
62.         line_at += 1
    
63.     #
    
64.     if (lines[line_at][0].lower() == 'c' or lines[line_at][0].lower() == 'k'):
    
65.         is_scaled = False
    
66.     else:
    
67.         is_scaled = True
    
68.     #
    
69.     line_at += 1
    
70.     #
    
71.     positions = []
    
72.     for i in range(line_at, line_at + nat):
    
73.         pos = [float(x) for x in lines[i].split()[:3]]
    
74.         #
    
75.         if is_scaled:
    
76.             pos = MAT_m_VEC(T(b), pos)
    
77.         #
    
78.         positions.append(pos)
    
79.     #
    
80.     poscar_header = ''.join(lines[1:line_at-1]) # will add title and 'Cartesian' later
    
81.     return nat, vol, b, positions, poscar_header
    
82. 
    
83. 
    
84. def parse_env_params(params):
    
85.     tmp = params.strip().split('_')
    
86.     if len(tmp) != 4:
    
87.         print "[parse_env_params]: ERROR there should be exactly four parameters"
    
88.         sys.exit(1)
    
89.     #
    
90.     [first, last, nderiv, step_size] = [int(tmp[0]), int(tmp[1]), int(tmp[2]), float(tmp[3])]
    
91.     #
    
92.     return first, last, nderiv, step_size
    
93. 
    
94. 
    
95. #### subs for the output from VTST tools
    
96. def parse_freqdat(freqdat_fh, nat):
    
97.     freqdat_fh.seek(0) # just in case
    
98.     #
    
99.     eigvals = [ 0.0 for i in range(nat*3) ]
    
100.     #
     
101.     for i in range(nat*3): # all frequencies should be supplied, regardless of requested to calculate
     
102.         tmp = freqdat_fh.readline().split()
     
103.         eigvals[i] = float(tmp[0])
     
104.     #
     
105.     return eigvals
     
106. #
     
107. def parse_modesdat(modesdat_fh, nat):
     
108.     from math import sqrt
     
109.     modesdat_fh.seek(0) # just in case
     
110.     #
     
111.     eigvecs = [ 0.0 for i in range(nat*3) ]
     
112.     norms =   [ 0.0 for i in range(nat*3) ]
     
113.     #
     
114.     for i in range(nat*3): # all frequencies should be supplied, regardless of requested to calculate
     
115.         eigvec = []
     
116.         for j in range(nat):
     
117.             tmp = modesdat_fh.readline().split()
     
118.             eigvec.append([ float(tmp[x]) for x in range(3) ])
     
119.         #
     
120.         modesdat_fh.readline().split() # empty line
     
121.         eigvecs[i] = eigvec
     
122.         norms[i] = sqrt( sum( [abs(x)**2 for sublist in eigvec for x in sublist] ) )
     
123.     #
     
124.     return eigvecs, norms
     
125. #### end subs for VTST
     
126. #
     
127. def get_modes_from_OUTCAR(outcar_fh, nat):
     
128.     from math import sqrt
     
129.     eigvals = [ 0.0 for i in range(nat*3) ]
     
130.     eigvecs = [ 0.0 for i in range(nat*3) ]
     
131.     norms   = [ 0.0 for i in range(nat*3) ]
     
132.     #
     
133.     outcar_fh.seek(0) # just in case
     
134.     while True:
     
135.         line = outcar_fh.readline()
     
136.         if not line:
     
137.             break
     
138.         #
     
139.         if "Eigenvectors after division by SQRT(mass)" in line:
     
140.             outcar_fh.readline() # empty line
     
141.             outcar_fh.readline() # Eigenvectors and eigenvalues of the dynamical matrix
     
142.             outcar_fh.readline() # ----------------------------------------------------
     
143.             outcar_fh.readline() # empty line
     
144.             #
     
145.             for i in range(nat*3): # all frequencies should be supplied, regardless of those requested to calculate
     
146.                 outcar_fh.readline() # empty line
     
147.                 p = re.search(r'^\s*(\d+).+?([\.\d]+) cm-1', outcar_fh.readline())
     
148.                 eigvals[i] = float(p.group(2))
     
149.                 #
     
150.                 outcar_fh.readline() # X         Y         Z           dx          dy          dz
     
151.                 eigvec = []
     
152.                 #
     
153.                 for j in range(nat):
     
154.                     tmp = outcar_fh.readline().split()
     
155.                     #
     
156.                     eigvec.append([ float(tmp[x]) for x in range(3,6) ])
     
157.                     #
     
158.                 eigvecs[i] = eigvec
     
159.                 norms[i] = sqrt( sum( [abs(x)**2 for sublist in eigvec for x in sublist] ) )
     
160.             #
     
161.             return eigvals, eigvecs, norms
     
162.         #
     
163.     print "[get_modes_from_OUTCAR]: ERROR Couldn't find 'Eigenvectors after division by SQRT(mass)' in OUTCAR. Use 'NWRITE=3' in INCAR. Exiting..."
     
164.     sys.exit(1)
     
165. #
     
166. def get_epsilon_from_OUTCAR(outcar_fh):
     
167.     epsilon = []
     
168.     #
     
169.     outcar_fh.seek(0) # just in case
     
170.     while True:
     
171.         line = outcar_fh.readline()
     
172.         if not line:
     
173.             break
     
174.         #
     
175.         if "MACROSCOPIC STATIC DIELECTRIC TENSOR" in line:
     
176.             outcar_fh.readline()
     
177.             epsilon.append([float(x) for x in outcar_fh.readline().split()])
     
178.             epsilon.append([float(x) for x in outcar_fh.readline().split()])
     
179.             epsilon.append([float(x) for x in outcar_fh.readline().split()])
     
180.             return epsilon
     
181.     #
     
182.     raise RuntimeError("[get_epsilon_from_OUTCAR]: ERROR Couldn't find dielectric tensor in OUTCAR")
     
183.     return 1
     
184. #
     
185. if __name__ == '__main__':
     
186.     from math import pi
     
187.     from shutil import move
     
188.     import os
     
189.     import datetime
     
190.     import time
     
191.     #import argparse
     
192.     import optparse
     
193.     #
     
194.     print ""
     
195.     print "    Raman off-resonant activity calculator,"
     
196.     print "    using VASP as a back-end."
     
197.     print ""
     
198.     print "    Contributors: Alexandr Fonari  (Georgia Tech)"
     
199.     print "                  Shannon Stauffer (UT Austin)"
     
200.     print "    MIT License, 2013"
     
201.     print "    URL: http://raman-sc.github.io"
     
202.     print "    Started at: "+datetime.datetime.now().strftime("%Y-%m-%d %H:%M")
     
203.     print ""
     
204.     #
     
205.     description  = "Before run, set environment variables:\n"
     
206.     description += "    VASP_RAMAN_RUN='mpirun vasp'\n"
     
207.     description += "    VASP_RAMAN_PARAMS='[first-mode]_[last-mode]_[nderiv]_[step-size]'\n\n"
     
208.     description += "bash one-liner is:\n"
     
209.     description += "VASP_RAMAN_RUN='mpirun vasp' VASP_RAMAN_PARAMS='1_2_2_0.01' python vasp_raman.py"
     
210.     #
     
211.     parser = optparse.OptionParser(description=description)
     
212.     parser.add_option('-g', '--gen', help='Generate POSCAR only', action='store_true')
     
213.     parser.add_option('-u', '--use_poscar', help='Use provided POSCAR in the folder, USE WITH CAUTION!!', action='store_true')
     
214.     (options, args) = parser.parse_args()
     
215.     #args = vars(parser.parse_args())
     
216.     args = vars(options)
     
217.     #
     
218.     VASP_RAMAN_RUN = os.environ.get('VASP_RAMAN_RUN')
     
219.     if VASP_RAMAN_RUN == None:
     
220.         print "[__main__]: ERROR Set environment variable 'VASP_RAMAN_RUN'"
     
221.         print ""
     
222.         parser.print_help()
     
223.         sys.exit(1)
     
224.     print "[__main__]: VASP_RAMAN_RUN='"+VASP_RAMAN_RUN+"'"
     
225.     #
     
226.     VASP_RAMAN_PARAMS = os.environ.get('VASP_RAMAN_PARAMS')
     
227.     if VASP_RAMAN_PARAMS == None:
     
228.         print "[__main__]: ERROR Set environment variable 'VASP_RAMAN_PARAMS'"
     
229.         print ""
     
230.         parser.print_help()
     
231.         sys.exit(1)
     
232.     print "[__main__]: VASP_RAMAN_PARAMS='"+VASP_RAMAN_PARAMS+"'"
     
233.     #
     
234.     first, last, nderiv, step_size = parse_env_params(VASP_RAMAN_PARAMS)
     
235.     assert first >= 1,    '[__main__]: First mode should be equal or larger than 1'
     
236.     assert last >= first, '[__main__]: Last mode should be equal or larger than first mode'
     
237.     if args['gen']: assert last == first, "[__main__]: '-gen' mode -> only generation for the one mode makes sense"
     
238.     assert nderiv == 2,   '[__main__]: At this time, nderiv = 2 is the only supported'
     
239.     disps = [-1, 1]      # hardcoded for
     
240.     coeffs = [-0.5, 0.5] # three point stencil (nderiv=2)
     
241.     #
     
242.     try:
     
243.         poscar_fh = open('POSCAR.phon', 'r')
     
244.     except IOError:
     
245.         print "[__main__]: ERROR Couldn't open input file POSCAR.phon, exiting...\n"
     
246.         sys.exit(1)
     
247.     #
     
248.     # nat, vol, b, poscar_header = parse_poscar_header(poscar_fh)
     
249.     nat, vol, b, pos, poscar_header = parse_poscar(poscar_fh)
     
250.     print pos
     
251.     #print poscar_header
     
252.     #sys.exit(0)
     
253.     #
     
254.     # either use modes from vtst tools or VASP
     
255.     if os.path.isfile('freq.dat') and os.path.isfile('modes_sqrt_amu.dat'):
     
256.         try:
     
257.             freqdat_fh = open('freq.dat', 'r')
     
258.         except IOError:
     
259.             print "[__main__]: ERROR Couldn't open freq.dat, exiting...\n"
     
260.             sys.exit(1)
     
261.         #
     
262.         eigvals = parse_freqdat(freqdat_fh, nat)
     
263.         freqdat_fh.close()
     
264.         #
     
265.         try:
     
266.             modes_fh = open('modes_sqrt_amu.dat' , 'r')
     
267.         except IOError:
     
268.             print "[__main__]: ERROR Couldn't open modes_sqrt_amu.dat, exiting...\n"
     
269.             sys.exit(1)
     
270.         #
     
271.         eigvecs, norms = parse_modesdat(modes_fh, nat)
     
272.         modes_fh.close()
     
273.     #
     
274.     elif os.path.isfile('OUTCAR.phon'):
     
275.         try:
     
276.             outcar_fh = open('OUTCAR.phon', 'r')
     
277.         except IOError:
     
278.             print "[__main__]: ERROR Couldn't open OUTCAR.phon, exiting...\n"
     
279.             sys.exit(1)
     
280.         #
     
281.         eigvals, eigvecs, norms = get_modes_from_OUTCAR(outcar_fh, nat)
     
282.         outcar_fh.close()
     
283.     #
     
284.     else:
     
285.         print "[__main__]: Neither OUTCAR.phon nor freq.dat/modes_sqrt_amu.dat were found, nothing to do, exiting..."
     
286.         sys.exit(1)
     
287.     #
     
288.     output_fh = open('vasp_raman.dat', 'w')
     
289.     output_fh.write("# mode    freq(cm-1)    alpha    beta2    activity\n")
     
290.     for i in range(first-1, last):
     
291.         eigval = eigvals[i]
     
292.         eigvec = eigvecs[i]
     
293.         norm = norms[i]
     
294.         #
     
295.         print ""
     
296.         print "[__main__]: Mode #%i: frequency %10.7f cm-1; norm: %10.7f" % ( i+1, eigval, norm )
     
297.         #
     
298.         ra = [[0.0 for x in range(3)] for y in range(3)]
     
299.         for j in range(len(disps)):
     
300.             disp_filename = 'OUTCAR.%04d.%+d.out' % (i+1, disps[j])
     
301.             #
     
302.             try:
     
303.                 outcar_fh = open(disp_filename, 'r')
     
304.                 print "[__main__]: File "+disp_filename+" exists, parsing..."
     
305.             except IOError:
     
306.                 if args['use_poscar'] != True:
     
307.                     print "[__main__]: File "+disp_filename+" not found, preparing displaced POSCAR"
     
308.                     poscar_fh = open('POSCAR', 'w')
     
309.                     poscar_fh.write("%s %4.1e \n" % (disp_filename, step_size))
     
310.                     poscar_fh.write(poscar_header)
     
311.                     poscar_fh.write("Cartesian\n")
     
312.                     #
     
313.                     for k in range(nat):
     
314.                         pos_disp = [ pos[k][l] + eigvec[k][l]*step_size*disps[j]/norm for l in range(3)]
     
315.                         poscar_fh.write( '%15.10f %15.10f %15.10f\n' % (pos_disp[0], pos_disp[1], pos_disp[2]) )
     
316.                         #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])
     
317.                     poscar_fh.close()
     
318.                 else:
     
319.                     print "[__main__]: Using provided POSCAR"
     
320.                 #
     
321.                 if args['gen']: # only generate POSCARs
     
322.                     poscar_fn = 'POSCAR.%+d.out' % disps[j]
     
323.                     move('POSCAR', poscar_fn)
     
324.                     print "[__main__]: '-gen' mode -> "+poscar_fn+" with displaced atoms have been generated"
     
325.                     #
     
326.                     if j+1 == len(disps): # last iteration for the current displacements list
     
327.                         print "[__main__]: '-gen' mode -> POSCAR files with displaced atoms have been generated, exiting now"
     
328.                         sys.exit(0)
     
329.                 else: # run VASP here
     
330.                     print "[__main__]: Running VASP..."
     
331.                     os.system(VASP_RAMAN_RUN)
     
332.                     try:
     
333.                         move('OUTCAR', disp_filename)
     
334.                     except IOError:
     
335.                         print "[__main__]: ERROR Couldn't find OUTCAR file, exiting..."
     
336.                         sys.exit(1)
     
337.                     #
     
338.                     outcar_fh = open(disp_filename, 'r')
     
339.             #
     
340.             try:
     
341.                 eps = get_epsilon_from_OUTCAR(outcar_fh)
     
342.                 outcar_fh.close()
     
343.             except Exception, err:
     
344.                 print err
     
345.                 print "[__main__]: Moving "+disp_filename+" back to 'OUTCAR' and exiting..."
     
346.                 move(disp_filename, 'OUTCAR')
     
347.                 sys.exit(1)
     
348.             #
     
349.             for m in range(3):
     
350.                 for n in range(3):
     
351.                     ra[m][n]   += eps[m][n] * coeffs[j]/step_size * norm * vol/(4.0*pi)
     
352.             #units: A^2/amu^1/2 =         dimless   * 1/A         * 1/amu^1/2  * A^3
     
353.         #
     
354.         alpha = (ra[0][0] + ra[1][1] + ra[2][2])/3.0
     
355.         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
     
356.         print ""
     
357.         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)
     
358.         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))
     
359.         output_fh.flush()
     
360.     #
     
361.     output_fh.close()
     

复制代码

sky

还不如去github上issues里给作者反馈错误

muxijiao

sky 发表于 2017-12-6 23:20
还不如去github上issues里给作者反馈错误

写邮件了

小梅的龙猫

我用了楼主所下载的那个程序，运行没有问题，就是计算有点慢。建议查看自己前面几步计算是否有错误，我最初用的时候也会报错，后面修改INCAR里面的一些参数，仔细核对以后是没有问题的。

exwhite

你好，我最近在看vasp计算Raman,试着算了 https://github.com/raman-sc/VASP里面  Sibulk-VASP 例子，根据自己的情况，修改了 里面 Raman.sub 的 export VASP_RAMAN_RUN='aprun -B /u/afonari/vasp.5.3.2/vasp.5.3/vasp &> job.out' 和 python /u/afonari/vasp_raman.py > vasp_raman.out， qsub Raman.sub，根据 README.md，qsub Raman.sub 之后计算完后，应该生成很多 OUTCAR.000* 和  vasp_raman.dat, 但是我并不能得到这些文件。刚开始接触这方面，可能犯了很傻的错误，可否指教我一下？非常感谢！

muxijiao

exwhite 发表于 2018-1-30 16:56
你好，我最近在看vasp计算Raman,试着算了 https://github.com/raman-sc/VASP里面  Sibulk-VASP 例子，根据 ...

这个。。我是自己的服务器做的，，可能是aprun那里调用问题，看看路径啊，命令啥的有没问题。不过你看看有没有一个OUTCAR生成，或者OSZICAr生成，是不是报错了啥的

exwhite

muxijiao 发表于 2018-2-8 16:23
这个。。我是自己的服务器做的，，可能是aprun那里调用问题，看看路径啊，命令啥的有没问题。不过你看看 ...

嗯，好的。谢谢~~~

FANYI

楼主好，请问下，VASP可否计算拉曼张量呀？从而计算拉曼强度？还是直接用vasp计算拉曼就行了呀？

mlzhao

你好，请问有vasp_raman的Python库吗

maoxinxina

exwhite 发表于 2018-1-30 16:56
**** 作者被禁止或删除 内容自动屏蔽 ****

你好，我在论坛上看见您在计算raman谱的时候没有出险OUTCAR.001等一系列文件，后来是怎么解决的呢。