拷打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

自动进行差分电荷分析

生成差分电荷密度图和数据文件