拷打DeepSeek實現自動生成差分電荷計算文件和后處理

差分電荷（charge density difference）是材料模擬中分析電子結構變化的直觀工具。
它把成鍵后的真實電荷密度減去成鍵前各碎片疊加的電荷密度，得到一張“電子遷移地圖”

于是可以一眼看出化學鍵形成時電子從哪里來到哪里去，表面吸附、缺陷、界面、催化活性位點等問題的本質都能用這張圖說話。

AI Promt

差分電荷的計算邏輯是把POSCAR拆分為好幾個部分再進行自洽計算，然后再通過vaspkit進行數據處理和輸出，

這樣，幫我寫一個腳本，可以識別POSCAR中的各類元素以及每個元素與多少原子，并把他們的信息輸出出來，提示用戶去選擇，可以選擇拆分某一個原子，或者某一類元素的所有原子或部分，或者某一層（z坐標在某個范圍內的）的原子，也可以自定義幾個原子的基團等等，總之這個地方需要足夠好的交互界面，

然后，依據用戶選定的結果，將POSCAR拆分為選定原子的結構和排除掉選定原子之后剩余原子的結構，拆分后的結構的晶格常數與拆之前保持完全一致，所有結構原子的坐標也保持不變，然后將未拆分之前的整體的結構，選定原子的結構，未選定原子的結構分別存放在三個文件夾中?

?POSCAR的模版如下

?POSCAR ? ? ? ? ? ??

1.00000000000000 ? ? ? ? ?4.8476366827364270 ? ?0.0000000000000000 ? ?0.0000000000000000 ? ?-2.4238183414183960 ? ?4.1981765155659660 ? ?0.0000000000000000 ? ? 0.0000000000000000 ? ?0.0000000000000000 ? 25.0000000000000000 ? Mn ? Pb ? Na ? ? 1 ? ? 1 ? ? 2Direct -0.0000000000000000 -0.0000000000000000 ?0.5000000000000000 ?0.6666666870000029 ?0.3333333429999996 ?0.5000000000000000 ?0.3333333429999996 ?0.6666666870000029 ?0.4252117344687524 ?0.3333333429999996 ?0.6666666870000029 ?0.5747882955312500

腳本中代碼，注釋只輸出英文字符

增加功能：生成計算文件后，進入到未拆分的結構所在的文件夾，執行vaspkit的命令（echo 102；echo 2；echo 0.04）|vaspkit生成精度為0.04（讓執行者確認并可調整）的KPOINTS文件，同時使用（echo 101；echo ST）|vaspkit生成靜態自洽計算的INCAR，同時，詢問用戶是否需要在INCAR中+U或者+SOC，對于需要+U，則執行（echo 101；echo PUST）|vaspkit對于需要+SOC，則執行（echo 101；echo SOCST）|vaspkit對于都加，則（echo 101；echo PUSOCST）|vaspkit執行完畢后請用戶核對INCAR的內容并確認再進行下一步然后，檢索INCAR中關于 ENCUT的關鍵詞，如果存在則檢查是否被注釋掉（前面有#），如果被注釋掉則刪去前面的符合#或者直接在INCAR最開頭再加一行ENCUT的參數ENCUT?= 500然后將KPOINTS和INCAR復制到選定原子結構文件夾中和未選定原子文件夾中，然后，分別在這三個文件夾中執行（echo 103）|vaspkit用于生成POTCAR文件最后，檢查核對，三個文件夾中，POSCAR里面所含有的原子元素，尤其是順序，和POTCAR中所含有的元素對應，對于PBE泛函，POT的開頭格式為 ?PAW_PBE Mn_pv 02Aug2007 ? ? ? ? ? ? ? ? ? 13.0000000000000 ? ? ?parameters from PSCTR are: ? VRHFIN =Mn: 3p4s3d ? LEXCH ?= PE ? EATOM ?= ?2025.2664 eV, ?148.8529 Ry ? TITEL ?= PAW_PBE Mn_pv 02Aug2007 ? LULTRA = ? ? ? ?F ? ?use ultrasoft PP ?可通過查找 PBE 或者TITEL 或者元素名詞關鍵詞來找到元素標記，最后再檢查所有文件夾的INCAR中,ENCUT以及SIGMA等參數是否一致,如果有不一致的的參數需要提醒用戶

新增功能:在完成vasp計算后，調用vaspkit進行數據后處理，即使用vaspkit的314功能，命令為(echo 314 ; echo 未拆分結構文件夾/CHGCAR 選定原子結構文件夾/CHGCAR 未選定原子結構文件夾/CHGCAR)|vaspkit這一步會生成文件CHGDIFF.vasp ，如果未生成，則提示用戶檢查每一個文件夾中的計算過程和輸出隨后將CHGDIFF.vasp 拷貝為CHGCAR在當然目錄，并通過vaspkit繪制面平均差分電荷密度(echo 316 ;echo 1;echo 3 ?)|vaspkit

代碼

#!/usr/bin/env python3
import os
import numpy as np
import shutil
import subprocess
from pathlib import Pathclass POSCARProcessor:def __init__(self, poscar_file='POSCAR'):self.poscar_file = poscar_fileself.lattice_constant = 1.0self.lattice_vectors = []self.element_names = []self.element_counts = []self.coordinate_type = 'Direct'self.atom_coordinates = []self.total_atoms = 0self.element_indices = {}  # Mapping from element to atom indicesdef read_poscar(self):"""Read POSCAR file"""with open(self.poscar_file, 'r') as f:lines = f.readlines()# Read titleself.title = lines[0].strip()# Read lattice constantself.lattice_constant = float(lines[1].strip())# Read lattice vectorsself.lattice_vectors = []for i in range(2, 5):vec = list(map(float, lines[i].split()))self.lattice_vectors.append(vec)# Read element namesself.element_names = lines[5].split()# Read element countsself.element_counts = list(map(int, lines[6].split()))self.total_atoms = sum(self.element_counts)# Read coordinate typeself.coordinate_type = lines[7].strip()# Read atom coordinatesself.atom_coordinates = []for i in range(8, 8 + self.total_atoms):coords = list(map(float, lines[i].split()[:3]))self.atom_coordinates.append(coords)# Build element to atom index mappingself._build_element_indices()def _build_element_indices(self):"""Build mapping from element to atom indices"""self.element_indices = {}atom_index = 0for element, count in zip(self.element_names, self.element_counts):self.element_indices[element] = list(range(atom_index, atom_index + count))atom_index += countdef display_atom_info(self):"""Display atom information"""print("\n" + "="*50)print("POSCAR File Information")print("="*50)print(f"Title: {self.title}")print(f"Lattice constant: {self.lattice_constant}")print(f"Elements: {self.element_names}")print(f"Atom counts per element: {self.element_counts}")print(f"Total atoms: {self.total_atoms}")print(f"Coordinate type: {self.coordinate_type}")print("\nAtom index ranges for each element:")for element, indices in self.element_indices.items():print(f"  {element}: {indices[0] + 1}-{indices[-1] + 1}")print("\nAtom coordinates (first 5 atoms):")for i in range(min(5, self.total_atoms)):print(f"  Atom {i+1}: {self.atom_coordinates[i]}")if self.total_atoms > 5:print("  ...")def get_user_selection(self):"""Get user selection"""print("\n" + "="*50)print("Select atoms to split")print("="*50)print("1. Select single atom")print("2. Select all atoms of an element")print("3. Select partial atoms of an element")print("4. Select atoms in z-coordinate range")print("5. Custom atom selection")print("6. Finish selection")selected_indices = set()while True:try:choice = input("\nSelect operation (1-6): ").strip()if choice == '6':breakif choice == '1':atom_idx = int(input("Enter atom number (1-based): ")) - 1if 0 <= atom_idx < self.total_atoms:selected_indices.add(atom_idx)print(f"Selected atom {atom_idx + 1}")else:print("Invalid atom number")elif choice == '2':element = input("Enter element symbol: ").strip()if element in self.element_indices:selected_indices.update(self.element_indices[element])print(f"Selected all {element} atoms")else:print("Element not found")elif choice == '3':element = input("Enter element symbol: ").strip()if element in self.element_indices:indices = self.element_indices[element]print(f"Available {element} atom numbers: {[i+1 for i in indices]}")selected = input("Enter atom numbers to select (space separated): ").split()for idx in selected:try:atom_idx = int(idx) - 1if atom_idx in indices:selected_indices.add(atom_idx)except ValueError:passprint(f"Selected specified {element} atoms")else:print("Element not found")elif choice == '4':try:z_min = float(input("Enter minimum z-coordinate: "))z_max = float(input("Enter maximum z-coordinate: "))for i, coords in enumerate(self.atom_coordinates):if z_min <= coords[2] <= z_max:selected_indices.add(i)print(f"Selected atoms with z-coordinate in [{z_min}, {z_max}]")except ValueError:print("Invalid input")elif choice == '5':try:custom_indices = input("Enter atom numbers to select (space separated): ").split()for idx in custom_indices:atom_idx = int(idx) - 1if 0 <= atom_idx < self.total_atoms:selected_indices.add(atom_idx)print("Selected specified atoms")except ValueError:print("Invalid input")else:print("Invalid selection")except KeyboardInterrupt:print("\nOperation cancelled")return Noneexcept Exception as e:print(f"Error occurred: {e}")return sorted(selected_indices)def write_poscar(self, filename, selected_indices, title_suffix):"""Write POSCAR file"""# Calculate element counts in new structureelement_counts_new = [0] * len(self.element_names)for idx in selected_indices:# Determine which element this atom belongs tocumulative = 0for i, count in enumerate(self.element_counts):if idx < cumulative + count:element_counts_new[i] += 1breakcumulative += count# Remove elements with zero countelement_names_new = []element_counts_final = []for i, count in enumerate(element_counts_new):if count > 0:element_names_new.append(self.element_names[i])element_counts_final.append(count)with open(filename, 'w') as f:f.write(f"{self.title} {title_suffix}\n")f.write(f"  {self.lattice_constant:.16f}\n")for vec in self.lattice_vectors:f.write("  " + " ".join(f"{x:.16f}" for x in vec) + "\n")f.write("  " + " ".join(element_names_new) + "\n")f.write("  " + " ".join(map(str, element_counts_final)) + "\n")f.write(f"{self.coordinate_type}\n")for idx in selected_indices:coords = self.atom_coordinates[idx]f.write("  " + " ".join(f"{x:.16f}" for x in coords) + "\n")def run_vaspkit_command(self, commands, folder='.'):"""Run vaspkit command with given input"""try:process = subprocess.Popen(['vaspkit'],stdin=subprocess.PIPE,stdout=subprocess.PIPE,stderr=subprocess.PIPE,text=True,cwd=folder)stdout, stderr = process.communicate(commands)return stdout, stderr, process.returncodeexcept Exception as e:print(f"Error running vaspkit: {e}")return None, None, -1def generate_kpoints(self, folder='.'):"""Generate KPOINTS file with user-specified precision"""kspacing = input("Enter KPOINTS precision (default 0.04): ").strip()if not kspacing:kspacing = "0.04"commands = f"102\n2\n{kspacing}\n0\n"stdout, stderr, returncode = self.run_vaspkit_command(commands, folder)if returncode == 0:print(f"KPOINTS generated successfully in {folder}")return Trueelse:print(f"Failed to generate KPOINTS in {folder}")return Falsedef generate_incar(self, folder='.'):"""Generate INCAR file with optional +U and +SOC"""print("\nSelect INCAR options:")print("1. Standard static calculation (ST)")print("2. +U calculation (PUST)")print("3. +SOC calculation (SOCST)")print("4. +U and +SOC calculation (PUSOCST)")choice = input("Enter choice (1-4, default 1): ").strip()if not choice:choice = "1"option_map = {"1": "ST","2": "PUST","3": "SOCST","4": "PUSOCST"}option = option_map.get(choice, "ST")commands = f"101\n{option}\n"stdout, stderr, returncode = self.run_vaspkit_command(commands, folder)if returncode == 0:print(f"INCAR generated successfully in {folder}")# Set ENCUT to 500self.set_encut(folder)# Show INCAR content for verificationself.show_incar_content(folder)confirm = input("Confirm INCAR content? (y/n, default y): ").strip().lower()if confirm in ['', 'y', 'yes']:return Trueelse:print("Please modify INCAR manually and continue")return Trueelse:print(f"Failed to generate INCAR in {folder}")return Falsedef set_encut(self, folder='.'):"""Set ENCUT to 500 in INCAR"""incar_file = os.path.join(folder, 'INCAR')if not os.path.exists(incar_file):return Falsewith open(incar_file, 'r') as f:lines = f.readlines()encut_found = Falsenew_lines = []for line in lines:if line.strip().startswith('ENCUT') or line.strip().startswith('#ENCUT'):# Remove comment and set ENCUTclean_line = line.replace('#', '').replace('!', '').strip()if '=' in clean_line:new_line = 'ENCUT = 500\n'else:new_line = 'ENCUT = 500\n'new_lines.append(new_line)encut_found = Trueelse:new_lines.append(line)# If ENCUT not found, add it at the beginningif not encut_found:new_lines.insert(0, 'ENCUT = 500\n')with open(incar_file, 'w') as f:f.writelines(new_lines)print(f"ENCUT set to 500 in {folder}/INCAR")return Truedef show_incar_content(self, folder='.'):"""Show INCAR content for verification"""incar_file = os.path.join(folder, 'INCAR')if os.path.exists(incar_file):print(f"\nINCAR content in {folder}:")print("-" * 40)with open(incar_file, 'r') as f:content = f.read()print(content)print("-" * 40)def generate_potcar(self, folder='.'):"""Generate POTCAR file"""commands = "103\n"stdout, stderr, returncode = self.run_vaspkit_command(commands, folder)if returncode == 0:print(f"POTCAR generated successfully in {folder}")return Trueelse:print(f"Failed to generate POTCAR in {folder}")return Falsedef check_element_consistency(self, folder='.'):"""Check consistency between POSCAR and POTCAR elements"""poscar_file = os.path.join(folder, 'POSCAR')potcar_file = os.path.join(folder, 'POTCAR')if not os.path.exists(poscar_file) or not os.path.exists(potcar_file):print(f"Missing files in {folder}")return False# Read elements from POSCARwith open(poscar_file, 'r') as f:poscar_lines = f.readlines()poscar_elements = poscar_lines[5].split()# Read elements from POTCARwith open(potcar_file, 'r', errors='ignore') as f:potcar_content = f.read()potcar_elements = []lines = potcar_content.split('\n')for i, line in enumerate(lines):if 'TITEL' in line or 'PAW_PBE' in line:# Extract element from line like: "PAW_PBE Mn_pv 02Aug2007"parts = line.split()if len(parts) >= 2:element_part = parts[1]element = element_part.split('_')[0]  # Extract Mn from Mn_pvpotcar_elements.append(element)# Check consistencyif poscar_elements == potcar_elements:print(f"? Element consistency check passed in {folder}")print(f"  POSCAR elements: {poscar_elements}")print(f"  POTCAR elements: {potcar_elements}")return Trueelse:print(f"? Element inconsistency in {folder}")print(f"  POSCAR elements: {poscar_elements}")print(f"  POTCAR elements: {potcar_elements}")return Falsedef check_incar_consistency(self):"""Check consistency of INCAR parameters across folders"""folders = ['original', 'selected', 'remaining']incar_params = {}for folder in folders:incar_file = os.path.join(folder, 'INCAR')if os.path.exists(incar_file):with open(incar_file, 'r') as f:lines = f.readlines()params = {}for line in lines:if '=' in line and not line.strip().startswith('#'):key, value = line.split('=', 1)key = key.strip()value = value.split('!')[0].split('#')[0].strip()params[key] = valueincar_params[folder] = params# Compare parametersall_keys = set()for params in incar_params.values():all_keys.update(params.keys())inconsistencies = []for key in all_keys:values = {}for folder, params in incar_params.items():values[folder] = params.get(key, 'NOT_SET')if len(set(values.values())) > 1:inconsistencies.append((key, values))if inconsistencies:print("\n??  INCAR parameter inconsistencies found:")for key, values in inconsistencies:print(f"  {key}:")for folder, value in values.items():print(f"    {folder}: {value}")else:print("\n? All INCAR parameters are consistent across folders")return not inconsistenciesdef check_vasp_calculation_completion(self, folder='.'):"""Check if VASP calculation completed successfully"""outcar_file = os.path.join(folder, 'OUTCAR')oszicar_file = os.path.join(folder, 'OSZICAR')if not os.path.exists(outcar_file) or not os.path.exists(oszicar_file):return False# Check if OUTCAR contains completion messagewith open(outcar_file, 'r', errors='ignore') as f:outcar_content = f.read()# Check if OSZICAR contains reasonable number of stepswith open(oszicar_file, 'r') as f:oszicar_lines = f.readlines()# Check for completion indicatorscompletion_indicators = ['Voluntary context switches','General timing and accounting','Total CPU time used']has_completion = any(indicator in outcar_content for indicator in completion_indicators)has_reasonable_steps = len(oszicar_lines) > 5  # At least some SCF stepsreturn has_completion and has_reasonable_stepsdef check_chgcar_exists(self, folder='.'):"""Check if CHGCAR file exists"""chgcar_file = os.path.join(folder, 'CHGCAR')return os.path.exists(chgcar_file)def perform_charge_difference_analysis(self):"""Perform charge difference analysis using vaspkit"""print("\n" + "="*50)print("Performing Charge Difference Analysis")print("="*50)# Check if VASP calculations completed successfullyfolders = ['original', 'selected', 'remaining']all_calculations_completed = Truemissing_chgcar = []for folder in folders:if not self.check_vasp_calculation_completion(folder):print(f"??  VASP calculation in {folder} may not have completed successfully")all_calculations_completed = Falseif not self.check_chgcar_exists(folder):print(f"??  CHGCAR not found in {folder}")missing_chgcar.append(folder)if not all_calculations_completed or missing_chgcar:print("\n? Cannot perform charge difference analysis:")if not all_calculations_completed:print("  - Some VASP calculations did not complete successfully")if missing_chgcar:print(f"  - Missing CHGCAR in: {missing_chgcar}")print("Please check the calculation outputs and run VASP calculations first.")return Falseprint("All VASP calculations completed successfully and CHGCAR files are available")# Generate charge difference using vaspkitprint("\nGenerating charge difference (CHGDIFF.vasp)...")commands = "314\noriginal/CHGCAR\nselected/CHGCAR\nremaining/CHGCAR\n"stdout, stderr, returncode = self.run_vaspkit_command(commands)if returncode == 0 and os.path.exists('CHGDIFF.vasp'):print("? CHGDIFF.vasp generated successfully")# Copy CHGDIFF.vasp to CHGCARshutil.copy2('CHGDIFF.vasp', 'CHGCAR')print("? Copied CHGDIFF.vasp to CHGCAR")# Plot plane-averaged charge differenceprint("\nPlotting plane-averaged charge difference...")plot_commands = "316\n1\n3\n"plot_stdout, plot_stderr, plot_returncode = self.run_vaspkit_command(plot_commands)if plot_returncode == 0:print("? Plane-averaged charge difference plot generated successfully")# Check if plot files were createdplot_files = ['PLANE-AVERAGE.dat', 'PLANE-AVERAGE.png']for plot_file in plot_files:if os.path.exists(plot_file):print(f"  Generated: {plot_file}")return Trueelse:print("? Failed to generate plane-averaged charge difference plot")return Falseelse:print("? Failed to generate CHGDIFF.vasp")if stdout:print(f"vaspkit stdout: {stdout}")if stderr:print(f"vaspkit stderr: {stderr}")return Falsedef process_selection(self, selected_indices):"""Process selection and create three folders"""if not selected_indices:print("No atoms selected")return# Create foldersfolders = ['original', 'selected', 'remaining']for folder in folders:if os.path.exists(folder):shutil.rmtree(folder)os.makedirs(folder)# Save original POSCARshutil.copy2(self.poscar_file, 'original/POSCAR')# Save selected atomsself.write_poscar('selected/POSCAR', selected_indices, 'Selected')# Save remaining atomsremaining_indices = [i for i in range(self.total_atoms) if i not in selected_indices]self.write_poscar('remaining/POSCAR', remaining_indices, 'Remaining')print(f"\nStructure splitting completed!")print(f"Original structure saved in: original/")print(f"Selected atoms structure saved in: selected/")print(f"Remaining atoms structure saved in: remaining/")print(f"Number of selected atoms: {len(selected_indices)}")print(f"Number of remaining atoms: {len(remaining_indices)}")# Generate calculation filesprint("\n" + "="*50)print("Generating calculation files")print("="*50)# Generate files in original folder firstprint("\nGenerating files in original folder...")self.generate_kpoints('original')self.generate_incar('original')# Copy files to other foldersfor folder in ['selected', 'remaining']:print(f"\nCopying files to {folder} folder...")shutil.copy2('original/KPOINTS', f'{folder}/')shutil.copy2('original/INCAR', f'{folder}/')# Generate POTCAR in all foldersprint("\nGenerating POTCAR files...")for folder in folders:print(f"Generating POTCAR in {folder}...")self.generate_potcar(folder)# Check consistencyprint("\n" + "="*50)print("Checking consistency")print("="*50)for folder in folders:self.check_element_consistency(folder)self.check_incar_consistency()print("\n" + "="*50)print("Setup completed!")print("="*50)print("Next steps:")print("1. Run VASP calculations in all three folders:")print("   cd original && mpirun -np 4 vasp_std && cd ..")print("   cd selected && mpirun -np 4 vasp_std && cd ..")print("   cd remaining && mpirun -np 4 vasp_std && cd ..")print("2. After calculations complete, run this script again with --analyze option")print("   python script_name.py --analyze")def analyze_results(self):"""Analyze VASP calculation results"""return self.perform_charge_difference_analysis()def main():import argparseparser = argparse.ArgumentParser(description='POSCAR Processor for Charge Difference Analysis')parser.add_argument('--analyze', action='store_true', help='Analyze VASP calculation results')args = parser.parse_args()# Check if POSCAR file existsif not os.path.exists('POSCAR'):print("Error: POSCAR file not found in current directory")return# Check if vaspkit is availabletry:subprocess.run(['which', 'vaspkit'], check=True, capture_output=True)except subprocess.CalledProcessError:print("Error: vaspkit not found. Please install vaspkit first.")returnprocessor = POSCARProcessor()try:if args.analyze:# Directly analyze resultsprocessor.analyze_results()else:# Normal processing flow# Read POSCAR fileprocessor.read_poscar()# Display atom informationprocessor.display_atom_info()# Get user selectionselected_indices = processor.get_user_selection()if selected_indices is not None:# Process selectionprocessor.process_selection(selected_indices)except Exception as e:print(f"Error during processing: {e}")import tracebacktraceback.print_exc()if __name__ == "__main__":main()

使用流程

第一步：準備階段

python script.py

交互式選擇要拆分的原子

自動生成三個文件夾的計算文件

第二步：運行計算

# 在三個文件夾中分別運行VASP計算

第三步：后處理分析

python script.py --analyze

自動進行差分電荷分析

生成差分電荷密度圖和數據文件