MD键长键角分析

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

PREC = Accurate
SMASS = -3   #micro canonical ensemble
LREAL = A   #real space
IBRION = 0   #molecular dynamics
NBLOCK = 1   #update XDATCAR every x steps
TEBEG = 300   #start temperature
TEEND = 300   #final temperature
ISIF = 2   #2, ions change; 3, shape and ions change
ISYM = 0   #symmetry
NSW = 1000   #max ionic steps
POTIM = 1   #step width scaling
EDIFFG = -0.01   #ionic relaxation
EDIFF = 1e-5   #electronic SC-loop
LCHARG = .F.   #not save CHGCAR CHG
LWAVE = .F.   #not save WAVECAR
ISMEAR = 0   #gaussian smearing
SIGMA = 0.05   #the width of the smearing in eV
ALGO = Normal   #electronic minimisation algorithm
IVDW = 12   #DFT-D3 method with Becke-Jonson damping, van der Waals

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

import os
import re
import math
import shutil

def print_intro():
    print(
        "Convert OUTCAR from MD simulation to multiple pfiles. Supports generating 0-9999 pfiles.\n\n"
        
        "Usage:\n"
        "Place OUTCAR and CONTCAR from MD simulation in the same directory.\n"
        "The script will generate multiple pfiles in this directory, matching the number of MD steps (NSW in INCAR).\n\n"
        
        'copyright by misaraty (misaraty@163.com) last update: 2025-06-26\n'
    )

def extract_nsw(outcar_path):
    # Extract NSW value from OUTCAR
    with open(outcar_path, 'r') as f:
        for line in f:
            if '   NSW' in line:
                parts = line.split()
                return int(parts[2])
    raise ValueError("NSW keyword not found in OUTCAR.")

def get_atom_count(contcar_path):
    # Calculate total number of atoms from CONTCAR
    with open(contcar_path, 'r') as f:
        lines = f.readlines()
    atom_counts = list(map(int, lines[6].split()))
    return sum(atom_counts)

def extract_forces(outcar_path, atom_count):
    # Extract atomic positions from all MD steps in OUTCAR
    forces = []
    with open(outcar_path, 'r') as f:
        lines = f.readlines()

    for idx, line in enumerate(lines):
        if 'POSITION                                       TOTAL-FORCE' in line:
            block_start = idx + 2
            block = lines[block_start:block_start + atom_count]
            for atom_line in block:
                coords = atom_line.split()[:3]
                forces.append("  " + "  ".join(coords))
    return forces

def write_pfiles(nsw, atom_count, contcar_path, forces, output_dir="pfiles"):
    # Write pfiles based on extracted coordinates
    if os.path.exists(output_dir):
        shutil.rmtree(output_dir)
    os.makedirs(output_dir)

    # Read first 7 lines of CONTCAR (lattice & atom info)
    with open(contcar_path, 'r') as f:
        contcar_lines = f.readlines()[:7]

    for step in range(nsw):
        filename = os.path.join(output_dir, f"p{step+1:04d}")
        with open(filename, 'w') as f:
            f.writelines(contcar_lines)
            f.write("C\n")  # atom label (can be adjusted if necessary)
            start = step * atom_count
            end = start + atom_count
            f.writelines(line + '\n' for line in forces[start:end])

def main():
    # Set working directory to script location and execute conversion
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    outcar = 'OUTCAR'
    contcar = 'CONTCAR'

    print_intro()
    nsw = extract_nsw(outcar)
    atom_count = get_atom_count(contcar)
    forces = extract_forces(outcar, atom_count)
    write_pfiles(nsw, atom_count, contcar, forces)
    print(f"{nsw} pfiles have been created.\n")

if __name__ == "__main__":
    main()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

import os
import re
import shutil

def print_intro():
    print(
        "Convert XDATCAR to multiple pfiles\n\n"
        
        "Usage:\n"
        "Place the XDATCAR file in the current directory.\n"
        "This script will automatically create a 'pfiles' directory and output one file per step.\n\n"
        
        'copyright by misaraty (misaraty@163.com) last update: 2025-06-26\n'
    )

def read_header_lines(xdatcar_path):
    # Read the first 7 lines from XDATCAR: system name, scaling factor, 3x3 lattice, element symbols, and atom counts
    with open(xdatcar_path, 'r') as f:
        header = [next(f) for _ in range(7)]
    return header

def count_atoms(header_lines):
    # Count the total number of atoms from the 7th line of the header
    atom_counts = list(map(int, header_lines[6].split()))
    return sum(atom_counts)

def parse_xdatcar_blocks(xdatcar_path, atom_count):
    """
    Parse XDATCAR and return a list of snapshots.
    Each snapshot is a list of atomic coordinates for one MD step.
    """
    snapshots = []
    with open(xdatcar_path, 'r') as f:
        lines = f.readlines()

    config_indices = [i for i, line in enumerate(lines) if line.startswith('Direct configuration=')]

    for idx in range(len(config_indices)):
        start = config_indices[idx] + 1
        end = start + atom_count
        coords = lines[start:end]
        coords = ['  '.join(line.strip().split()) + '\n' for line in coords]
        snapshots.append(coords)

    return snapshots

def write_pfiles(header_lines, snapshots, output_dir='pfiles'):
    # Write each snapshot to a separate file in the specified output directory
    if os.path.exists(output_dir):
        shutil.rmtree(output_dir)
    os.makedirs(output_dir)

    for i, coords in enumerate(snapshots, 1):
        fname = os.path.join(output_dir, f"p{i:04d}")
        with open(fname, 'w') as f:
            f.writelines(header_lines)
            f.write("Direct\n")
            f.writelines(coords)

def main():
    print_intro()
    xdatcar = 'XDATCAR'
    if not os.path.exists(xdatcar):
        raise FileNotFoundError("XDATCAR file not found in the current directory.")

    header = read_header_lines(xdatcar)
    atom_count = count_atoms(header)
    snapshots = parse_xdatcar_blocks(xdatcar, atom_count)
    write_pfiles(header, snapshots)
    print(f"Generated {len(snapshots)} pfiles in the './pfiles' directory.")

if __name__ == '__main__':
    main()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

import os
import numpy as np
import warnings
warnings.filterwarnings('ignore')

# Set current working directory to the script location
os.chdir(os.path.dirname(os.path.abspath(__file__)))
print('copyright by misaraty (misaraty@163.com)\n' + 'last update: 2025-06-26\n')

# Parameters
num = 1000  # Number of files
name = [[2, 4], [2, 5]]  # Atom pairs (indices start from 1)

# Read lattice information and coordinate type from p0001
with open('./pfiles/p0001') as f:
    lines = f.readlines()

# Read 3×3 lattice matrix
lattice = np.array([
    list(map(float, lines[2].split())),
    list(map(float, lines[3].split())),
    list(map(float, lines[4].split()))
])

# Calculate inverse lattice matrix
inv_lattice = np.linalg.inv(lattice)

# Check if coordinates are in Direct mode
is_direct = lines[7].strip().startswith('Direct')

# Total number of atoms
natoms = sum(map(int, lines[6].split()))
coord_start = 8  # Starting line of coordinates

def get_cart_coords(filepath):
    # Read coordinates from pfile and convert to Cartesian if needed
    with open(filepath) as f:
        lines = f.readlines()
    coords = np.array([
        list(map(float, lines[coord_start + i].split()))
        for i in range(natoms)
    ])
    if is_direct:
        coords = coords @ lattice  # Convert from Direct to Cartesian
    return coords

def minimum_image(vec, lattice, inv_lattice):
    # Apply minimum image convention to a displacement vector
    frac = vec @ inv_lattice           # Convert to fractional
    frac -= np.round(frac)            # Apply periodic boundary condition
    return frac @ lattice             # Convert back to Cartesian

# Main loop to compute distances
list_sum = []

for pair in name:
    idx_a = pair[0] - 1
    idx_b = pair[1] - 1
    dist_list = []

    for i in range(1, num + 1):
        coords = get_cart_coords(f'./pfiles/p{i:04d}')
        delta = coords[idx_a] - coords[idx_b]
        delta = minimum_image(delta, lattice, inv_lattice)
        dist = np.linalg.norm(delta)
        dist_list.append(dist)

    list_sum.append(dist_list)

# Save results to file
np.savetxt('bond_length_time.dat', np.array(list_sum).T, header=str(name))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

import os
import numpy as np
import warnings
warnings.filterwarnings('ignore')

# Set the path to the current script directory
os.chdir(os.path.dirname(os.path.abspath(__file__)))
print('copyright by misaraty (misaraty@163.com)\n' + 'last update: 2025-06-26\n')

# Parameter settings
num = 1000  # Number of frames
name = [[3, 2, 5], [4, 2, 5]]  # Angles defined by triplets A-B-C, B is the vertex

# Read lattice info and coordinate type from p0001
with open('./pfiles/p0001') as f:
    lines = f.readlines()

# Lattice matrix and its inverse
lattice = np.array([
    list(map(float, lines[2].split())),
    list(map(float, lines[3].split())),
    list(map(float, lines[4].split()))
])
inv_lattice = np.linalg.inv(lattice)

# Check if coordinates are in Direct mode
is_direct = lines[7].strip().startswith('Direct')

# Total number of atoms
natoms = sum(map(int, lines[6].split()))
coord_start = 8  # Line where atomic coordinates start

# Function to read coordinates and convert to Cartesian if needed
def get_cart_coords(filepath):
    with open(filepath) as f:
        lines = f.readlines()
    coords = np.array([
        list(map(float, lines[coord_start + i].split()))
        for i in range(natoms)
    ])
    if is_direct:
        coords = coords @ lattice
    return coords

# Apply the minimum image convention
def minimum_image(vec, lattice, inv_lattice):
    frac = vec @ inv_lattice
    frac -= np.round(frac)
    return frac @ lattice

# Main loop to compute angles over all frames
list_sum = []

for triplet in name:
    idx_a, idx_b, idx_c = [x - 1 for x in triplet]
    angle_list = []

    for i in range(1, num + 1):
        coords = get_cart_coords(f'./pfiles/p{i:04d}')
        posA = coords[idx_a]
        posB = coords[idx_b]
        posC = coords[idx_c]

        vecBA = minimum_image(posA - posB, lattice, inv_lattice)
        vecBC = minimum_image(posC - posB, lattice, inv_lattice)

        # Calculate angle in degrees
        dot_product = np.dot(vecBA, vecBC)
        norm_BA = np.linalg.norm(vecBA)
        norm_BC = np.linalg.norm(vecBC)
        cos_theta = np.clip(dot_product / (norm_BA * norm_BC), -1.0, 1.0)
        angle_deg = np.degrees(np.arccos(cos_theta))

        angle_list.append(angle_deg)

    list_sum.append(angle_list)

# Save angle vs. time data
np.savetxt('angle_time.dat', np.array(list_sum).T, header=str(name))