根据OUTCAR画出原子磁矩

准备文件OURCAR和POSCAR，使用下面脚本：

from tkinter import W
import numpy as np

def cart2frac(cart, latt):
    	# mt=np.transpose(latt)
	mt = np.linalg.inv(latt)
	mt = np.transpose(mt)
	frac = np.matmul(mt, cart)
	return frac


def frac2cart(frac, latt):
	mt = np.transpose(latt)
	cart = np.matmul(mt, frac)
	return cart

def readOUTCAR(nAtoms):
    
    MAGMOMlist=np.zeros((nAtoms, 3), dtype=np.float)	
    outcarfile = 'OUTCAR'
    fIn = open(outcarfile, 'r')
    lines = fIn.readlines()


    b=0
    for line in lines:
        b=b+1
        if ("magnetization (x)" in line):
            print(line+str(b)+'\n')                      #找出所有关键字所在的行
            aa=b                                          #把最后一次找到的值赋给aa
    print(aa)

    cc=0
    for i in range(aa+3,aa+3+nAtoms):
        temp=lines[i].rstrip().split()
        MAGMOMlist[cc,0]=temp[4]

        cc=cc+1


    print(MAGMOMlist,cc)

    b=0
    for line in lines:
        b=b+1
        if ("magnetization (y)" in line):
            print(line+str(b)+'\n')                      #找出所有关键字所在的行
            aa=b                                          #把最后一次找到的值赋给aa
    print(aa)

    cc=0
    for i in range(aa+3,aa+3+nAtoms):
        temp=lines[i].rstrip().split()
        MAGMOMlist[cc,1]=temp[4]

        cc=cc+1


    print(MAGMOMlist,cc)

    b=0
    for line in lines:
        b=b+1
        if ("magnetization (z)" in line):
            print(line+str(b)+'\n')                      #找出所有关键字所在的行
            aa=b                                          #把最后一次找到的值赋给aa
    print(aa)

    cc=0
    for i in range(aa+3,aa+3+nAtoms):
        temp=lines[i].rstrip().split()
        MAGMOMlist[cc,2]=temp[4]

        cc=cc+1


    print(MAGMOMlist,cc)
    return MAGMOMlist


def readPOSCAR(filename='POSCAR'):
	print('Open file:', filename, 'to read')
	f = open(filename, 'r')
	# first line 
	cellname=f.readline().strip()
	# second line 
	scale = float(f.readline().strip().split()[0])
	## line 3-5
	lat=np.zeros((3,3),dtype=np.float) 
	for i in range(3):
		s = f.readline().strip().split()
		for j in range(3):
			lat[i, j] = float(s[j]) * scale
	## line 6, namelist
	namelist = f.readline().strip().split()
	nTypes = len(namelist)
	## line 7, numlist
	numlist = [0] * nTypes
	s = f.readline().strip().split()
	for i in range(nTypes):
		numlist[i] = int(s[i])
	nAtoms = sum(numlist)
	## line 8, C/K or D
	s = f.readline().strip()
	if s[0] == 'C' or s[0] == 'c' or s[0] == 'K' or s[0] == 'k':
		cart=True
	else:
		cart=False
	## line 9-, coordinates
	fraclist = np.zeros((nAtoms, 3), dtype=np.float)
	cartlist = np.zeros((nAtoms, 3), dtype=np.float)		
	for i in range(nAtoms):
		s = f.readline().strip().split()
		for j in range(3):
			fraclist[i, j] = float(s[j])
		if cart is True:
			cartlist[i, :] = fraclist[i, :] * scale
			fraclist[i, :] = cart2frac(cartlist[i, :], lat)
		else:
			cartlist[i, :] = frac2cart(fraclist[i, :], lat)	

	return lat, nAtoms, namelist, numlist, cartlist, fraclist

def writePOSCAR(MAGMOMlist, lat, nAtoms, namelist, numlist, cartlist, fraclist, cart=True):
    f = open("./"+"POSCAR.xsf", 'w')
    f.write('CRYSTAL\n')
    f.write('PRIMVEC\n')
    # f.write(nAtoms,'\n')
    for i in range(3):
        f.write('%23.16f\t%23.16f\t%23.16f\n' % (lat[i, 0],lat[i, 1],lat[i, 2]))
    f.write('PRIMCOORD'+'\n')
    f.write(str(nAtoms)+'\n')
    a=0
    for i in range(len(namelist)):
        for j in range(numlist[i]):
            # f.write(str(namelist[i])+' '+str(cartlist[a])+'\n')
            f.write('%s\t%f\t%f\t%f\t%f\t%f\t%f\t\n'%(namelist[i],cartlist[a,0],cartlist[a,1],cartlist[a,2],MAGMOMlist[a,0],MAGMOMlist[a,1],MAGMOMlist[a,2]))
            a=a+1
	# for i in namelist:
	# 	f.write('%4s ' % i)
	# f.write('\n')
	# for i in numlist:
	# 	f.write('%4d ' % i)
	# f.write('\n')
	# if cart is True:
	# 	f.write('Cart\n')
	# 	for i in range(nAtoms):
	# 		f.write('%16.8f\t%16.8f\t%16.8f\n' % (cartlist[i,0],cartlist[i,1],cartlist[i,2]))
	# else:
	# 	f.write('Direct\n')
	# 	for i in range(nAtoms):
	# 		f.write('%16.8f\t%16.8f\t%16.8f\n' % (fraclist[i,0],fraclist[i,1],fraclist[i,2]))
	# f.close()
    
lat, nAtoms, namelist, numlist, cartlist, fraclist  = readPOSCAR()
MAGMOMlist=readOUTCAR(nAtoms)
writePOSCAR(MAGMOMlist, lat, nAtoms, namelist, numlist, cartlist, fraclist, cart=True)

print(lat, nAtoms, namelist, numlist )


#1、写入CRYSTAL PRIMVEC，这两个都是不变的
#2、读取POSCAR，把POSCAR分成3部分，分别为晶格常数，原子的名称和个数，原子的位置（笛卡尔坐标和分数坐标）
#3、把上面读取的晶格常数写进xsf文件的第三行
#4、写入PRIMCOORD
#5、写入所有原子的个数，上面POSCAR已经读取了
#6、读取outcar,找到最后一个magnetization出现的行数，读取MAGMOM
#7、把读取的整理出来放在xsf中

运行完之后生成POSCAR.sxf文件，拖入VESTA中显示，效果如下：