计算化学公社

标题: cp2k怎么批量提出动力学过程中每帧的结构 [打印本页]

作者
Author: zhangs 时间: 2024-10-31 12:47
标题: cp2k怎么批量提出动力学过程中每帧的结构
各位老师，请教下，我用cp2k跑了一个AIMD，现在想计算每一步的单点能，有什么方法可以提出每一帧的轨迹文件吗？

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-10-31 13:13
本帖最后由 Uus/pMeC6H4-/キ于 2024-11-1 11:10 编辑

CP2K的轨迹文件<proj>-pos-1.xyz就是动力学中输出的每一帧的结构写成的.xyz拼成的，那要提取只需很简单的Linux文件裁切命令就行：

split ./"<filename>" "<prefix>" -l <xyz_lines> -d -a <digits> --additional-suffix="<suffix>"

复制代码

比如一个叫md-pos-1.xyz的轨迹文件，结构自始至终有128个原子（这样每一个提取出来的xyz文件都有128+2=130行内容，末尾还有一个额外的空行），打算提取出的文件名以"frame"开头，中间是代表帧号的4位数字，最后以"-cut.xyz"结尾，那切换到轨迹文件所在目录后裁切命令就写作

split ./"md-pos-1.xyz" "frame" -l 130 -d -a 4 --additional-suffix="-cut.xyz"

复制代码

编辑：稍微解释一下，split指令默认用字母给输出文件命名编号排序，故加-d改成用数字。-a 4表示帧号数字一共四位，不够4位就用先导0补齐；具体要几位取决于轨迹文件记录的帧数，按CP2K输入文件&MD/STEPS和&PRINT/&TRAJECTORY/&EACH/MD字段的设置自己算，-a 4适合几千帧的情况，如有上万帧就用-a 5，等等。若-a设大了会有多余先导0，-a设小了则会裁切中途提示"split: output file suffixes exhausted"。另外，由于MD记录的第一帧是编号i=0，上面命令的写法会按默认设置从0000开始命名编号；但结构优化记录的第一帧编号是i=1，上面命令的-d建议改成--numeric-suffixes=1以从0001开始。

作者
Author: sobereva 时间: 2024-10-31 16:24
利用&MD里的&REFTRAJ，直接就可以从xyz轨迹里读取各帧结构并计算能量/受力/属性，不需要特意提出来白走趟弯路

作者
Author: chenzhe 时间: 2024-11-5 09:22
本帖最后由 chenzhe 于 2024-11-5 09:58 编辑

sobereva 发表于 2024-10-31 16:24
利用&MD里的&REFTRAJ，直接就可以从xyz轨迹里读取各帧结构并计算能量/受力/属性，不需要特意提出来白走趟弯 ...

老师，具体&REFTRAJ命令怎么用的？
之前老版本是这样。
&REFTRAJ
   TRAJ_FILE_NAME xxx.xyz
   EVAL_ENERGY_FORCES .TRUE.
   EVAL_FORCES .TRUE.
   FIRST_SNAPSHOT 1
   LAST_SNAPSHOT 50
   STRIDE 1
&END REFTRAJ

但是新版本cp2k手册上写这EVAL_ENERGY_FORCES和EVAL_FORCES已经弃用，不过我没看到新的相应案例。

作者
Author: sobereva 时间: 2024-11-8 21:51

chenzhe 发表于 2024-11-5 09:22
老师，具体&REFTRAJ命令怎么用的？
之前老版本是这样。
&REFTRAJ

在CP2K自带的测试文件里搜例子

作者
Author: exsolution 时间: 2024-11-28 20:47

chenzhe 发表于 2024-11-5 09:22
老师，具体&REFTRAJ命令怎么用的？
之前老版本是这样。
&REFTRAJ

你好请问解决了吗这个cp2k2022版本的字段还和2023不一样吗

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-11-28 20:58

exsolution 发表于 2024-11-28 20:47
你好请问解决了吗这个cp2k2022版本的字段还和2023不一样吗

此字段相关的代码是在2023年9月这个commit中改的，在2024.1版之后才有新格式。

作者
Author: exsolution 时间: 2024-11-28 21:26

Uus/pMeC6H4-/キ发表于 2024-11-28 20:58
此字段相关的代码是在2023年9月这个commit中改的，在2024.1版之后才有新格式。

老师好我在第一遍跑aimd的时候忘记开启输出受力信息了
那么我在已有的轨迹基础上，使用&REFTRAJ 字段重新计算受力信息是不是和重进计算新的aimd 成本差不多呢？

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-11-28 21:53
本帖最后由 Uus/pMeC6H4-/キ于 2024-11-28 22:45 编辑

exsolution 发表于 2024-11-28 21:26
老师好我在第一遍跑aimd的时候忘记开启输出受力信息了
那么我在已有的轨迹基础上，使用&REFTRAJ 字段重 ...

取决于上次分子动力学输出轨迹时记录体系坐标的频率、重算的跨步和计算级别。如果计算级别不变，需要重算的结构数量只有原先轨迹的几分之一，那成本肯定明显减少。

编辑：有余裕的话甚至可以尝试提高计算级别，看这帖开始的讨论。

作者
Author: ljh123 时间: 2024-12-13 21:16

sobereva 发表于 2024-11-8 21:51
在CP2K自带的测试文件里搜例子

sob老师，我是用NPT跑出来的AIMD轨迹，想用REFTRAJ命令重算指定帧的能量，但是发现REFTRAJ重算过程中无法改变盒子尺寸，只能是输入文件中的尺寸。但是NPT每一帧的盒子尺寸都在波动，请问这个有办法解决吗，还是只能提取结构来跑单点能呢

作者
Author: sobereva 时间: 2024-12-14 02:04

ljh123 发表于 2024-12-13 21:16
sob老师，我是用NPT跑出来的AIMD轨迹，想用REFTRAJ命令重算指定帧的能量，但是发现REFTRAJ重算过程中无法 ...

只能提取结构来跑单点能

作者
Author: ljh123 时间: 2024-12-17 01:05

sobereva 发表于 2024-12-14 02:04
只能提取结构来跑单点能

sob老师我在计算单点能时出现了问题。我是在自旋多重度为7时跑的aimd轨迹，并提取结构分别计算了自旋多重度为5和7时的单点能（分别对应电子转移前后，用于计算重组能）。但是当自旋多重度设置为5计算重组能时，大部分结构的 Integrated absolute spin density仍然接近7，导致算出来的重组能几乎为0。请问有什么办法可以让scf收敛到Integrated absolute spin density为5的状态吗

作者
Author: sobereva 时间: 2024-12-17 02:30

ljh123 发表于 2024-12-17 01:05
sob老师我在计算单点能时出现了问题。我是在自旋多重度为7时跑的aimd轨迹，并提取结构分别计算了自旋多重 ...

似乎你开了smearing或者RELAX_MULTIPLICITY？此时自旋多重度没法严格指定，下文提到的ppt里说了
驳网上流传的对CP2K缺点的不实描述
http://sobereva.com/729（http://bbs.keinsci.com/thread-50094-1-1.html）

作者
Author: ljh123 时间: 2024-12-17 02:36

sobereva 发表于 2024-12-17 02:30
似乎你开了smearing或者RELAX_MULTIPLICITY？此时自旋多重度没法严格指定，下文提到的ppt里说了
驳网上 ...

sob老师我没开这两个关键词。后面我又重算了该结构（在自旋多重度为7下跑MD得到结构）自旋多重度设置为5时的单点能。一共算了三次，三次的单点能值都不同，其中两次的 Integrated absolute spin density为4.3左右，一次为6左右

作者
Author: sobereva 时间: 2024-12-17 03:05

ljh123 发表于 2024-12-17 02:36
sob老师我没开这两个关键词。后面我又重算了该结构（在自旋多重度为7下跑MD得到结构）自旋多重度设置为5 ...

没具体文件没法说
没用上述关键词的话，实际得到的波函数的自旋多重度理应和你设的MULTIPLICITY一致

作者
Author: ljh123 时间: 2024-12-17 13:18
本帖最后由 ljh123 于 2024-12-17 13:21 编辑

sobereva 发表于 2024-12-17 03:05
没具体文件没法说
没用上述关键词的话，实际得到的波函数的自旋多重度理应和你设的MULTIPLICITY一致

#Generated by Multiwfn (http://sobereva.com/multiwfn)
&GLOBAL
  PROJECT 3
  PRINT_LEVEL LOW
  RUN_TYPE ENERGY
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep
  &SUBSYS
&CELL
PERIODIC  XYZ
CELL_FILE_FORMAT  CIF
CELL_FILE_NAME  ../CIF/3.cif
&END CELL
&TOPOLOGY
   COORDINATE CIF
   COORD_FILE ../CIF/3.cif
&END TOPOLOGY
&KIND O
   ELEMENT O
   BASIS_SET DZVP-MOLOPT-SR-GTH-q6
   POTENTIAL GTH-PBE
&END KIND
&KIND H
   ELEMENT H
   BASIS_SET DZVP-MOLOPT-SR-GTH-q1
   POTENTIAL GTH-PBE
&END KIND
&KIND Fe
   ELEMENT Fe
   BASIS_SET DZVP-MOLOPT-SR-GTH-q16
   POTENTIAL GTH-PBE
&END KIND
&KIND N
   ELEMENT N
   BASIS_SET DZVP-MOLOPT-SR-GTH-q5
   POTENTIAL GTH-PBE
&END KIND
&KIND C
   ELEMENT C
   BASIS_SET DZVP-MOLOPT-SR-GTH-q4
   POTENTIAL GTH-PBE
&END KIND
&KIND Br
   ELEMENT Br
   BASIS_SET DZVP-MOLOPT-SR-GTH-q7
   POTENTIAL GTH-PBE
&END KIND
  &END SUBSYS

  &DFT
BASIS_SET_FILE_NAME  BASIS_MOLOPT
POTENTIAL_FILE_NAME  POTENTIAL
#  WFN_RESTART_FILE_NAME 2-RESTART.cif
CHARGE 0 #Net charge
MULTIPLICITY 5 #Spin multiplicity
UKS
&QS
   EPS_DEFAULT 1.0E-11 #Set all EPS_xxx to values such that the energy will be correct up to this value
&END QS
&POISSON
   PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
   PSOLVER PERIODIC #The way to solve Poisson equation
&END POISSON
&XC
   &XC_FUNCTIONAL PBE
   &END XC_FUNCTIONAL
   &VDW_POTENTIAL
      POTENTIAL_TYPE PAIR_POTENTIAL
      &PAIR_POTENTIAL
      PARAMETER_FILE_NAME dftd3.dat
      TYPE DFTD3(BJ)
      REFERENCE_FUNCTIONAL PBE
      #CALCULATE_C9_TERM T #Calculate C9-related three-body term, more accurate for large system
      &END PAIR_POTENTIAL
   &END VDW_POTENTIAL
&END XC
&MGRID
   CUTOFF  350
   REL_CUTOFF  50
&END MGRID
&SCF
   MAX_SCF 25 #Maximum number of steps of inner SCF
   EPS_SCF 5.0E-06 #Convergence threshold of density matrix of inner SCF
#    SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
#    IGNORE_CONVERGENCE_FAILURE #Continue calculation even if SCF not converged, works for version >= 2024.1
   &OT
      PRECONDITIONER FULL_KINETIC #FULL_SINGLE_INVERSE is also worth to try. FULL_ALL is better but quite expensive for large system
      MINIMIZER DIIS #CG is worth to consider in difficult cases
      LINESEARCH 2PNT #1D line search algorithm for CG. 2PNT is default. 3PNT is more expensive but may be better. GOLD is best but very expensive
      ALGORITHM STRICT #Algorithm of OT. Can be STRICT (default) or IRAC
   &END OT
   &OUTER_SCF
      MAX_SCF 20 #Maximum number of steps of outer SCF
      EPS_SCF 5.0E-06 #Convergence threshold of outer SCF
   &END OUTER_SCF
   &PRINT
      &RESTART #Note: Use "&RESTART OFF" can prevent generating .wfn file
      BACKUP_COPIES 0 #Maximum number of backup copies of wfn file. 0 means never
      &END RESTART
   &END PRINT
&END SCF
  &END DFT
&END FORCE_EVALnull

sob老师这是我的输入文件。结构是在自旋多重度为7的设定下的跑MD得到的。但是这个单点能计算完后，输出文件里的总自旋多重度变成了7

作者
Author: sobereva 时间: 2024-12-18 02:54

ljh123 发表于 2024-12-17 13:18
#Generated by Multiwfn (http://sobereva.com/multiwfn)
&GLOBAL
PROJECT 3

Integrated absolute spin density对应的是下图的M_abs。仅当体系中的自旋密度都>=0，它才等于M_tot，而M_tot才与自旋多重度相对应。所以你看到两次Integrated absolute spin density明显不同时，应该看自旋密度图，4.3的那个应当自旋密度有正有负。两次得到不同Integrated absolute spin density结果说明偶然性地收敛到了不同波函数。

(, 下载次数 Times of downloads: 6)

作者
Author: JCenter 时间: 2025-5-24 03:33

sobereva 发表于 2024-10-31 16:24
利用&MD里的&REFTRAJ，直接就可以从xyz轨迹里读取各帧结构并计算能量/受力/属性，不需要特意提出来白走趟弯 ...

老师，我的轨迹输出导出的是dcd格式，再用REFTRAJ重跑轨迹时，就提示报错Invalid variable name ${} in file.是不是REFTRAJ任务不支持输入dcd格式的轨迹文件啊

作者
Author: Uus/pMeC6H4-/キ 时间: 2025-5-24 11:19

JCenter 发表于 2025-5-24 03:33
老师，我的轨迹输出导出的是dcd格式，再用REFTRAJ重跑轨迹时，就提示报错Invalid variable name ${} in f ...

似乎REFTRAJ明确支持的轨迹文件只有xyz格式，如果要重跑dcd格式的轨迹可能得先自行转换为xyz（以及如果盒子尺寸改变，还有cell）。

作者
Author: JCenter 时间: 2025-5-24 14:22

Uus/pMeC6H4-/キ发表于 2025-5-24 11:19
似乎REFTRAJ明确支持的轨迹文件只有xyz格式，如果要重跑dcd格式的轨迹可能得先自行转换为xyz（以及如果盒 ...

昨晚已经试了，确实是

。我直接用vmd将dcd格式转换成了xyz格式后，再进行REFTRAJ任务时就能直接运行了。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)