脚本合集

# shell脚本
##  1.交批作业
 ```shell
 #!/bin/sh
 #PBS -N cp2k-test
 #PBS -l nodes=1:ppn=32,walltime=999:00:00
 ##PBS -q workq 
 pwd
 cd $PBS_O_WORKDIR
 echo $PBS_O_WORKDIR
 workdir=$PWD
 workpath=/home/zhou/work/haq/cp2k/bazrs/ba3zr2s7_scan/1x2/hole/diffnew
 
 source /opt/intel/bin/compilervars.sh intel64
 source /opt/intel/impi/2018.4.274/bin64/mpivars.sh intel64
 export PATH=$PATH:/opt/
 
 CP2K_EXEC="mpirun -n 32 cp2k_centos.7.6.1810_intel18.0.0_glibc.2.17-260.popt cp2k.inp 1>cp2k.out 2>cp2k.err"
 
 START=300 # start configuration
 END=1200 # end configuration
 STEP=300
 #for i in `seq ${START} ${END}`
 for((i=60;i<3001;i+=60));
 do
     ii=`printf "%d" $i`
         mkdir $workpath/run/${ii}
         cp $workpath/cp2k.inp $workpath/run/${ii}/cp2k.inp
         cp $workpath/${ii}.cif $workpath/run/${ii}/cell.cif
         cp $workpath/GTH-SCAN-POTENTIAL $workpath/BASIS_MOLOPT $workpath/rVV10_kernel_table.dat $workpath/run/${ii}
         cp $workpath/run.std $workpath/run/${ii}/run
         cd $workpath/run/${ii}
         qsub run
         cd ../..
 done
 
 rm -rf /tmp/nodefile.$$
 rm -rf /tmp/nodes.$$
 
 
 ```
## 2.批量下载
 ```shell
 
 workpath=/home/zhou/work/haq/cp2k/bazrs/ba3zr2s7_scan/1x2/hole/diffnew
 for((i=60;i<3001;i+=60));
 do
     ii=`printf "%d" $i`
         cd $workpath/run/${ii}
         mv BaZrS_aimd-cube-ELECTRON_DENSITY-1_0.cube ${ii}_0.cube
         sz ${ii}_0.cube
         rm *RESTART.wfn*
         cd ../..
 done
 
 
 ```
## Python脚本
#### 1.读取一帧xyz文件计算平均键长
```python
from ase.io import read
from ase.neighborlist import NeighborList

# 读取xyz文件
atoms = read('fs.xyz')

# 设置键长的阈值
cutoff = 3.0

# 创建NeighborList对象
nl = NeighborList([cutoff / 2] * len(atoms), self_interaction=False, bothways=True)

# 构建邻居列表
nl.update(atoms)

# 初始化计数器和键长之和
count = 0
bond_length_sum = 0.0

# 遍历所有原子对，计算键长
for i, atom in enumerate(atoms):
    for j in nl.get_neighbors(i)[0]:
        neighbor = atoms[j]
        if atom.symbol == 'Zr' and neighbor.symbol == 'S':
            bond_length = atoms.get_distance(i, j)
            count += 1
            bond_length_sum += bond_length

# 计算平均键长
if count > 0:
    average_bond_length = bond_length_sum / count
    print(f'Average Zr-S bond length: {average_bond_length:.3f} Å')
else:
    print('No Zr-S bonds found.')
```
#### 2.读取一帧xyz文件计算平均键长(只计算某一高度范围内的原子)
```python
from ase.io import read
from ase.neighborlist import NeighborList

# 读取xyz文件
atoms = read('fs.xyz')

# 设置键长的阈值
cutoff = 3.0

# 创建NeighborList对象
nl = NeighborList([cutoff / 2] * len(atoms), self_interaction=False, bothways=True)

# 构建邻居列表
nl.update(atoms)

# 初始化计数器和键长之和
count = 0
bond_length_sum = 0.0

# 遍历所有原子对，计算键长
for i, atom in enumerate(atoms):
    # 限制z坐标在8-12之间的原子对
    if 1 <= atom.position[2] <= 4:
        for j in nl.get_neighbors(i)[0]:
            neighbor = atoms[j]
            if atom.symbol == 'Zr' and neighbor.symbol == 'S':
                bond_length = atoms.get_distance(i, j)
                count += 1
                bond_length_sum += bond_length

# 计算平均键长
if count > 0:
    average_bond_length = bond_length_sum / count
    print(f'Average Zr-S bond length: {average_bond_length:.3f} Å')
else:
    print('No Zr-S bonds found.')

```
## 3.读取多帧xyz文件计算平均键长

```python
import sys
import numpy as np

# 定义计算键长函数
def calc_distance(coords, atom1, atom2):
    # 计算向量
    v = coords[atom1] - coords[atom2]
    # 计算键长
    distance = np.linalg.norm(v)
    return distance

# 读取XYZ文件
filename = sys.argv[1]
with open(filename, 'r') as f:
    lines = f.readlines()
# 解析XYZ文件
num_atoms = int(lines[0])
num_frames = int(len(lines[2:]) / (num_atoms + 2)) + 1
coords = np.zeros((num_frames, num_atoms, 3))
for i in range(num_frames):
    start = i * (num_atoms + 2) + 2
    end = start + num_atoms
    for j in range(start, end):
        line = lines[j].split()
        coords[i, j-start] = np.array([float(line[1]), float(line[2]), float(line[3])])

# 找到所有Zr-S键的键长
zr_indices = []  # 存储Zr原子的索引
s_indices = []  # 存储S原子的索引
for i in range(num_atoms):
    if lines[i+2].split()[0] == 'Zr':
        zr_indices.append(i)
    elif lines[i+2].split()[0] == 'S':
        s_indices.append(i)

# 计算每一帧结构的平均键长（小于3Å）
avg_distances = np.zeros(num_frames)
for k in range(num_frames):
    sum_distance = 0.0
    count = 0
    for i in zr_indices:
        for j in s_indices:
            if i != j:
                distance = calc_distance(coords[k], i, j)
                if distance < 3.0:
                    sum_distance += distance
                    count += 1
    if count > 0:
        avg_distances[k] = sum_distance / count

# 将结果写入文件
with open('zr_s_avg_distances.dat', 'w') as f:
    for i in range(num_frames):
        f.write('{} {}\n'.format(i, avg_distances[i]))
```
## 4.读取cif文件计算键长
```python
from ase.io import read
from ase.build import make_supercell
from ase.neighborlist import NeighborList
import numpy as np

# 读取cif文件并创建超胞
structure = read('cell.cif')
supercell = make_supercell(structure, [[2, 0, 0], [0, 2, 0], [0, 0, 2]])

# 初始化计数器和键长之和
count = 0
bond_length_sum = 0.0

# 设置截断半径
cutoff = 3.0

# 计算邻居列表
cutoffs = [cutoff] * len(supercell)
nl = NeighborList(cutoffs, self_interaction=False, bothways=True)
nl.update(supercell)

# 遍历所有原子对，计算键长
for i, site1 in enumerate(supercell):
    for j, site2 in enumerate(supercell):
        if i >= j:
            continue
        if site1.symbol == 'Zr' and site2.symbol == 'S':
            indices, offsets = nl.get_neighbors(i)
            if j in indices:
                index = np.where(indices == j)[0][0]
                distance = np.linalg.norm(supercell.positions[i] - supercell.positions[j] + np.dot(offsets[index], supercell.get_cell()))
                if distance <= cutoff:
                    bond_length_sum += distance
                    count += 1

# 计算MD结果的平均键长
if count > 0:
    average_bond_length = bond_length_sum / count
    print(f'Average Zr-S bond length: {average_bond_length:.3f} Å')
else:
    print('No Zr-S bonds found.')
```
## 5.读取cif文件计算键长（只计算某一高度范围内的原子）
```python

from ase.io import read
from ase.build import make_supercell
from ase.neighborlist import NeighborList
import numpy as np

# 读取cif文件并创建超胞
structure = read('cell.cif')
supercell = make_supercell(structure, [[2, 0, 0], [0, 2, 0], [0, 0, 2]])

# 初始化计数器和键长之和
count = 0
bond_length_sum = 0.0

# 设置截断半径
cutoff = 3.0

# 计算邻居列表
cutoffs = [cutoff] * len(supercell)
nl = NeighborList(cutoffs, self_interaction=False, bothways=True)
nl.update(supercell)

# 遍历所有原子对，计算键长
for i, site1 in enumerate(supercell):
    if 1 <= site1.position[2] <= 4:
        for j, site2 in enumerate(supercell):
            if i >= j:
                continue
            if site1.symbol == 'Zr' and site2.symbol == 'S':
                indices, offsets = nl.get_neighbors(i)
                if j in indices:
                    index = np.where(indices == j)[0][0]
                    distance = np.linalg.norm(supercell.positions[i] - supercell.positions[j] + np.dot(offsets[index], supercell.get_cell()))
                    if distance <= cutoff:
                        bond_length_sum += distance
                        count += 1

# 计算MD结果的平均键长
if count > 0:
    average_bond_length = bond_length_sum / count
    print(f'Average Zr-S bond length: {average_bond_length:.3f} Å')
else:
    print('No Zr-S bonds found.')

```

## 6.计算MD每一帧的键角
执行命令`python3 **.py **.xyz 156 114 126`

![](/media/202306/2023-06-12_113752_6537810.2947825757748309.png)
```python
import sys
import numpy as np

# 读取命令行参数
atom1 = int(sys.argv[2]) - 1
atom2 = int(sys.argv[3]) - 1
atom3 = int(sys.argv[4]) - 1

# 定义键角计算函数
def calc_angle(coords):
    # 计算向量
    v1 = coords[atom1] - coords[atom2]
    v2 = coords[atom3] - coords[atom2]
    # 计算夹角
    cos_angle = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2))
    angle = np.arccos(cos_angle) * 180 / np.pi
    # 判断夹角的正负
 #   cross = np.cross(v1, v2)
 #   if np.dot(cross, coords[atom2]) < 0:
    #    angle = -angle
    return angle

# 读取XYZ文件
filename = sys.argv[1]
with open(filename, 'r') as f:
    lines = f.readlines()

# 解析XYZ文件
num_atoms = int(lines[0])
num_frames = int(len(lines[2:]) / (num_atoms + 2)) + 1
coords = np.zeros((num_frames, num_atoms, 3))
for i in range(num_frames):
    start = i * (num_atoms + 2) + 2
    end = start + num_atoms
    for j in range(start, end):
        line = lines[j].split()
        coords[i, j-start] = np.array([float(line[1]), float(line[2]), float(line[3])])

# 计算键角
angles = np.zeros(num_frames)
for i in range(num_frames):
    angles[i] = calc_angle(coords[i])

# 将结果写入文件
with open('angles.dat', 'w') as f:
    for i in range(num_frames):
        f.write('{} {}\n'.format(i, angles[i]))
```

## 7.把XDATCATR文件转变为xyz格式
（xyz文件可用MS打开，没有周期性，可以放进VMD加周期性，XDATCAR可用OVITO打开）
```python
#lipai@mail.ustc.edu.cn
#convert XDATCAR to unwraped xyz file
import numpy as np
from copy import deepcopy

xdatcar = open('XDATCAR', 'r')
xyz = open('XDATCAR.xyz', 'w')

system = xdatcar.readline()
scale = float(xdatcar.readline().rstrip('\n'))
print(scale)

#get lattice vectors
a1 = np.array([ float(s)*scale for s in xdatcar.readline().rstrip('\n').split() ])
a2 = np.array([ float(s)*scale for s in xdatcar.readline().rstrip('\n').split() ])
a3 = np.array([ float(s)*scale for s in xdatcar.readline().rstrip('\n').split() ])

comment='Lattice=\"'+str(a1[0])+' '+str(a1[1])+' '+str(a1[2])
comment=comment+str(a2[0])+' '+str(a2[1])+' '+str(a2[2])
comment=comment+str(a3[0])+' '+str(a3[1])+' '+str(a3[2])+'\"'

#Read xdatcar
element_names = xdatcar.readline().rstrip('\n').split()

element_dict = {}
element_numbers = xdatcar.readline().rstrip('\n').split()

Natom = 0
Ntype = len(element_names)
Nname=[]
for t in range(Ntype):
    Natom += int(element_numbers[t])
    for i in range(int(element_numbers[t])):
        Nname.append(element_names[t])

print(Ntype,Natom)

f_prev=np.zeros([Natom,3])
f_next=np.zeros([Natom,3])
while True:
    line = xdatcar.readline()
    if len(line) == 0:
        break
    xyz.write(str(Natom) + "\n"+comment+"\n")

for atom in range(Natom):
        p = xdatcar.readline().rstrip('\n').split()
        f_next[atom,:] = np.array([ float(s) for s in p ])
        for x in range(3):
            if(f_next[atom,x]-f_prev[atom,x]<-0.5):
                f_next[atom,x]+=1
            elif(f_next[atom,x]-f_prev[atom,x]>0.5):
                f_next[atom,x]-=1

c_coords=f_next[atom,0]*a1+f_next[atom,1]*a2+f_next[atom,2]*a3
        xyz.write(Nname[atom]+" "+str(c_coords[0])+" "+str(c_coords[1])+" "+str(c_coords[2])+"\n")

f_prev=deepcopy(f_next)

xdatcar.close()
xyz.close()
```
## 8.从OUTCAR中提取轨迹文件
虽然说有XDATCAR,但是打开了机器学习的话可能输出的结构不全

```
def find_positions(outcar_lines):
    positions = []
    for idx, line in enumerate(outcar_lines):
        if "POSITION                                       TOTAL-FORCE" in line:
            # 计算原子坐标的最后一行
            last_coord_line = idx + 2 + sum(count for _, count in elements_and_counts)
            # 检查原子坐标最后一行下面的三十行内是否包含 "--------------------------------------- Ionic step"
            for i in range(1, 31):  # 检查下面的30行
                if last_coord_line + i < len(outcar_lines) and "--------------------------------------- Ionic step" in outcar_lines[last_coord_line + i]:
                    positions.append(idx)
                    break  # 找到后立即退出循环
    return positions

def extract_structure(outcar_lines, start_idx, total_atoms):
    start_data = start_idx + 2
    data = outcar_lines[start_data : start_data + total_atoms]
    return data

def process_atoms(data, element_counts):
    atoms = []
    index = 0
    for element, count in element_counts:
        for _ in range(count):
            if index >= len(data):
                break
            coord = data[index].split()[:3]  # 只提取前三列坐标
            atoms.append(f"{element} {' '.join(coord)}")
            index += 1
    return atoms

def main(outcar_path, elements_and_counts):
    total_atoms = sum(count for _, count in elements_and_counts)
    with open(outcar_path, 'r') as f:
        outcar_lines = f.readlines()
    positions = find_positions(outcar_lines)
    with open('XDATCAR_MD', 'w') as output_file:
        for idx, pos in enumerate(positions):
            data = extract_structure(outcar_lines, pos, total_atoms)
            if len(data) < total_atoms:
                output_file.write(f"Frame {idx + 1}: Incomplete data, skipping.\n")
                continue
            atoms = process_atoms(data, elements_and_counts)
            output_file.write(f"Frame {idx + 1}:\n")
            for atom in atoms:
                output_file.write(f"  {atom}\n")
            output_file.write("\n")  # 添加空行以分隔帧

if __name__ == "__main__":
    import sys
    outcar_path = 'OUTCAR'  # 根据实际情况修改
    elements_and_counts = [('Sn', 49), ('F', 98), ('O', 49)]  # 根据实际情况修改
    main(outcar_path, elements_and_counts)

```
## 9.从轨迹文件中读取特定原子对的x方向差值的平均
承接上一个代码得到的轨迹文件，求的M与O原子的x方向的差值，可以作为极性/偶极矩等的指标

```
import numpy as np

def readlatticeconstant(outcarpath):
    with open(outcarpath, 'r') as f:
        for line in f:
            if "length of vectors" in line:
                next_line = next(f).strip()  # 跳过当前行并读取下一行
                parts = next_line.split()  # 将行分割成列表
                if len(parts) > 1:  # 确保列表至少有两个元素
                    lattice_constant_b = float(parts[1])  # 读取下一行的第二个数字
                    return lattice_constant_b
                else:
                    raise ValueError(f"无法从行中读取晶格常数: '{next_line}'")
    raise ValueError("未能从OUTCAR文件中读取晶格常数b")

def parsetext(textpath,elements):
    frames = []
    with open(textpath, 'r') as f:
        lines = f.readlines()
        frame = []
        for line in lines:
            if "Frame" in line:
                if frame:
                    frames.append(frame)
                frame = []
            elif any(element in line for element in elements):
                parts = line.split()
                atom = [parts[0], float(parts[1]), float(parts[2]), float(parts[3])]
                frame.append(atom)
        if frame:
            frames.append(frame)
    return frames

def findnearestoxygen(M_atoms, O_atoms):
    x_diffs = []
    for M_atom in M_atoms:
        M_y = M_atom[2]
        min_distance = float('inf')
        nearest_O = None
        for O_atom in O_atoms:
            O_y = O_atom[2]
            if O_y < M_y:
                distance = np.sqrt((M_atom[1] - O_atom[1])**2 + (M_atom[2] - O_atom[2])**2 + (M_atom[3] - O_atom[3])**2)
                if distance < min_distance:
                    min_distance = distance
                    nearest_O = O_atom
        if nearest_O:
            x_diffs.append(M_atom[1] - nearest_O[1])
    return x_diffs

def calculatexdifferenceaverage(frames, Melement, Oelement, lattice_constant):
    all_frame_averages = []
    output_data = []
    for i, frame in enumerate(frames):
        M_atoms = [atom for atom in frame if atom[0] == Melement]
        O_atoms = [atom for atom in frame if atom[0] == Oelement]
        max_M_y = max([atom[2] for atom in M_atoms])
        O_atoms = [[atom[0], atom[1], atom[2] - lattice_constant if atom[2] > max_M_y else atom[2], atom[3]] for atom in O_atoms]
        x_diffs = findnearestoxygen(M_atoms, O_atoms)
        if x_diffs:
            frame_average = np.mean(x_diffs)
            all_frame_averages.append(frame_average)
            output_data.append(f"Frame {i}: Average x difference = {frame_average}\n")
    if all_frame_averages:
        overall_average = np.mean(all_frame_averages)
        output_data.append(f"Overall average x difference: {overall_average}\n")
        return overall_average, output_data
    else:
        return None, output_data

if __name__ == "__main__":
    textpath = 'XDATCAR_MD'  # 根据实际情况修改
    outcarpath = 'OUTCAR'  # 根据实际情况修改
    Melement = 'Sn'  # 根据实际情况修改
    Oelement = 'O'  # 根据实际情况修改
    outputpath = 'output.txt'  # 输出文件路径
    lattice_constant = readlatticeconstant(outcarpath)
    frames = parsetext(textpath, [Melement, Oelement])
    overall_average, output_data = calculatexdifferenceaverage(frames, Melement, Oelement, lattice_constant)
    with open(outputpath, 'w') as f:
        f.write(f"读取的晶格常数b: {lattice_constant}\n")
        f.writelines(output_data)

```

脚本合集

【温馨提示：本站文档可配置可见范围，如登录后可见、对特定群组可见等，看不到就是没权限】

注册码获取邮箱
work@liuyaoze.com

脚本合集

【温馨提示：本站文档可配置可见范围，如登录后可见、对特定群组可见等，看不到就是没权限】

注册码获取邮箱 work@liuyaoze.com

注册码获取邮箱
work@liuyaoze.com