读取cp2k分子动力学轨迹自动生成输入文件脚本

准备*-pos-1.xyz文件，和脚本放在同一目录，运行python pos-syz.py

import os
import re

# 设置时间步数，例如MD跑了5000步，step设置4000，那么将会生成1001到5000的轨迹，命名为0001fs到4000fs
step = 100

# 自洽的inp模板
def inp_template(atom_position):
    """自洽的inp模板，需要根据自己的体系手动修改"""
    return f"""
#Generated by Multiwfn
&GLOBAL
  PROJECT MoS2
  PRINT_LEVEL LOW
  RUN_TYPE ENERGY
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep
  &SUBSYS
    &CELL
      A     9.54729600     0.00000000     0.00000000
      B    -4.77363588     8.26820787     0.00000000
      C     0.00000000     0.00000000    19.98652600
      PERIODIC XYZ #Direction(s) of applied PBC (geometry aspect)
    &END CELL
    &COORD
{atom_position}
    &END COORD
    &KIND S
      ELEMENT S
      BASIS_SET DZVP-MOLOPT-SR-GTH-q6
      POTENTIAL GTH-PBE
    &END KIND
    &KIND Mo
      ELEMENT Mo
      BASIS_SET DZVP-MOLOPT-SR-GTH-q14
      POTENTIAL GTH-PBE
    &END KIND
  &END SUBSYS

  &DFT
  
    &PRINT
      &MO_CUBES
        WRITE_CUBE .TRUE.           ! 启用写Cube文件
        NHOMO 1                      ! 不自动输出HOMO轨道
        NLUMO 1                      ! 不自动输出LUMO轨道
        STRIDE 1 1 1                 ! Cube文件的分辨率，1 1 1 表示最高分辨率
      &END MO_CUBES
    &END PRINT
  
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  POTENTIAL
#   WFN_RESTART_FILE_NAME MoS2-RESTART.wfn
    CHARGE    0 #Net charge
    MULTIPLICITY    1 #Spin multiplicity
    &QS
      EPS_DEFAULT 1.0E-11 #Set all EPS_xxx to values such that the energy will be correct up to this value
    &END QS
    &POISSON
      PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
      PSOLVER PERIODIC #The way to solve Poisson equation
    &END POISSON
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
    &MGRID
      CUTOFF  350
      REL_CUTOFF  50
    &END MGRID
    &SCF
      MAX_SCF 128
      EPS_SCF 5.0E-06 #Convergence threshold of density matrix of inner SCF
#     SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
#     IGNORE_CONVERGENCE_FAILURE #Continue calculation even if SCF not converged, works for version >= 2024.1
      &DIAGONALIZATION
        ALGORITHM STANDARD #Algorithm for diagonalization
      &END DIAGONALIZATION
      &MIXING #How to mix old and new density matrices
        METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
        ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one
        NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme
      &END MIXING
      &PRINT
        &RESTART #Note: Use "&RESTART OFF" can prevent generating .wfn file
          BACKUP_COPIES 0 #Maximum number of backup copies of wfn file. 0 means never
        &END RESTART
      &END PRINT
    &END SCF
  &END DFT
&END FORCE_EVAL
"""

# 后面的一般不需要修改
########################################################################################
# 基础目录
base_dir = './run'

# 确保基础目录存在
if not os.path.exists(base_dir):
    os.makedirs(base_dir)

def extract_atom_number(file_path):
    """从文件中提取原子数量"""
    with open(file_path, 'r') as file:
        return int(file.readline().strip())

def extract_last_i(lines):
    """从文件行中提取最后一次出现的 "i =" 的值"""
    for line in reversed(lines):
        match = re.search(r"i =\s*(\d+)", line)
        if match:
            return int(match.group(1))
    return None

def read_atom_positions(file_path, start_line, atom_number):
    """
    读取指定起始行后的原子位置信息
    
    参数:
    file_path -- 文件路径
    start_line -- 起始行号
    atom_number -- 原子数量
    
    返回:
    atoms -- 原子位置信息的列表
    """
    atoms = []
    with open(file_path, 'r') as file:
        for current_line, line in enumerate(file, start=1):
            if current_line >= start_line and current_line < start_line + atom_number:
                atoms.append(line.strip())  
            elif current_line >= start_line + atom_number:
                break  
    return atoms

def create_folder(folder_name):
    """创建文件夹"""
    os.makedirs(folder_name, exist_ok=True)
    print(f"Folder '{folder_name}' created")


# 文件路径
file_path = './MoS2-pos-1.xyz'

# 读取文件内容
with open(file_path, 'r') as file:
    lines = file.readlines()

# 提取所需值
atom_number = extract_atom_number(file_path)
last_i_value = extract_last_i(lines)

# 打印提取的值
print(f"Last 'i' value: {last_i_value}, Atom Number: {atom_number}")


# 创建步骤文件夹
for i in range(1, step + 1):
    folder_name = os.path.join(base_dir, f"{i:04d}fs")
    create_folder(folder_name)

# 遍历每个文件夹
for i in range(1, step + 1):
    folder_name = os.path.join(base_dir, f"{i:04d}fs")
    # 创建文件路径
    file_path = os.path.join(folder_name, "cp2k.inp")
    print((last_i_value - step + i) * (atom_number + 2) + 3)
    atoms = read_atom_positions('./MoS2-pos-1.xyz', (last_i_value - step + i) * (atom_number + 2) + 3, atom_number)
    atom_position = "\n".join(atoms)
    # 使用特定内容替换文本中的占位符
    text_content = inp_template(atom_position)
    # 写入文件
    with open(file_path, "w") as file:
        file.write(text_content)
    print(f"File '{file_path}' has been created in '{folder_name}' with modified content.")