1. 随机位移的方法
给整个结构中的原子引入随机位移，代替md过程，增加采样
核心脚本
random_distance.py
需要的文件
sub-grade.sh POSCAR inv2
import numpy as np
from sys import argv
from ase.atoms import Atoms
from ase.io import read
import os
import subprocess
import shutil
def write2my(file_path, ene_i, lat_i, ele_i, coo_i, foc_i, vir_i=None):
lat_i = lat_i.reshape(3, 3)
coo_i = coo_i.reshape(-1, 3)
foc_i = foc_i.reshape(-1, 3)
with open(file_path, 'w') as file:
file.write(f"# total energy = {ene_i} eV\n\n")
if vir_i is not None:
vir_i = np.array(vir_i).reshape(-1)
file.write(f"VIRIAL\n")
for i in vir_i:
file.write(f'{i:20.8f}')
file.write("\n")
file.write("CRYSTAL\n")
file.write("PRIMVEC\n")
for j in lat_i:
for k in j:
file.write(f'{k:20.8f}')
file.write('\n')
file.write("PRIMCOORD\n")
file.write(f"{len(coo_i)} 1\n")
for j in range(len(coo_i)):
file.write(f'{ele_i[j]:2}')
# coo
for k in coo_i[j]:
file.write(f"{k:20.8f}")
# force
for k in foc_i[j]:
file.write(f"{k:20.8f}")
file.write("\n")
pass
def read_from_ase_atoms(atoms: Atoms):
try:
ene = atoms.get_potential_energy()
except:
ene = 0.0
lat = atoms.get_cell()
pos = atoms.get_positions()
try:
foc = atoms.get_forces()
except:
foc = np.zeros_like(pos)
try:
sts = atoms.get_stress()
xx, yy, zz, yz, xz, xy = - sts * atoms.get_volume()
vir = np.array(
[[xx, xy, xz],
[xy, yy, yz],
[xz, yz, zz]]).reshape(-1)
except:
vir = np.zeros([3, 3])
ele = atoms.get_chemical_symbols()
return ene, lat, pos, foc, vir, ele
def read_poscar(filename):
with open(filename,'r') as file:
lines =file.readlines()
title =lines[0].strip()
scale =float(lines[1].strip())
#for i in range(2,5):
# lattice_list.append(list(map(float,line[i].split())))
#lattice_array=np.array(lattice_list)
lattice_vectors = np.array([list(map(float,lines[i].split())) for i in range(2,5)])
element_type=lines[5].split()
atomnumbers=np.array(list(map(int,lines[6].split())))
atompositions=np.array([list(map(float,lines[i].split()[:3])) for i in range(8,8+sum(atomnumbers))])
cor_type = lines[7].strip()
if cor_type.lower() == "direct" or cor_type.lower() == "d":
atompositions[atompositions < 0]+=1
atompositions[atompositions >=1]-=1
elif cor_type.lower() == "cartesian" or cor_type.lower() == "c":
pass
return(title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions)
def write_poscar(filename,title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions):
with open(filename,'w') as file:
file.write(f"{title}\n")
file.write(f"{scale}\n")
for vec in lattice_vectors:
file.write(f"{vec[0]:.8f} {vec[1]:.8f} {vec[2]:.8f}\n")
#" ".join用于将可迭代对象的元素以空格连接起来
file.write(" ".join(element_type)+ "\n")
file.write(" ".join(map(str,atomnumbers))+ "\n")
file.write(f"{cor_type}\n")
for pos in atompositions:
file.write(f"{pos[0]:.8f} {pos[1]:.8f} {pos[2]:.8f}\n")
def change_dir_to_car(filename):
title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions= read_poscar(filename)
if cor_type.lower() == "cartesian" or cor_type.lower() == "c":
pass
elif cor_type.lower() == "direct" or cor_type.lower() == "d":
atompositions = np.dot(atompositions,lattice_vectors)
cor_type = "Cartesian"
return title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions
def atom_pos(atomnumbers,atompositions):
atom_pos = np.split(atompositions,np.cumsum(atomnumbers)[:-1])
return atom_pos
parent_dir = os.path.dirname(os.getcwd())
oos_dir = os.path.join(parent_dir, "oos")
os.makedirs(oos_dir, exist_ok=True)
counter = 0
counter2 = 0
counter3 = 0
exit_flag = False
for itera in range(1, 100):
if exit_flag:
break
mignitude1 = 0.0018*itera
mignitude2 = 0.0009*itera
mignitude3 = 0.0017*itera
mignitude4 = 0.0022*itera
title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions=change_dir_to_car("POSCAR")
atom_pos = atom_pos(atomnumbers,atompositions)
for i in range(1,10001):
counter3 +=1
if exit_flag:
break
counter +=1
rd1=np.random.randn(atomnumbers[0],3)*mignitude1
rd2=np.random.randn(atomnumbers[1],3)*mignitude1
rd3=np.random.randn(atomnumbers[2],3)*mignitude1
rd4=np.random.randn(atomnumbers[3],3)*mignitude1
atom_pos[0]+=rd1
atom_pos[1]+=rd2
atom_pos[2]+=rd3
atom_pos[3]+=rd4
atompositions = np.concatenate(atom_pos)
write_poscar(f"{counter}.vasp",title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions)
datas = read(f"{counter}.vasp", index=":")
ene, lat, pos, foc, vir, ele = read_from_ase_atoms(datas[0])
write2my(
f"{counter}.xsf",
ene_i=ene,
lat_i=lat,
ele_i=ele,
coo_i=pos,
foc_i=foc,
vir_i=vir)
os.remove(f"{counter}.vasp")
print(counter3)
if counter3 >= 100:
subprocess.run(["sh", "sub-grade.sh"], check=True)
with open("grade-log", 'r') as file:
lines = file.readlines()
for line in lines:
if exit_flag:
break
if 'grad:' in line:
parts = line.split()
filename = parts[0]
grade_value = float(parts[-1])
if 3 < grade_value < 8:
counter2 += 1
if counter2 >= 50:
exit_flag = True
break
src = os.path.join(os.getcwd(), filename)
dst = os.path.join(oos_dir, filename)
shutil.move(src, dst)
else:
os.remove(filename)
counter3 = 0
iter_scheduling_remote.sh 中修改的内容
echo "Simulation ..."
# run local >>>>>>>>>>>>>
mkdir random_vasp
mv sub-grade.sh random_vasp
mv Ap_inv* random_vasp
mv inv2 random_vasp
mv random_distance.py random_vasp
mv POSCAR random_vasp
cd random_vasp
python random_distance.py
2. 切割的方法
切割后的结构grade比较大，不太好用
"""
输入可以有多种形式,输出目前只有writ_poscar可以
"""
import numpy as np
import os
from ase.atoms import Atoms
from ase.io import read
import subprocess
import shutil
def write2my(file_path, ene_i, lat_i, ele_i, coo_i, foc_i, vir_i=None):
lat_i = lat_i.reshape(3, 3)
coo_i = coo_i.reshape(-1, 3)
foc_i = foc_i.reshape(-1, 3)
with open(file_path, 'w') as file:
file.write(f"# total energy = {ene_i} eV\n\n")
if vir_i is not None:
vir_i = np.array(vir_i).reshape(-1)
file.write(f"VIRIAL\n")
for i in vir_i:
file.write(f'{i:20.8f}')
file.write("\n")
file.write("CRYSTAL\n")
file.write("PRIMVEC\n")
for j in lat_i:
for k in j:
file.write(f'{k:20.8f}')
file.write('\n')
file.write("PRIMCOORD\n")
file.write(f"{len(coo_i)} 1\n")
for j in range(len(coo_i)):
file.write(f'{ele_i[j]:2}')
# coo
for k in coo_i[j]:
file.write(f"{k:20.8f}")
# force
for k in foc_i[j]:
file.write(f"{k:20.8f}")
file.write("\n")
pass
def read_from_ase_atoms(atoms: Atoms):
try:
ene = atoms.get_potential_energy()
except:
ene = 0.0
lat = atoms.get_cell()
pos = atoms.get_positions()
try:
foc = atoms.get_forces()
except:
foc = np.zeros_like(pos)
try:
sts = atoms.get_stress()
xx, yy, zz, yz, xz, xy = - sts * atoms.get_volume()
vir = np.array(
[[xx, xy, xz],
[xy, yy, yz],
[xz, yz, zz]]).reshape(-1)
except:
vir = np.zeros([3, 3])
ele = atoms.get_chemical_symbols()
return ene, lat, pos, foc, vir, ele
def read_poscar(filename):
with open(filename,'r') as file:
lines =file.readlines()
title =lines[0].strip()
scale =float(lines[1].strip())
#for i in range(2,5):
# lattice_list.append(list(map(float,line[i].split())))
#lattice_array=np.array(lattice_list)
lattice_vectors = np.array([list(map(float,lines[i].split())) for i in range(2,5)])
element_type=lines[5].split()
atomnumbers=np.array(list(map(int,lines[6].split())))
atompositions=np.array([list(map(float,lines[i].split()[:3])) for i in range(8,8+sum(atomnumbers))])
cor_type = lines[7].strip()
if cor_type.lower() == "direct" or cor_type.lower() == "d":
atompositions[atompositions < 0]+=1
atompositions[atompositions >=1]-=1
elif cor_type.lower() == "cartesian" or cor_type.lower() == "c":
pass
return(title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions)
def write_poscar(filename,title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions):
with open(filename,'w') as file:
file.write(f"{title}\n")
file.write(f"{scale}\n")
for vec in lattice_vectors:
file.write(f"{vec[0]:.8f} {vec[1]:.8f} {vec[2]:.8f}\n")
#" ".join用于将可迭代对象的元素以空格连接起来
file.write(" ".join(element_type)+ "\n")
file.write(" ".join(map(str,atomnumbers))+ "\n")
file.write(f"{cor_type}\n")
for pos in atompositions:
file.write(f"{pos[0]:.8f} {pos[1]:.8f} {pos[2]:.8f}\n")
#输出不同类型原子的坐标
def atom_pos(atomnumbers,atompositions):
atom_pos = np.split(atompositions,np.cumsum(atomnumbers)[:-1])
return atom_pos
"""
功能区
"""
#计算体积
def cal_volume(lattice_vectors):
volume = np.abs(np.dot(lattice_vectors[0],np.cross(lattice_vectors[1],lattice_vectors[2])))
return volume
#转换分数坐标为笛卡尔坐标
def change_dir_to_car(filename):
title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions= read_poscar(filename)
if cor_type.lower() == "cartesian" or cor_type.lower() == "c":
pass
elif cor_type.lower() == "direct" or cor_type.lower() == "d":
atompositions = np.dot(atompositions,lattice_vectors)
cor_type = "Cartesian"
return title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions
#笛卡尔坐标转换为分数坐标
def change_car_to_dir(filename):
title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions = read_poscar(filename)
if cor_type.lower() == "cartesian" or cor_type.lower() == "c":
atompositions = np.dot(atompositions,np.linalg.inv(lattice_vectors))
cor_type = "Direct"
atompositions[atompositions < 0]+=1
atompositions[atompositions >=1]-=1
elif cor_type.lower() == "direct" or cor_type.lower() == "d":
pass
return title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions
#晶格向量的变换
def lattice_change(filename,cell_type):
title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions= change_car_to_dir(filename)
if cell_type == '211':
pass
elif cell_type == "121":
middle_value = lattice_vectors[0].copy()
lattice_vectors[0] = lattice_vectors[1]
lattice_vectors[1] = lattice_vectors[2]
lattice_vectors[2] = middle_value
middle_value2 = atompositions[:,0].copy()
atompositions[:,0] = atompositions[:,1]
atompositions[:,1] = atompositions[:,2]
atompositions[:,2] = middle_value2
elif cell_type == "112":
middle_value = lattice_vectors[0].copy()
lattice_vectors[0] = lattice_vectors[2]
lattice_vectors[2] = lattice_vectors[1]
lattice_vectors[1] = middle_value
middle_value2 = atompositions[:,0].copy
atompositions[:,0] = atompositions[:,2]
atompositions[:,2] = atompositions[:,1]
atompositions[:,1] = middle_value2
return title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions
#将211晶胞分解为两个111晶胞
def split_211(input_file,cell_type,output_file1,output_file2):
title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions= lattice_change(input_file,cell_type)
atom_type_pos = atom_pos(atomnumbers,atompositions)
total_atoms1 = []
total_atoms_number1 = []
total_atoms2 = []
total_atoms_number2 = []
if cor_type.lower() == "direct" or cor_type.lower() == "d":
for sub_atom_type_pos in atom_type_pos:
print(sub_atom_type_pos)
coord_1=sub_atom_type_pos[sub_atom_type_pos[:,0]<=0.48].copy()
total_atoms1.append(coord_1)
total_atoms_number1.append(len(coord_1))
coord_2=sub_atom_type_pos[sub_atom_type_pos[:,0]>0.52].copy()
total_atoms2.append(coord_2)
total_atoms_number2.append(len(coord_2))
# elif cor_type.lower() == "cartesian" or cor_type.lower() == "c":
# for sub_atom_type_pos in atom_type_pos:
# sub_atom_type_pos=np.dot(sub_atom_type_pos,np.linalg.inv(lattice_vectors))
# sub_atom_type_pos[sub_atom_type_pos < 0]+=1
# sub_atom_type_pos[sub_atom_type_pos >=1]-=1
# coord_1=sub_atom_type_pos[sub_atom_type_pos[:,0]<=0.5].copy()
# total_atoms1.append(coord_1)
# total_atoms_number1.append(len(coord_1))
# coord_2=sub_atom_type_pos[sub_atom_type_pos[:,0]>0.5].copy()
# total_atoms2.append(coord_2)
# total_atoms_number2.append(len(coord_2))
# cor_type = "Direct"
total_atoms1=np.vstack(total_atoms1)
total_atoms2=np.vstack(total_atoms2)
total_atoms1=np.dot(total_atoms1,lattice_vectors)
total_atoms2=np.dot(total_atoms2,lattice_vectors)
lattice_vectors[0]/=2
cor_type ="Cartesian"
# print(total_atoms1)
# print(total_atoms_number1)
write_poscar(output_file1,title,scale,lattice_vectors,element_type,total_atoms_number1,cor_type,total_atoms1)
write_poscar(output_file2,title,scale,lattice_vectors,element_type,total_atoms_number2,cor_type,total_atoms2)
#POSCAR坐标的修改
def coords_change(axis,value,cor_type,atompositions):
if cor_type.lower() == "Cartesian" or cor_type.lower() == "c":
print("use change_car_to_dir to input")
elif cor_type.lower() == "direct" or cor_type.lower() == "d":
if axis == 'x':
atompositions[:,0] += value
elif axis == 'y':
atompositions[:,1] += value
elif axis == 'z':
atompositions[:,2] += value
else:
raise ValueError(f"invalid axis")
else:
raise ValueError(f"invalid cor_type")
return atompositions
#切割x*y*z倍的晶胞为x*y*z个晶胞
def split_cell(input_file,x,y,z):
input_file_without_ext = os.path.splitext(input_file)[0]
title,scale,lattice_vectors,element_type,atomnumbers,cor_type,atompositions= change_car_to_dir(input_file)
atom_type_pos = atom_pos(atomnumbers,atompositions)
total_cells = [[[{elem:[] for elem in element_type} for _ in range(z)] for _ in range(y)] for _ in range(x)]
total_atom_numbers = [[[{elem: 0 for elem in element_type} for _ in range(z)] for _ in range(y)] for _ in range(x)]
if cor_type.lower() == "direct" or cor_type.lower == "d":
for elem_index,sub_atom_type_pos in enumerate(atom_type_pos):
elem = element_type[elem_index]
for atom in sub_atom_type_pos:
ix = int(atom[0]*x)
iy = int(atom[1]*y)
iz = int(atom[2]*z)
total_cells[ix][iy][iz][elem].append(atom)
total_atom_numbers[ix][iy][iz][elem] += 1
cor_type ="Cartesian"
lattice_vectors_changed = np.zeros((3,3))
lattice_vectors_changed[0] = lattice_vectors[0]/x
lattice_vectors_changed[1] = lattice_vectors[1]/y
lattice_vectors_changed[2] = lattice_vectors[2]/z
os.makedirs("xsf_output")
for ix in range(x):
for iy in range(y):
for iz in range(z):
output_file =f"{input_file_without_ext}_{ix}{iy}{iz}.vasp"
sub_atoms = np.vstack([atom for elem in element_type for atom in total_cells[ix][iy][iz][elem]])
sub_atoms = np.dot(sub_atoms, lattice_vectors)
subcell_numbers = [total_atom_numbers[ix][iy][iz][elem] for elem in element_type]
write_poscar(output_file,title,scale,lattice_vectors_changed,element_type,subcell_numbers,cor_type,sub_atoms)
datas = read(output_file, index=":")
ene, lat, pos, foc, vir, ele = read_from_ase_atoms(datas[0])
output_file_path = os.path.join("xsf_output", f"{output_file}.xsf")
write2my(
output_file_path,
ene_i=ene,
lat_i=lat,
ele_i=ele,
coo_i=pos,
foc_i=foc,
vir_i=vir)
split_cell("POSCAR",2,2,2)
#os.remove("POS")
subprocess.run(["dsub","-s","sub-grade.sh"], check=True)
counter = 0
with open("grade-log", 'r') as file:
lines = file.readlines()
for line in lines:
if 'grad:' in line:
parts = line.split()
filename = parts[0]
filename_first = os.path.splitext(filename)[0]
grade_value = float(parts[-1])
if 3 < grade_value < 8:
os.makedirs(filename_first)
shutil.move(f"{filename_first}", os.path.join (filename_first, "POSCAR"))
else:
os.remove(f"{filename_first}")
转载请注明来源 有问题可通过github提交issue