KineticPy：长时间动力学计算能量景观中的破碎遍历性的开源软件

软件X 11（2020）100393KineticPy：一种计算具有破碎遍历性的能量景观中的长时间动力学的工具一红Z.科林·毛罗Wilkinson，John C.毛罗·萨奇宾夕法尼亚州立大学材料科学与工程系，University Park，PA 16802，USAar t i cl e i nf o文章历史记录：接收7十一月2019收到修订版2019年12月18日接受2019年12月18日保留字：能源景观焓景观动力学主方程组破态遍历统计力学a b st ra ct能量景观是描述原子尺度复杂系统动力学的有力工具。然而，由于动力学过程的时间尺度差异很大，为了计算性质在长时间内的演变（例如，实验的）时间尺度。KineticPy是一个基于Python的程序，旨在计算能源景观（或焓景观）中每个流域的职业概率的演变，利用不同的时间尺度，该系统利用破遍历的框架来加速计算。考虑到能源前景的相关参数（即，一系列固有的结构，过渡点，和相应的属性值）和所需的温度路径，KineticPy计算任意时间尺度上景观的动力学。KineticPy提供了一个统一的平台，用于探索能源景观和计算任何感兴趣的系统的属性值的演变。KineticPy是开源的，以模块化的方式结构，以方便将来对代码的扩展©2019作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0永久链接https://github.com/ElsevierSoftwareX/SOFTX_2019_346法律代码许可证GNU GPL v3代码版本控制系统[无]软件语言和模块Python：numpy和ABC要求https://github.com/Mauro-Glass-Group/KineticPy链接到手册[无]技术支持邮箱：yzm158@psu.edu1. 介绍Stillinger提出了使用能量景观方法解耦化学系统的振动和构型动力学的想法[1能量景观是对应于原子系统的每个可用结构的能量的表示。配置被分组为盆地，其中盆地中的每个配置通过最陡下降到相同的能量最小值而最小化，称为固有结构。连续的能源景观，然后可以映射到一个离散的固有结构和相关的过渡*通讯作者。电子邮件地址：jcm426@psu.eduJ.C. Mauro）。https://doi.org/10.1016/j.softx.2019.100393连接这些固有结构的点，并且系统的动力学和热力学可以从这组固有结构和转变点计算。Stillinger根据尝试克服跃迁障碍的次数描述了振动贡献，而构型变化则由实际克服障碍以探索景观中新盆地的概率决定[5能源景观方法提供了对各种系统的深入了解，例如，[10][11] [12][13][14][15这种方法可以从势能景观（其中体积是固定的）扩展到焓景观，以考虑在等压条件下系统的动力学[11]。焓景观的一个特别成功的应用是在玻璃弛豫的研究中，2352-7110/©2019作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx2Y.Z. Mauro，C.J.威尔金森和JCMauro/SoftwareX 11（2020）100393∑（）th−1KTDTJexpQJ我KT两个盆地之间的过渡速率由过渡态理论给出，结合玻尔兹曼概率因子，JKT在非等松弛理论的发展、假想温度的重新解释以及理解非单调松弛行为的起源方面都是至关重要的[12在能量景观框架中有两种计算系统动力学的一般方法：随机技术，例如动力学蒙特卡罗和具有内置统计平均的确定性技术，例如本文实现的主方程求解器。动力学蒙特卡罗（KMC）是一种常用的技术，它依赖于绘制出许多可能的振动路径，然后使用每个路径的激活势垒和玻尔兹曼统计来确定将遵循哪条路径[4，17，18]。这是一个强大的技术，用于动态探索能量景观，但它只能沿着一个可能的实现在同一时间传播系统另一种方法是通过求解主方程系统来确定性地计算景观中流域占用概率的演变[13，19]。如果系统的相关固有结构和转变点已经使用当前可用于此目的的软件包之一绘制，则用于求解能量景观动力学的主方程方法特别强大[4，20然而，没有通用的软件，有效地执行动态计算中映射的景观。KineticPy为解决这个问题提供了一个开源解决方案另一种方法是离散路径采样技术，它依赖于简化的两态近似[24]。KineticPy使用随时间变化的盆地占用概率集pi计算系统的任意属性（Φ）的演变，Φ（t）=pi（t）φi，（ 1.1）我其中φi是与每个盆地i相关联的属性值。为了计算所得到的概率，软件迭代地求解每个成对盆地集合的主方程dpi=∑（pjKji−piKij）。（一、二）j=i如果概率是玻尔兹曼分布的，则通常用于求解主方程的数值迭代是冗余的，并且可以被跳过以节省计算成本。当系统冷却时，最高能量势垒不再允许快速平衡，即，概率不再是玻尔兹曼分布，因为实验时间尺度小于在最高激活势垒上平衡所需的时间。然而，由低激活势垒连接的景观区域将保持玻尔兹曼分布在限于那些特定微观状态的尽管整个系统是非遍历的，但它仍将具有满足遍历性的局部区域。这种由非遍历（慢）跃迁隔离的能量景观的划分允许我们定义一个“元基”然而，与实验时间尺度相比，元基之间的转换是缓慢的，因此这些较慢的转换导致元基之间的遍历性的损失。因此，如果我们定义一个实验的时间尺度，我们可以通过求解一组减少的主方程，捕获的元基间的动力学，同时分布的概率元基内的玻尔兹曼分布的元基的限制合奏减少计算复杂性的数值解。因此，元基内的转变明显快于作为温度的函数而变化的阈值速率（Kth）。该阈值速率可以写为，例如，作为DTK=dt·10 K（1.5）其中T是温度，t是时间，阈值的特定值（与温度的时间导数成比例，以10 K−1为尺度，其中K是开尔文）来自参考文献[13]。Kth的值的单位为s−1（秒的倒数）。这使我们可以写出在元基j中占据盆地i的概率作为玻尔兹曼分布和占据给定元基概率的乘积，f[U]这里，Kij是盆地i到盆地j的过渡速率，pi是当前占据每个盆地的概率，因此，∑p i（t）= 1。（一、三）其中f j是占据给定元基j的概率，U i是变形基j中固有结构的能量，Qj是元基的受限配分函数由下式给出iQ=∑exp[−Ui]（1.7）Kij =νi exp[− Uij −Ui].（1.4）由方程式（1.7）对给定元基j内的所有固有结构i求和。 这使我们能够解决一个简化集元基之间的主方程。完整的数学-其中T是绝对温度，k是玻尔兹曼该方法的优点，如Mauro等人[13]所提出的，是利用了一些过渡的局部遍历性质[12，16，25]。当一个系统被淬灭时，遍历性会连续崩溃，因为随着温度的降低，跃迁会逐渐变慢[12，13，25]。这种控制能量景观动态的破态遍历性概念，特别是作为玻璃化转变温度起源的描述，可以追溯到Palmer的早期工作[25]，后来由Mauro等人扩展和推广。在高温下，跃迁比相关的实验时间尺度更快，因此产生的occu-occupational概率是玻尔兹曼分布的。如果所得Mauro等人[13]先前已经发表了数学描述，使得能够有效地计算任何任意系统的动力学。2. 软件描述KineticPy是一个使用Scipy，Numpy和ABC开发的Python脚本，它被设计用于运行带有输入脚本和Python解释器的主软件模块。它需要与KineticPy.py相同的工作目录中的5个输入文件，所有文件都具有相同的文件名，但具有不同的扩展名（所有能量都以电子伏特为单位*.min文件包含两列，即，盆地数和盆地处的能量prp文件包含两列：盆地编号和与盆地对应的属性（如摩尔体积··pi=−（1.6）i∈j.Y.Z. Mauro，C.J.威尔金森和JCMauro/SoftwareX 11（2020）1003933（（））−−◦−deg文件包含两列：盆地编号和对数简并度*.trs文件包含三列：转换点/状态（第1列）、转换点处的能量（第2列）、连接到转换点的左侧盆地（第3列）和右侧盆地（第4列*.eig文件包含三列：从一个盆地（第1列）到另一个盆地（第2列）的特征值（第3列）补充信息中列出了这些文件的示例集。a. 软件构架KineticPy程序中有四个模块、子程序或Python脚本，即Constant.py、Tempera-ture.pySimulation.py、www.example.com和KineticPy.py。以下是这些模块的简要说明。网址：Constant.py该模块为用户提供了一些常用的物理常数，如阿伏伽德罗Temperature.py:在该模块中，五个温度路径，即线性、指数、失超保持线性、固定和锯齿形，被定义为五个类别，即，类LinearTemperature、ExponentialTemperature、QuenchHoldLinearTemperature 、 FixedTemperature 和ZigzagTemperature，re-trigger。用户还可以从文件（类FromFileTemperature）中读取自己的温度路径。所有这些类都是从一个抽象类Temperature派生的。Temperature 类 继 承 自 ABC （ AbstractBaseClasses），abstractmethod装饰器用于声明每个方法的抽象。Simulation.py:此模块导入常量和温度模块，包括主要的计算。KineticPy.py:KineticPy.py导入Simulation模块。程序创建了Simulation类的一个实例，并要求用户提供输入参数文件的名称。然后，它调用类并要求输入输出文件的名称。当类的一个实例被创建时。传递给方法的参数用于初始化类的属性，这些属性可以在任何类方法中使用对于给定时间（inTime），温度路径类中的Evaluate方法计算给定时间的温度，SafeTimeStep方法返回安全时间步长。每个温度路径由其自身的参数确定/定义。参数和温度路径公式描述如下：线性温度：线性温度路径由三个参数定义，即，初始温度、最终温度和总冷却时间。在给定的时间，Evaluate的输出温度由下式给出：初始温度总冷却时间内的最终温度-初始温度《时代周刊》，w这里0<在<总冷却时间中的时间。指数温度：指 数 温 度 路 径 由 三 个 参 数 定 义 ， 即 ，inInitialTemperature、 in-FinalTemperature 和inDecayTime。在给定的时间，Evaluate的输出温度由下式给出：初始温度-（in最终温度 −inInitialTemperature）1 .一、0exp in Time，inDecayTime这里的时间> 0。淬火保持线性温度：淬火保持线性温度路径由四个参数定义，即，inQuenchTempera-ture、inHoldTime、inFinalTemperature和inCooling- ingTime。在给定的时间，Evaluate的输出温度为◦ inFinalTemperature，如果inTime>（inCoolingTime+inHoldTime）inQuenchTemperature+计算结果将保存到该文件中。InFinalTemperatureInQuenchTemperature冷却时间（inTimeb. 软件功能温度模块：在该模块中，五个温度路径，即线性、指数、失超保持线性、固定和锯齿形，被定义为五个类别，即，类LinearTemperature、ExponentialTemperature、QuenchHoldLinearTemperature 、 FixedTemperature 和ZigzagTemperature，re-trigger。用户还可以从文件（类FromFileTemperature）中读取自己的温度路径。所有这些类都是从一个抽象类Temperature派生的。Temperature 类 继 承 自 ABC （ AbstractBaseClasses），abstractmethod装饰器用于声明每个方法的抽象。每个类包含三个方法：init、Evaluate和SafeTimeStep。Python中类的init方法的作用类似于面向对象编程概念中类的构造函数。调用此方法−inHoldTime），如果inTime>inHoldTime◦ 淬火温度，其他固定温度：线性温度路径由param- eterinTemperature定义。在给定的时间inTime，Evaluate的输出温度由in-Temperature给出。曲折温度：Z 形 温 度 路 径 由 九 个 参 数 定 义 ， 即 ，InInitialTemperature、InLowTemperature、InHighTemperature、InFinalTemperature、InInitialCoolingTime、InReheatingTime、InFinalCoolingTime、InLowHoldTime和InHighHoldTime。在给定的时间内，Evaluate的输出温度为◦ inInitialTemperature，如果inTime0<········+·····4Y.Z. Mauro，C.J.威尔金森和JCMauro/SoftwareX 11（2020）100393+−−）−−++公司简介◦−=初始温度（初始温度中的低温度）初始冷却时间if inTime inInitialCoolingTime<《时代周刊》，TotalCoolingTime=600#9000;#ReheatingTime= 20;return20;◦ inLoww Temperature， if inTime（ inInitialCoolingTime+ inLowwHoldTime）◦ 低温度+return0;return0;（在高温和低温下）再热时间（inTime−模拟时间=600#9040;inInitialCoolingTime- inLowHoldTime，如果inTime（inInitialCoolingTime+低保持时间+再加热时间）◦ inHighTemperature， if inTime（ inInitialCoolingTime+ inLowHoldTimereturn1;return71;int count= 64;原子质量= 78。九十六;return0;+再加热时间+高保持时间）◦ 产品介绍HopDistance=2。0e−10;PrintPhaseDistributions=0;（最终温度）在高温下）inFinalCoolintTime（inTime−读取上述输入参数文件后初始冷却时间-低保持时间在中，温度路径（mTemper-−如果inReheatingTime（inHighHoldTime）、inTime（inInitialCoolingTimeaturePath）然后使用上述参数创建和初始化。在低保持时间在再加热时间在高保持时间在最终冷却时间）inFinalTemperature， other wise模拟模块：在仿真在五个输入文件中定义该系统的势能景观，其中文件名在输入参数文件中的变量MinimumFile、PropertyFile 、De generacyFile 、TransitionFile 、EigenvalueFile中指定。这五个文件的示例列在https：//github中。com/Mauro-Glas同学们， 也就是说，init方法，平衡代谢S-Group/KineticPy/Example。方 法 、 MetabasinProbability 方 法 、 CalculateInter-MetabasinRates 方 法 、 DetermineMetabasins 方 法 和MasterEquationSolver方法。init方法：语句创建Simulation类的实例时，程序读入这五个文件并初始化属性，例如：变量/ 数 组 / 矩 阵 ， 例 如 mMinimumNumbers 、mEnergyMatrix 、 mRecurtiveEnergyMatrix 、mPropertyVector、mFrequencyMatrix、mLogDegeneracy等。模拟（inParameterFileName），init模拟类。这些描述方法会自动调用。它以参数文件名作为参数。这个输入文件包含进行计算所需的参数.下面是一个输入文件Se_test5的例子（程序可以为输入文件使用任何用户定义的Python有效文件名）。，表示采用线性冷却路径（温度路径=线性）从（初始温度=）800 K到（最终温度=）200 K，总计（总冷却时间=）600 s。Se_test5SimulationName= Se_test5;#ExecutionMode= Production;MinimumFile= Se_test5.min;TransitionFile= Se_test5.trs;EigenvalueFile= Se_test5.eig;PropertyFile= Se_test5.prp;DegeneracyFile=Se_test5.deg;InitialTemperaturePath= Linear; #FileName=TemperaturePath.txt;#LowTemperature= 200;#高温=270;属性如下：mMinimumNumbers ： 此 numpy 数 组 由[0，mNumBasins）范围内的均匀间隔整数值组成，其中mNumBasins是系统中盆地的数量mEnergyMatrix：2-D对称numpy矩阵（mNumBasins X mNumBasins）。 对角元素是每个盆地的最小能量，非对角元素是从一个盆地转换到另一个盆地所需的能量，所有这些都乘以常数kElectronVolt。出于计算目的，所有这些值都向上移动所有盆能量的最小能量（mBiasEnergy）。因此，盆地能量的最小值为零。mLogDegeneracy：包含每个盆地的对数简并ln（g i）值的mNumBasins元素numpy数组mNumEnergy Matrix：包含有效自由能H ikTln（gi）的（mNumBasins X mNumBasins）numpy矩阵。mPropertyVector：包含每个盆地属性值的mNumBasins元素numpy数组mFrequencyMatrix ： （ mNumBasins XmNumBasins）numpy矩阵，包含由下式给出的变化频率：return200;◦ 其中b是本征值，B2int maximum= 9000;m=π原子质量。KAvogadros编号：1000◦··◦◦◦◦◦MY.Z. Mauro，C.J.威尔金森和JCMauro/SoftwareX 11（2020）1003935（）下一页∑fαexpαA∈αmPhaseSpaceVector：mNumBasins-元素numpy阵列，其包含占据由参考文献[ 13 ]中的修改的等式（32）给出的任何微观状态的概率。为例如，在时间0，初始温度为T0（开尔文），概率为根据[ 13 ]中的第4节中的算法在时间步的开始处重新分配元基内概率，并且根据参考文献[14]中的等式32重新分配元基内概率。[13]：pA=Q−，其中F A =H i∈A −kTln（gA）为由pi（0）=1 exp（−Hi−kT0ln（gi））给出，状态A的自由能和分配函数QkT0QA=∑exp-Hi∈A−kTln（gA）. 亲-w（此处为pa）约束函数Q=--iexpA∈αkT0H ikT0 ln（ gi）kT0巴辛岛◦Hi是Gram使用这些概率来计算-计算有效基间转换速率从metabasinα到metabasinβ，即，ωαβ=mNumMetabasins：数量变基初始值设置为1。nβ∑fαKA，B∈βPA.如果任何一种间基mMetabasin率：amNumMetabasins-元素列表，包含将能量景观划分为元基之后的转换速率。mNumMetabasins在程序中动态列表中的每个元素都是一个动态大小为（mNumBasins X mNumBasins）的2D numpy矩阵。mMetabasinRank ： mNumMetabasins 元素numpy数组，包含每个分区元库的转换速率的等级。mMetabasinIndices：一个mNumMetabasins元素列表，包含元数据库内/间的交换率.该列表中的每个整数元素由一个numpy数组组成，该数组具有每个分区元数据库中的盆地编号平衡代谢方法：该方法将温度（inTempera- ture）作为其参数，并计算占据每个元库的等概率。类Simulation的属 性 ， 即 ， mPhaseSpaceVector 已 更 新 到 位 。Metabasin概率方法：该方法将元库的索引（inIndex）作为其参数，并返 回 计 算 出 的 占 用 此 元 库 的 概 率 。CalculateInterMetabasinRates方法：该方法将温度（inTempera- ture）作为其参数并计算元基间速率。类Simulation的属性，即，mInterMetabasinRates已更新到位。测定代谢物方法：该方法以温度（inTemperature）和时间步长（inTimeStep）为自变量，根据参考文献[13]中描 述 的 划 分 算 法 ， 在 该 温 度 和 时 间 步 长 类Simulation 的 属 性 ， 即 ， mMetabasinRates 、mMetabasinRank和mMetabasinIndices 将在适当位置更新。MasterEquationSolver方法：如名称所示，该方法求解系统的主方程动力学，使用参考文献[13]第6节中描述的算法的预定义温度路径。程序在温度路径上循环，在每个时间步，它选择一个时间步dt，在该时间步上温度可以假定为常数。它将系统划分为多个元数据库转换速率低于阈值速率，将使用参考文献[13]第5节中描述的近似解。重复此划分过程，并在时间步长上对主方程的缩减集进行积分。整个算法在每一个新的步骤中重复，结果保存在输出文件中，其名称由用户提供。c. 常见错误用 户 可 以 下 载 四 个 .py 文 件 ， 即 ， Constants.py 、 Tem-perature.py 、 Simulation.py 和 KineticPy.py 文 件 。 主 程 序KineticPy可以通过点击Run（绿色箭头）按钮在Spyder（通过Anaconda3的科学Python开发环境）中运行。用户还可以右键单击KineticPy.py文件并选择“空闲编辑”。在在弹出的IDLE窗口中选择Run/Run Module来运行KineticPy.py脚本。一个常见的错误是无法导入numpy。用户需要将numpy导入到运行代码的Python中，可以通过以下操作完成：打开命令窗口类型在KineticPy IDLE窗口中选择Run/Run Module以运行程序3. 说明性示例在软件存储库的示例文件夹中提供了先前报告的硒景观的示例输入脚本和相关文件[26]。 该模拟使用线性淬火模块在700 s的时间内将系统从800 K线性淬火到100 K（图1）。1.一、当系统被淬灭时，同时计算在相同温度下的变基数目、当前焓和平衡系统取焓曲线相对于温度的数值导数，可以记录相应的该模型显示了对玻璃化转变和量热法之间关系的见解，量热法是用于量化玻璃化转变的最常用技术。在参考文献中可以找到探索该软件的更多[3、5、6、16、19、264. 影响这个开源软件提供了探索和计算能源景观中占据盆地的概率演变的能力，通过利用破碎遍历性和可变元基计数的概念的加速算法代码的输出包括任何感兴趣的任意属性的演化。该代码将允许精确模拟玻璃化转变/弛豫，◦◦◦◦········FAKT6Y.Z. Mauro，C.J.威尔金森和JCMauro/SoftwareX 11（2020）100393Fig. 1. （A）平衡和非平衡系统在玻璃熔化期间的系统能量。(B)平衡态和非平衡态能量的热容和（C）淬火过程中的变基数目蛋白质折叠，也可用于预测催化反应行为。玻璃弛豫和蛋白质折叠都是复杂的扩散驱动过程，为了能够预测它们的物理意义，必须探索景观;该软件为获得这些物理见解提供了桥梁该软件的模块化设计允许在发现新物理和开发新景观方法时轻松扩展该代码已被广泛用于玻璃转变和弛豫行为的研究，包括开发玻璃非平衡粘度的焓变模型，改进对临界温度概念的物理理解，以及揭示玻璃中波动弛豫的固有非单调性质[125. 结论我们报告了一个开放源代码，探索能源景观（或焓景观）在很长的时间尺度。KineticPy允许在任何任意的能源景观中显式计算长时间动力学，利用破碎遍历系统的不同时间尺度特性。感兴趣的具体应用包括蛋白质折叠和玻璃转变/弛豫。该软件的简单性和可扩展性将使能源景观动力学的研究更加深入。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作附录A. 补充数据与本文相关的补充材料可以在https://doi.org/10.1016/j.softx.2019.100393上找到。引用[1]斯蒂林杰 化学物理杂志1988;88：7818-25.[2]Stillinger FH，Weber TA. Phys Rev A1982;25：978-89.[3]Mauro JC，Loucks RJ，Balakrishnan J，Varshneya AK. J Non CrystSolids2007;353：1268-73.[4]MousseauN，Barkema GT. Phys Rev E 1998;57：2419[5]Li J，Kushima A，Eapen J，Lin X，Qian X，Mauro JC，Diep P，Yip S. PLoSOne2011;6：1-7.[6]Prada-Gracia D，Gómez-Gardeñes J，Echenique P，Falo F. PLoS ComputBiol2009;5.[7]瑙米斯 Phys Rev E 2005;71：1-7.[8]El-MellouhiF，Mousseau N，Lewis LJ.Phys Rev B 2008;78：1[9]豪雅A 20.第20章. 373101。[10]放大图片作者：Benjamin SJ，Hammes GG，Hammes-Schiffer S.生物化学2008;47：3317 - 21.[11]MauroJC，Loucks RJ，Balakrishnan J. J Phys Chem B 2006;110：5005[12]MauroJC，Gupta PK，Loucks RJ.第126章. 184511[13]MauroJC，Loucks RJ，Gupta PK.J Phys Chem A 2007;111：7957[14]MauroJC，Loucks RJ，Gupta PK.J Am Ceram Soc 2009;92：75[15]MauroJC，Uzun SS，Bras W，Sen S.102.第一百零二章：155506[16]放大图片作者：MauroJC，Allan DC，Potuzak M.物理学修订版B2009;80.094204。[17]Weinkamer R，Fratzl P，Gupta HS，Penrose O，Lebowitz JL. PhaseTransit2004;77：433-56.[18]Fichthorn KA，Weinberg WH. J Chem Phys1991;95：1090-6.[19]PollakE，Auerbach A，Talknerz P. Biophys J 2008;95：4258[20]吴晓波，李晓波，李晓波.J Chem Phys 2000;113：9901[21]Barettin D，Sibani P. Phys Rev E2011;84：1[22]NeelamrajuS，Oligschleger C，Schön JC.第147章.[23]威尔士DJ，Scheraga HA。科学（80-.）1999;285：1368-72.[24]第二十四话Mol Phys 2002;100：3285[25日]Palmer RG. Adv Phys 1982;31：669[26]MauroJC，Loucks RJ.物理学修订版B 2007;76.174202.[27]瑙米斯 GG. J 非 Cryst 固体2006;352：4865-70.[28]SmolinN，Daggett V. J Phys Chem B 2008;112：6193