氢水建模

import numpy as np
from ase.build import molecule
from ase import Atoms
from ase.io import read
from ase.io.vasp import write_vasp

def rotation_matrix_from_euler(alpha, beta, gamma, degrees=True):
    """
    生成绕 x, y, z 轴依次旋转 alpha, beta, gamma 的旋转矩阵。
    Rotation applied as: first Rx(alpha), then Ry(beta), then Rz(gamma).
    返回一个 3x3 的 numpy 矩阵 R，使得 r' = R @ r (列向量)。
    """
    if degrees:
        a = np.deg2rad(alpha)
        b = np.deg2rad(beta)
        g = np.deg2rad(gamma)
    else:
        a, b, g = alpha, beta, gamma

    Rx = np.array([[1, 0, 0],
                   [0, np.cos(a), -np.sin(a)],
                   [0, np.sin(a),  np.cos(a)]])
    Ry = np.array([[ np.cos(b), 0, np.sin(b)],
                   [0, 1, 0],
                   [-np.sin(b), 0, np.cos(b)]])
    Rz = np.array([[np.cos(g), -np.sin(g), 0],
                   [np.sin(g),  np.cos(g), 0],
                   [0, 0, 1]])
    # 总旋转：先 Rx 然后 Ry 然后 Rz -> R = Rz @ Ry @ Rx
    R = Rz @ Ry @ Rx
    return R

def add_water(atoms, target_pos, orientation=(0, 0, 0),
              anchor='COM', use_fractional=True, cell=None):

    if cell is None:
        cell = atoms.get_cell()

    # 计算目标笛卡尔坐标
    if use_fractional:
        if cell is None or np.linalg.norm(cell) == 0:
            raise ValueError("没有可用的 cell 来把分数坐标转换为笛卡尔。传入 cell 或将 use_fractional=False。")
        target_cart = np.dot(np.asarray(target_pos), cell)  # frac @ cell
    else:
        target_cart = np.asarray(target_pos, dtype=float)

    # 生成水分子
    water = molecule("H2O")  # ASE 内置
    # 选择锚点并将其移动到原点
    if anchor.upper() == 'COM':
        ref = water.get_center_of_mass()
    elif anchor.upper() == 'O':
        syms = water.get_chemical_symbols()
        try:
            o_idx = syms.index('O')
            ref = water.positions[o_idx].copy()
        except ValueError:
            # 退回到 COM
            ref = water.get_center_of_mass()
    else:
        raise ValueError("anchor 必须是 'COM' 或 'O'。")

    # 把参考点移动到原点
    water.translate(-ref)

    # 应用旋转（关于原点）
    R = rotation_matrix_from_euler(*orientation, degrees=True)
    # positions 为 (N,3)，对每个原子应用 R: pos' = (R @ pos.T).T
    water.positions = (R @ water.positions.T).T

    # 平移到目标笛卡尔位置
    water.translate(target_cart)

    # 合并并返回
    new_atoms = atoms + water
    return new_atoms

from ase.io import read
prim = read('../../CONTCAR')
new_atoms = add_water(prim, target_pos=(0.5, 0.5, 0.75), orientation=(10, 20, 90))
from ase.io import write
write_vasp("0Vo/ori2/POSCAR", new_atoms, direct=True,sort=True)

import numpy as np
from ase.build import molecule
from ase import Atoms
from ase.io import write, read
from ase.io.vasp import write_vasp

def rotation_matrix_from_euler(alpha, beta, gamma, degrees=True):
    """生成绕 x, y, z 轴依次旋转的旋转矩阵。"""
    if degrees:
        a = np.deg2rad(alpha)
        b = np.deg2rad(beta)
        g = np.deg2rad(gamma)
    else:
        a, b, g = alpha, beta, gamma

    Rx = np.array([[1, 0, 0],
                   [0, np.cos(a), -np.sin(a)],
                   [0, np.sin(a),  np.cos(a)]])
    Ry = np.array([[ np.cos(b), 0, np.sin(b)],
                   [0, 1, 0],
                   [-np.sin(b), 0, np.cos(b)]])
    Rz = np.array([[np.cos(g), -np.sin(g), 0],
                   [np.sin(g),  np.cos(g), 0],
                   [0, 0, 1]])
    return Rz @ Ry @ Rx


def add_h2(atoms, target_pos, orientation=(0, 0, 0),
           anchor='COM', use_fractional=True, cell=None):
    """
    在 atoms 中添加一分子 H2。

    参数:
      atoms         : ASE Atoms 对象
      target_pos    : 目标位置。如果 use_fractional=True，传入分数坐标 (fx,fy,fz)，否则笛卡尔坐标 (Å)。
      orientation   : (alpha, beta, gamma)，绕 x,y,z 轴的欧拉角（度）。
      anchor        : 'COM' 或 'H1'。决定氢分子哪个点对齐到 target_pos。
                      'COM' -> 质心对齐； 'H1' -> 把第一个氢原子对齐到目标。
      use_fractional: 是否把 target_pos 当作分数坐标，默认 True。
      cell          : 可选晶胞。如果 None，则用 atoms.get_cell()。

    返回:
      new_atoms : 含有氢分子的 ASE Atoms 对象
    """
    if cell is None:
        cell = atoms.get_cell()

    # 计算目标笛卡尔坐标
    if use_fractional:
        target_cart = np.dot(np.asarray(target_pos), cell)  # frac -> cart
    else:
        target_cart = np.asarray(target_pos, dtype=float)

    # 生成氢分子
    h2 = molecule("H2")

    # 锚点
    if anchor.upper() == 'COM':
        ref = h2.get_center_of_mass()
    elif anchor.upper() == 'H1':
        ref = h2.positions[0].copy()
    else:
        raise ValueError("anchor 必须是 'COM' 或 'H1'。")

    # 移动锚点到原点
    h2.translate(-ref)

    # 旋转
    R = rotation_matrix_from_euler(*orientation, degrees=True)
    h2.positions = (R @ h2.positions.T).T

    # 平移到目标
    h2.translate(target_cart)

    # 合并
    new_atoms = atoms + h2
    return new_atoms


from ase.io import read
prim = read('../../CONTCAR')
new_atoms = add_h2(prim, target_pos=(0.5, 0.5, 0.25), orientation=(45, 45, 0))
from ase.io import write
write_vasp("ori2/POSCAR", new_atoms, direct=True,sort=True)