这里给出一个CP2K计算shell脚本:
主要功能是建立一个计算文件夹,在其下实现QS下GPW方法的晶格优化,脚本中可选择结构优化,或者MD,以及对角化或者OT方法,示例材料是六角结构单层InP3的计算,可是CP2K并不能得到和平面波方法一致的晶格结构,所以仅供脚本参考:
#!/bin/bash
#1文件夹名称,并创建此文件夹,同时文件夹名也用于后期计算文本文件输入和输出等标识
file_pre="cp2k_InP3"
file_comp="ubuntu"
file_date="210312"
workdir="${file_pre}_${file_comp}_${file_date}_cellOpt_diagNosmear"
if test ! -d ${workdir}
then
echo "${workdir} not exist and make!!!"
mkdir ${workdir}
else
echo "${workdir} already exist and not make!!!"
fi
cd ${workdir}
#2创建一个文本文件,用于计算结束后提取计算输出文件中的能量等信息
etot_file="${workdir}_etot"
if test -r "${etot_file}"
then
rm "${etot_file}"
fi
cat >> ${etot_file} << EOF2
#../${workdir}/${etot_file}
EOF2
grid_header=true
#3CP2K的运行环境变量,根据自己安装情况选择
ubuntu_set="/安装路径/cp2k/tools/toolchain/install/setup"
cp2k_set_name=${ubuntu_set}
source ${cp2k_set_name}
#openmpi等并行运算的前缀
pre_cp2k="mpirun -np 16 "
#cp2k运行名称:
cp2k_name="cp2k.popt"
#是否是计算中断后restart运算,是则t,否则f
res_j="f"
#4 cp2k的QS的GPW方法重要的收敛参数设置
cutoff_sh=350
#用于参数收敛测试循环
for rel_cutoff_sh in 60
do
file_count="cutoff${cutoff_sh}_relcut${rel_cutoff_sh}"
pre_fix="${file_pre}_${file_count}"
cp2k_filename_in="${file_pre}_${file_comp}_${file_date}_${file_count}.inp"
cp2k_filename_out="${file_pre}_${file_comp}_${file_date}_${file_count}.oup"
#${pre_cp2k} ${cp2k_name} -i ./${restart_cp2k_filename_in} -o ./${restart_cp2k_filename_out}
restart_cp2k_filename_in="${pre_fix}-1.restart"
restart_cp2k_filename_out="${pre_fix}_restart.out"
cat > ${cp2k_filename_in} << EOF0
&GLOBAL
PROJECT ${pre_fix}
RUN_TYPE CELL_OPT
#ENERGY ENERGY_FORCE CELL_OPT GEO_OPT MD
PRINT_LEVEL LOW
!#MEDIUM LOW
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&SUBSYS
&KIND Sn
ELEMENT Sn
BASIS_SET DZVP-MOLOPT-SR-GTH
!#DZVP-MOLOPT-SR-GTH SZV-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
!# GTH-PBE GTH-PADE
&END KIND
&KIND Ge
ELEMENT Ge
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
!# GTH-PBE GTH-PADE-q4
&END KIND
&KIND Ga
ELEMENT Ga
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
!# (GTH-PBE-q13 GTH-PBE) GTH-PBE-q3 GTH-PADE-q3 (GTH-PADE-q13 GTH-LDA-q13 GTH-PADE GTH-LDA)
&END KIND
&KIND In
ELEMENT In
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE
!# (GTH-PBE-q13 GTH-PBE) GTH-PBE-q3 GTH-PADE-q3 (GTH-PADE-q13 GTH-LDA-q13 GTH-PADE GTH-LDA)
&END KIND
&KIND P
ELEMENT P
BASIS_SET DZVP-MOLOPT-SR-GTH
!#BASIS_MOLOPT:DZVP-MOLOPT-SR-GTH TZVP-MOLOPT-GTH TZV2P-MOLOPT-GTH
!#BASIS_pob:pob-TZVP plus-pob-TZVP
!#BASIS_MOLOPT_UCL:
POTENTIAL GTH-PBE
!#GTH-PADE GTH-PBE
&END KIND
&KIND O
ELEMENT O
BASIS_SET DZVP-MOLOPT-GTH
!#BASIS_MOLOPT:DZVP-MOLOPT-GTH TZVP-MOLOPT-GTH TZV2P-MOLOPT-GTH TZV2PX-MOLOPT-GTH
!#BASIS_pob:
!#BASIS_MOLOPT_UCL:
POTENTIAL GTH-PBE
!#GTH-PADE GTH-PBE
&END KIND
&CELL
!#InP3_HEX60_QEdeg0.02relax
ABC [angstrom] 7.539600 7.539600 20
ALPHA_BETA_GAMMA 90 90 60
PERIODIC XYZ
SYMMETRY HEXAGONAL
!# CUBIC HEXAGONAL ORTHORHOMBIC TETRAGONAL_AB
MULTIPLE_UNIT_CELL 1 1 1
&END CELL
&COORD
P 0.644488351 0.211022298 0.093732097
P 0.644397199 0.644449544 0.093704383
P 0.211153257 0.644449544 0.093704383
In 0.833294676 0.833412648 0.126252974
In 0.166686801 0.166628397 0.120714080
P 0.355577966 0.355515472 0.153470304
P 0.788905562 0.355515472 0.153470304
P 0.355496188 0.789007625 0.153450475
#UNIT ANGSTROM
SCALED
&END COORD
&PRINT
&ATOMIC_COORDINATES
&END ATOMIC_COORDINATES
&END PRINT
&END SUBSYS
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT_UCL
BASIS_SET_FILE_NAME BASIS_MOLOPT
BASIS_SET_FILE_NAME BASIS_MOLOPT_LnPP1
#! BASIS_MOLOPT BASIS_def2_QZVP_RI_ALL BASIS_pob
POTENTIAL_FILE_NAME POTENTIAL
&QS
METHOD GPW
EPS_DEFAULT 1.0E-10
!#Default value: 1.00000000E-010
EXTRAPOLATION ASPC
!# ps ASPC
#EXTRAPOLATION_ORDER 3
&END QS
&MGRID
NGRIDS 5
CUTOFF ${cutoff_sh}
REL_CUTOFF ${rel_cutoff_sh}
&END MGRID
&XC
&XC_FUNCTIONAL PBE
!# PBE PADE
&END XC_FUNCTIONAL
&END XC
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-6
!#Default value: 1.00000000E-005
MAX_SCF 200
&DIAGONALIZATION .TRUE.
ALGORITHM STANDARD
&END DIAGONALIZATION
&OT .FALSE.
PRECONDITIONER FULL_ALL
!#FULL_KINETIC(use for very large systems-DEFAULT)
!#FULL_ALL (Most effective state selective preconditioner based on diagonalization, requires ENERGY_GAP to be an underestimate of the HOMO-LUMO gap. This preconditioner is recommended for almost all systems, except very large systems)
ENERGY_GAP 0.002
MINIMIZER DIIS
!#CG(DEFAULTmost reliable, use for difficult systems) DIIS(less reliable than CG, but sometimes about 50% faster) SD(not recommended) BROYDEN
&END OT
&MIXING T
METHOD BROYDEN_MIXING
!#BROYDEN_MIXING BROYDEN_MIXING_NEW DIRECT_P_MIXING PULAY_MIXING
ALPHA 0.4
NBROYDEN 8
&END MIXING
# ADDED_MOS 40
&SMEAR .FALSE.
METHOD FERMI_DIRAC
!# ENERGY_WINDOW:WINDOW_SIZE FERMI_DIRAC:ELECTRONIC_TEMPERATURE LIST
ELECTRONIC_TEMPERATURE [K] 300
&END SMEAR
&PRINT
&RESTART OFF
&END
&END PRINT
&END SCF
&POISSON
PERIODIC XYZ
#POISSON_SOLVER PERIODIC
&END POISSON
#PLUS_U_METHOD MULLIKEN
!# MULLIKEN LOWDIN MULLIKEN_CHARGES
&END DFT
#!force_eval
&PRINT
&FORCES HIGH
&END FORCES
&STRESS_TENSOR
&END STRESS_TENSOR
&END PRINT
STRESS_TENSOR ANALYTICAL
!#ANALYTICAL DIAGONAL_ANALYTICAL DIAGONAL_NUMERICAL NUMERICAL
&END FORCE_EVAL
&MOTION
&MD
ENSEMBLE NVT
!# NVE NVT
TEMPERATURE [K] 300
TIMESTEP [fs] 1
STEPS 6000
&THERMOSTAT
REGION GLOBAL
#TYPE NOSE
&END THERMOSTAT
&END MD
&GEO_OPT
TYPE MINIMIZATION
MAX_DR 3.0E-03
MAX_FORCE 4.5E-04
RMS_DR 1.5E-03
RMS_FORCE 3.0E-04
#MAX_ITER 200
OPTIMIZER BFGS
!#CG
!#Good choice for 'small' systems (use LBFGS for large systems)
&END GEO_OPT
&CELL_OPT
TYPE DIRECT_CELL_OPT
#!DIRECT_CELL_OPT(default) GEO_OPT MD
CONSTRAINT Z
!#YZ Z
EXTERNAL_PRESSURE 0 0 0 0 0 0 0 0 0
KEEP_SYMMETRY T
KEEP_ANGLES T
MAX_DR 3.0E-003
MAX_FORCE 4.5E-004
#RMS_DR 1.5E-03
#RMS_FORCE 3.0E-04
#MAX_ITER 200
OPTIMIZER BFGS
!#CG LBFGS BFGS(default)
PRESSURE_TOLERANCE 100
!#100 bar
&END CELL_OPT
#&CONSTRAINT
# &FIXED_ATOMS
# COMPONENTS_TO_FIX XYZ
# LIST 1
# &END FIXED_ATOMS
#&END CONSTRAINT
!# MOTION PRINT
&PRINT
&TRAJECTORY
&EACH
MD 2
!# 1
&END EACH
&END TRAJECTORY
&VELOCITIES
&EACH
MD 2
!# 1
&END EACH
&END VELOCITIES
&FORCES
&EACH
MD 2
!# 1
&END EACH
&END FORCES
&RESTART_HISTORY
&EACH
MD 500
!# 500
CELL_OPT 50
!# 1
GEO_OPT 50
!# 1
&END EACH
&END RESTART_HISTORY
&RESTART
BACKUP_COPIES 3
&EACH
MD 200
!# 20
CELL_OPT 20
!# 1
GEO_OPT 20
!# 1
&END EACH
&END RESTART
&END PRINT
&END MOTION
EOF0
if [ "${res_j}" = "t" ]
then
echo "restart the cp2k (${pre_cp2k} ${cp2k_name} -i ./${restart_cp2k_filename_in} -o ./${restart_cp2k_filename_out} ) calculation..."
${pre_cp2k} ${cp2k_name} -i ./${restart_cp2k_filename_in} -o ./${restart_cp2k_filename_out}
echo "done!!!"
else
echo "running the cp2k ( as "${pre_cp2k} ${cp2k_name} -i ./${cp2k_filename_in} -o ./${cp2k_filename_out}" )"
${pre_cp2k} ${cp2k_name} -i ./${cp2k_filename_in} -o ./${cp2k_filename_out}
echo "done!!!"
fi
total_energy=$(grep -e '^[ \t]*Total energy' ${cp2k_filename_out} | awk '{print $3}')
ngrids=$(grep -e '^[ \t]*QS| Number of grid levels:' ${cp2k_filename_out} | awk '{print $6}')
if ${grid_header} ; then
printf "# Cutoff (Ry) | Total Energy (Ha)" >> ${etot_file}
for ((igrid=1; igrid <= ngrids; igrid++)) ; do
printf "| NG on grid %d" $igrid >> ${etot_file}
done
printf "\n" >> ${etot_file}
grid_header=false
fi
printf "%10.2f %10.2f %15.10f" ${cutoff_sh} ${rel_cutoff_sh} $total_energy >> ${etot_file}
for ((igrid=1; igrid <= ngrids; igrid++)) ; do
grid=$(grep -e '^[ \t]*count for grid' ${cp2k_filename_out} | \
awk -v igrid=$igrid '(NR == igrid){print $5}')
printf " %6d" $grid >> ${etot_file}
done
printf "\n" >> ${etot_file}
done
更多推荐
CP2K计算脚本
发布评论