CN114220489A - 原子结构弛豫的非单调线搜索方法及装置 - Google Patents

原子结构弛豫的非单调线搜索方法及装置 Download PDF

Info

Publication number
CN114220489A
CN114220489A CN202111534901.5A CN202111534901A CN114220489A CN 114220489 A CN114220489 A CN 114220489A CN 202111534901 A CN202111534901 A CN 202111534901A CN 114220489 A CN114220489 A CN 114220489A
Authority
CN
China
Prior art keywords
line search
accept
tentative
determining
trial
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111534901.5A
Other languages
English (en)
Other versions
CN114220489B (zh
Inventor
高兴誉
刘歆
宋海峰
胡雨宽
方俊
杨真
赵亚帆
王丽芳
刘海风
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
INSTITUTE OF APPLIED PHYSICS AND COMPUTATIONAL MATHEMATICS
Academy of Mathematics and Systems Science of CAS
Original Assignee
INSTITUTE OF APPLIED PHYSICS AND COMPUTATIONAL MATHEMATICS
Academy of Mathematics and Systems Science of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by INSTITUTE OF APPLIED PHYSICS AND COMPUTATIONAL MATHEMATICS, Academy of Mathematics and Systems Science of CAS filed Critical INSTITUTE OF APPLIED PHYSICS AND COMPUTATIONAL MATHEMATICS
Priority to CN202111534901.5A priority Critical patent/CN114220489B/zh
Publication of CN114220489A publication Critical patent/CN114220489A/zh
Application granted granted Critical
Publication of CN114220489B publication Critical patent/CN114220489B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/20—Identification of molecular entities, parts thereof or of chemical compositions

Landscapes

  • Theoretical Computer Science (AREA)
  • Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种原子结构弛豫的非单调线搜索方法及装置，使用BB步长作为试探步长，从而能够充分利用其优越于SD的局部收敛速度，以及使用非单调线搜索策略来判定试探位置是否被接受，减少调用线搜索算法的次数，从而极大地减少线搜索带来的昂贵开销。

Description

原子结构弛豫的非单调线搜索方法及装置
技术领域
本发明涉及原子结构弛豫技术领域，具体而言，涉及一种原子结构弛豫的非单调线搜索方法及装置。
背景技术
给定某一原子体系，其原子弛豫的优化问题可以简单地表述为
Figure BDA0003412809360000011
其中，R:＝[R1,R2,…,RN]，Ri为第i个原子的三维Cartesian坐标，N>0为原子个数；A＝[a1,a2,a3]，aj为单胞的第j个基矢(下文中称A为单胞)；E为体系能量泛函；Ω为问题(1.1)的可行域，其取法由应用场景决定。
考虑求解固定单胞的问题(1.1)，假设给定单胞A恒等于某个A0，此时问题(1.1)变为
Figure BDA0003412809360000012
从形式上看，问题(1.2)(对变量R)是无约束优化问题。简记E(R):＝E(R,A0)，记E对R的负梯度为
Figure BDA0003412809360000013
也即原子的受力。
求解问题(1.2)的方法可分为两大类：基于力的方法(force based methods)与基于分子动力学的方法(molecular dynamics based methods)。基于力的方法可以分为一阶方法、二阶方法等，其中，梯度法、共轭梯度(conjugate gradient,CG)法、牛顿法、拟牛顿(quasi-Newton,QN)法均属于线搜索型算法。
现有的基于力的原子弛豫方法均存在以下缺点：
i)基于适当线搜索的一阶方法尽管具有全局收敛性，但在基于第一原理的原子弛豫优化的背景下，线搜索十分昂贵；另外，求解一般非凸问题时，一阶方法的在能量曲面极小点附近收敛速度(即，局部收敛速度)较慢；
ii)尽管预条件对CG具有潜在的加速效果，但设计普适的预条件子极其困难。早期的工作针对个别体系设计了预条件子，而近几年的工作则集中在设计合适的代理势；后者往往还需要拟合代理势的参数。这些都限制了预条件应用的范围；
iii)二阶方法虽然在一定条件下具有局部超线性(甚至二次)收敛速度，但往往以空间复杂度或计算复杂度为代价。例如，牛顿法需要显式地计算、存储Hessian及其逆，信赖域(trust-region,TR)则需要近似地求解子问题。另外，求解一般的非凸问题时，(拟)牛顿法一般不保证从远离平衡位置出发迭代的收敛性(即，全局收敛性)。这限制了它们在远离平衡位置体系上的应用。
发明内容
为解决上述问题，本发明提供一种原子结构弛豫的非单调线搜索方法，包括以下步骤：
步骤1，输入初始原子位置R0，并令k:＝0；
步骤2：判断原子受力Fk是否满足收敛准则。若是，则终止；若否，则执行步骤3；
步骤3，更新搜索方向Dk＝Fk；
步骤4，计算Barzilai-Borwein步长；
步骤5，根据所述Barzilai-Borwein步长计算试探位置
Figure BDA0003412809360000021
步骤6，根据非单调线搜索策略判断是否接受所述试探位置
Figure BDA0003412809360000022
若接受，则执行步骤8；若不接受，则执行步骤7；
步骤7，执行线搜索算法，在满足线搜索终止条件的情况下执行步骤8；
步骤8，更新原子位置
Figure BDA0003412809360000023
并令k:＝k+1；执行所述步骤2。
可选地，基于以下公式计算得到两个用于加速梯度法的Barzilai-Borwein步长：
Figure BDA0003412809360000031
其中，S为前后原子位置之差，Y为后前原子受力之差；
或者，基于以下公式确定Barzilai-Borwein步长：
Figure BDA0003412809360000032
可选地，所述根据非单调线搜索策略判断是否接受所述试探位置
Figure BDA0003412809360000033
包括：判断试探位置能量值的连续上升次数是否大于或等于预设次数阈值；若否则确定接受所述试探位置
Figure BDA0003412809360000034
若是则确定不接受所述试探位置
Figure BDA0003412809360000035
或者，根据预先储存的最近m个试探位置的能量值
Figure BDA0003412809360000036
判断当前的所述试探位置
Figure BDA0003412809360000037
的能量值是否充分小于
Figure BDA0003412809360000038
若是则确定接受所述试探位置
Figure BDA0003412809360000039
若否则确定不接受所述试探位置
Figure BDA00034128093600000310
或者，根据预先储存的试探位置的能量值，判断当前的所述试探位置
Figure BDA00034128093600000311
的能量值是否小于预先存储的所有能量值的凸组合，若是则确定接受所述试探位置
Figure BDA00034128093600000312
若否则确定不接受所述试探位置
Figure BDA00034128093600000313
可选地，所述执行线搜索算法，包括：计算试探步长及计算试探位置，以及判断是否满足线搜索终止条件；若是，则执行所述步骤8；若否，则再次执行所述步骤7。
可选地，若非单调线搜索策略采用判断试探位置能量值的连续上升次数大于或等于所述预设次数阈值，则在执行所述执行步骤7之前，所述方法还包括：将当前的原子位置替换为能量值连续上升之前的上一试探位置。
可选地，所述线搜索算法包括：Brent算法、插值方法或回溯方法。
本发明提供一种原子结构弛豫的非单调线搜索装置，包括：
初始输入模块，用于输入初始原子位置R0，并令k:＝0；
收敛判断模块，用于判断原子受力Fk是否满足收敛准则。若是，则终止；若否，则触发方向更新模块运行；
方向更新模块，用于更新搜索方向Dk＝Fk；
步长计算模块，用于计算Barzilai-Borwein步长；
试探位置计算模块，用于根据所述Barzilai-Borwein步长计算试探位置
Figure BDA0003412809360000041
非单调线搜索判断模块，用于根据非单调线搜索策略判断是否接受所述试探位置
Figure BDA0003412809360000042
若接受，则触发位置更新模块运行；若不接受，则触发线搜索执行模块运行；
线搜索执行模块，用于执行线搜索算法，在满足线搜索终止条件的情况下触发位置更新模块运行；
位置更新模块，用于更新原子位置
Figure BDA0003412809360000043
并令k:＝k+1；触发所述收敛判断模块运行。
可选地，基于以下公式计算得到两个用于加速梯度法的Barzilai-Borwein步长：
Figure BDA0003412809360000044
其中，S为前后原子位置之差，Y为后前原子受力之差；
或者，基于以下公式确定Barzilai-Borwein步长：
Figure BDA0003412809360000045
可选地，所述非单调线搜索判断模块，具体用于：判断试探位置能量值的连续上升次数是否大于或等于预设次数阈值；若否则确定接受所述试探位置
Figure BDA0003412809360000046
若是则确定不接受所述试探位置
Figure BDA0003412809360000047
或者，根据预先储存的最近m个试探位置的能量值
Figure BDA0003412809360000048
判断当前的所述试探位置
Figure BDA0003412809360000049
的能量值是否充分小于
Figure BDA00034128093600000410
若是则确定接受所述试探位置
Figure BDA00034128093600000411
若否则确定不接受所述试探位置
Figure BDA00034128093600000412
或者，根据预先储存的试探位置的能量值，判断当前的所述试探位置
Figure BDA0003412809360000051
的能量值是否小于预先存储的所有能量值的凸组合，若是则确定接受所述试探位置
Figure BDA0003412809360000052
若否则确定不接受所述试探位置
Figure BDA0003412809360000053
可选地，所述线搜索执行模块，具体用于：计算试探步长及计算试探位置，以及判断是否满足线搜索终止条件；若是，则触发所述位置更新模块运行；若否，则再次调用所述线搜索执行模块。
本发明实施例使用BB步长作为试探步长，从而能够充分利用其优越于最速下降法的局部收敛速度，以及使用非单调线搜索策略来判定试探位置是否被接受，减少调用线搜索算法的次数，从而极大地减少线搜索带来的昂贵开销。相对于现有基于力的原子弛豫方法，本发明实施例提供的非单调线搜索方法属于一阶方法，但是基于BB步长提高了局部收敛速度；相对于预条件对CG具有潜在的加速效果，其不需要克服设计普适预条件子的困难；相对于二阶方法，其空间复杂度或计算复杂度也较低。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1为现有的线搜索型算法流程示意图；
图2为本发明实施例提供的一种原子结构弛豫的非单调线搜索方法流程示意图；
图3为本发明实施例提供的CG与NBB在55个算例上的性能分析示意图；
图4为本发明实施例提供的一种原子结构弛豫的非单调线搜索装置的结构示意图。
具体实施方式
为使本发明的上述目的、特征和优点能够更为明显易懂，下面结合附图对本发明的具体实施例做详细的说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
以下首先介绍基于力的方法。在叙述算法时，采用上标表示迭代指标，用下标表示分量；例如
Figure BDA0003412809360000061
表示经过第k步迭代，第j个原子位置的第i个分量；相应地，Fk则表示F在Rk处的取值。
(1)一阶方法
一阶方法的迭代格式通常为
Rk+1:＝Rk+αkDk, (1.3)
其中αk>0为步长，Dk为搜索方向。求解问题(1.2)常用的一阶方法包括梯度法和共轭梯度CG法。
在(1.3)式中，对
Figure BDA0003412809360000062
梯度法取Dk:＝Fk；另外，若
Figure BDA0003412809360000063
则称为最速下降(steepest descent，SD)法；而若(1.4)式中的极小化问题被非精确求解，则称为梯度下降(gradient descent,GD)法。
CG在k＝0时与GD相同，而在k>0时，取
Dk:＝Fk+βkDk-1， (1.5)
Figure BDA0003412809360000064
其中βk为共轭参数。βk的常用选取方式包括以下四种：
Figure BDA0003412809360000065
Figure BDA0003412809360000066
Figure BDA0003412809360000071
Figure BDA0003412809360000072
这里
Figure BDA0003412809360000073
MijNij是Frobenius内积，
Figure BDA0003412809360000074
是Frobenius范数。
求解(1.4)、(1.6)式中极小化问题的过程称为线搜索；如果极小化问题被精确求解(如(1.4)式所表示)，则被称为精确线搜索，反之为非精确线搜索(如(1.6)式所表示)。在原子弛豫优化时，精确线搜索是极其昂贵的，甚至是不可能的。此时，往往依照一些停机条件做非精确线搜索。常用的线搜索方法包括Brent算法、插值(interpolation)方法与回溯(backtracking)方法；常用的停机条件包括Armijo条件、Goldstein条件与(强)Wolfe条件。
当目标函数是强凸二次函数时，可以证明SD与CG的收敛速度依赖于目标函数Hessian的条件数。此时，设计合适的预条件子(preconditioner)可以改善优化问题的条件数，进而改进算法的收敛速度。尽管问题(1.2)的目标函数E不是强凸二次的，但现有的一些工作表明，预条件能够明显地改善CG求解问题(1.2)时的数值表现。
现有的使用预条件共轭梯度法(PCG)求解问题(1.2)的工作从不同的角度近似E的Hessian，并以此作为CG的预条件子。下文中“加速比”指PCG与CG停机分别所需迭代步数的比值。
早期的工作根据所处理体系的特点直接构建Hessian的近似。例如，在纳米体系中，不论原子的数目多少，体系局部的成键环境(local bonding environments)一般较为相似，使用少数几个原子作为结构中的代表元(motif)，之后由代表元的Hessian构造整个体系Hessian的近似，最终以其作为CG的预条件子。数值实验表明，在一些纳米体系中，基于代表元的PCG相较于CG加速比在2～4之间。
近几年则多使用“代理势(surrogate potential)”近似E，并基于此来加速CG。
例1，使用Lennard-Jones(LJ)12-6势作为代理势。在迭代的过程中，有时需要基于F拟合LJ势的参数(这等价于求解一个线性最小二乘问题)，之后将拟合后LJ势的Hessian作为CG的预条件子。数值实验表明，在立方体材料Cd4Se4等体系中，基于拟合LJ势的PCG相较于CG加速比在2～5之间；而在Pt38、Pt54以及一些异质体系中，尽管PCG与CG均达到了最大迭代步数限制，但PCG相较于CG收敛到的能量值更低。
在金属团簇和磁性原子簇中，通往能量曲面局部极小的CG迭代路径往往非常蜿蜒曲折。因此，提出在每一步迭代中，在当前点(如Rk)基于F拟合一个代理势；之后，以代理势作为价值函数，做几步CG迭代得到若干试探点(如{Rk,j})，并依此生成分段线性的搜索曲线；最后，重新以E为价值函数，沿搜索曲线线搜索得到下一步的迭代点。将此类方法命名为CLS。对于金属体系，提出使用Gupta力场作为代理势，并使用启发式算法拟合其中的参数。数值实验表明，在Pt20、Co120、Cu20Au18等金属体系中，CLS相较于CG加速比在2～4之间。
例2，与例1中类似，在每一步迭代中，在当前点(如Rk)根据共价键和金属域对体系能量的贡献(unified valence force field+Gupta potential)，构建力场EFF(这里“FF”指“force field”)；之后，以EFF为价值函数，从当前点出发得到极小点(如
Figure BDA0003412809360000081
)；再以方向
Figure BDA0003412809360000082
按类似于(1.5)式的方式更新搜索方向；最后，重新以E为价值函数，沿搜索方向线搜索得到下一步迭代点。可以验证，当E为强凸二次函数时，该方法等价于直接使用Hessian作为CG的预条件子。其与例1的主要不同之处在于，力场EFF比Gupta更加广泛，并且不需要拟合力场的参数(均选取为默认数值)。大量数值实验表明，在多种体系中，基于EFF的PCG相较于CG加速比在2～9之间。
(2)二阶方法
求解问题(1.2)的二阶方法包括牛顿法、拟牛顿QN法和信赖域TR方法。
牛顿法的迭代格式为
Rk+1:＝Rk+αk(Hk)-1Fk, (1.7)
其中Hk为Rk处E的Hessian。
QN则是将(1.7)式的Hk换成某个近似
Figure BDA0003412809360000091
从
Figure BDA0003412809360000092
出发，QN通过某种格式更新Hessian的近似。最常用的更新格式是Broyden-Fletcher-Goldfarb-Shanno(BFGS)公式。记BFGS拟牛顿法在第k+1步迭代后更新的Hessian近似为
Figure BDA0003412809360000093
它总是满足
Figure BDA0003412809360000094
其中Sk＝Rk+1-Rk、Yk＝Fk-Fk+1分别为前后原子位置、后前原子受力之差。(1.8)式也被称为“割线方程”。为了按(1.7)式更新原子位置，拟牛顿法实际需要计算
Figure BDA0003412809360000095
可以验证，
Figure BDA0003412809360000096
具有显式表达式。另外，BFGS拟牛顿法还有对应的有限内存版本；后者在空间复杂度和每步的计算复杂度上均更优。现有文献中使用拟牛顿法弛豫氮化碳、氮化硼的多种多晶型体系，用以计算它们的结构、黏结与电光性质。
前文所述的方法均为线搜索型的方法。TR与上文所述的所有方法均不相同。在第k+1步迭代中，TR(近似)求解(1.2)的近似子问题
Figure BDA0003412809360000097
得到迭代步Dk，其中
Figure BDA0003412809360000098
为E在Rk处的二次近似，
Figure BDA0003412809360000099
为Hk的近似，‖·‖为某种矩阵范数，Δk>0为信赖域半径。之后，根据下降比
Figure BDA00034128093600000910
判断是否接受Dk以及修改信赖域半径。现有文献中将BFGS公式用于更新信赖域子问题(1.9)中的Hessian近似；另外，使用二次函数近似最近邻原子之间的相互作用势，并以此构建初始Hessian近似。在SiO2、LaCuOS、MgVO3、Bi4Ti3O12等体系进行了实验；数值结果表明，相较于已有的对角初始Hessian近似，设计的初始Hessian近似可以将求解Kohn-Sham方程的次数减少30％。
(3)线搜索型算法流程
现有技术中，梯度法、CG、牛顿法、QN均属于线搜索型算法。它们具有较为统一的流程，图1示出了现有的线搜索型算法流程示意图，包括以下步骤：
S101，输入初始原子位置R0，并令k:＝0。
S102，判断原子受力Fk是否满足收敛准则(例如||Fk||F小于某一充分小的量)。若是，则终止算法；若否，则进入S103。
S103，更新搜索方向Dk。
S104，计算初始试探步长
Figure BDA0003412809360000101
S105，更新试探位置
Figure BDA0003412809360000102
S106，判断是否满足线搜索终止条件。若满足线搜索终止条件，则记
Figure BDA0003412809360000103
并进入S107；否则，调用线搜索算法更新试探步长
Figure BDA0003412809360000104
再次进入S105。
S107，
Figure BDA0003412809360000105
并令k:＝k+1；进入S102。
线搜索型算法一般在搜索方向Dk、初始试探步长
Figure BDA0003412809360000106
线搜索算法及其终止条件的选取上有所不同。
值得注意的是，在图1展示的流程中，一旦初始化(S101)/更新了原子位置(S105)，就需要调用“电子步”算法求解Kohn-Sham方程，以获得当前位置处的能量值和原子受力。这在基于第一原理的原子弛豫优化的背景下是十分昂贵的。因此本发明实施例中将减少S106中线搜索过程调用“电子步”的次数，作为提高线搜索型算法效率的关键。这也意味着要改进初始试探步长(S104)和线搜索策略(S106)。
为克服现有技术求解问题(1.2)时的缺陷，本发明实施例基于技术1：使用Barzilai-Borwein步长的梯度法，以及技术2：一种非单调线搜索策略，对现有线搜索型算法流程进行改进。
本实施例的方法属于梯度法(见技术1)，因此是一阶算法；同时，本实施例的方法使用了线搜索(技术2)，因此也是线搜索型算法。技术1可以为线搜索型算法提供好的初始试探步长，并且使算法拥有较快的局部收敛速度；而技术2则着眼于减少线搜索(即图1中S106)的开销，同时保证算法的全局收敛性。
图2示出了本发明实施例提供的一种原子结构弛豫的非单调线搜索方法流程示意图，包括以下步骤：
步骤1，输入初始原子位置R0，并令k:＝0。
步骤2：判断原子受力Fk是否满足收敛准则。若是，则终止；若否，则执行步骤3。
步骤3，更新搜索方向Dk＝Fk。
步骤4，计算Barzilai-Borwein步长。
Barzilai和Borwein提出了如下两个用于加速梯度法的步长：
Figure BDA0003412809360000111
其中，S为前后原子位置之差，Y为后前原子受力之差。
基于以下公式确定交替Barzilai-Borwein步长：
Figure BDA0003412809360000112
Figure BDA0003412809360000113
统称为Barzilai-Borwein(BB)步长，它们各自分别满足
Figure BDA0003412809360000114
也即，
Figure BDA0003412809360000115
I和
Figure BDA0003412809360000116
I是近似满足割线方程(1.8)的纯量阵。这表明，使用BB步长的梯度法(后文简称为BB方法)虽然是一阶方法，但在某种意义上与QN具有相似性。进一步，还可以使用BB步长、交替BB步长的变体。
BB方法在许多应用中都表现出了相较于SD更快的局部收敛速度，本实施例提出将BB步长应用于原子弛豫优化领域。事实上，可以证明使用BB方法在极小化二维强凸二次函数时具有R-超线性收敛速度；在极小化一般的强凸二次函数时，则有R-线性收敛速度。
在研究中发现，交替BB步长相对于BB步长在数值表现中往往更好，因此在本实施例中可以选用交替BB步长加速梯度法。需要指出的是，在第一步中，BB步长无法定义，故可以使用GD；当<Sk-1,Yk-1><0时，取其绝对值进行后续计算。本实施例中使用交替BB步长作为试探步长，从而能够充分利用其优越于SD的局部收敛速度。
步骤5，根据Barzilai-Borwein步长计算试探位置
Figure BDA0003412809360000121
步骤6，根据非单调线搜索策略判断是否接受试探位置
Figure BDA0003412809360000122
若是，则执行步骤8；若否，则执行步骤7。
考虑到尽管交替BB方法在极小点附近具有较快的收敛速度，但它在全局的收敛性质并没有保障。特别地，BB方法产生的能量值序列{E(Rk)}一般是非单调的。如果从距离极小点较远的位置开始迭代，BB方法可能会发散。为此，本实施例引入了必要的非单调线搜索(non-monotone line search,NLS)策略。
可选地，本实施例中采用的NLS策略主要有两类：最大型(MNLS)与平均型(ANLS)。前者需要储存最近的m个能量值
Figure BDA0003412809360000123
并要求下一个能量值充分小于
Figure BDA0003412809360000124
后者则要求下一个能量值充分小于过往所有能量值的一个凸组合。前者对于m的选取较为敏感；后者在迭代过程中仍会调用较多的线搜索。上述充分小于是指下一个能量值与比较对象之间的差值满足预设条件，一般地，可以自适应选取差值比较条件。
本实施例提供了一种简单NLS(SNLS)：当能量值连续上升p次时，将原子位置移回并运行GD；否则，仍然运行交替BB方法。需要说明的是，非单调线搜索策略还可以采用上述MNLS、ANLS及SNLS的变体。
基于MNLS，上述判断是否接受试探位置
Figure BDA0003412809360000125
的步骤，可按照以下方式执行：根据预先储存的最近m个试探位置的能量值
Figure BDA0003412809360000131
判断当前的试探位置
Figure BDA0003412809360000132
的能量值是否充分小于
Figure BDA0003412809360000133
若是，则确定接受该试探位置
Figure BDA0003412809360000134
若否则确定不接受所述试探位置
Figure BDA0003412809360000135
基于ANLS，上述判断是否接受试探位置
Figure BDA0003412809360000136
的步骤，可按照以下方式执行：根据预先储存的试探位置的能量值，判断当前的试探位置
Figure BDA0003412809360000137
的能量值是否小于预先存储的所有能量值的凸组合；若是，则确定接受该试探位置
Figure BDA0003412809360000138
若否则确定不接受所述试探位置
Figure BDA0003412809360000139
上述MNLS及ANLS均可以根据试探位置的体系能量是否充分下降判断是否能够接受试探位置，若充分下降则能够接受，否则不能接受。需要说明的是，在初始输入参数时，根据MNLS的要求需要设置储存过往能量值数目m>0。
基于SNLS，上述判断是否接受试探位置
Figure BDA00034128093600001310
的步骤，可按照以下方式执行：判断试探位置能量值的连续上升次数是否大于或等于预设次数阈值；若否则确定接受该试探位置
Figure BDA00034128093600001311
若是则确定不接受所述试探位置
Figure BDA00034128093600001312
在初始输入参数时，根据SNLS的要求需要设置预设次数阈值p>0。
在试探位置的能量值连续上升的情况下，可以判定该试探位置不能被接受。该预设次数阈值可以根据原子体系不同灵活确定，与原子体系对能量上升的容忍程度相关，本实施例对此不作限定。示例性地，该预设次数阈值为2，则试探位置
Figure BDA00034128093600001313
体系能量
Figure BDA00034128093600001314
满足
Figure BDA00034128093600001315
即满足等于预设次数阈值条件。
上述SNLS，可以将上述MNLS及ANLS纳入框架之中，即在判断条件中也加入上述能量充分下降的判定条件。
通过NLS策略判定试探位置是否被接受，在被接受的条件下即可不调用线搜索算法(图2中虚线框中步骤)，极大地减少线搜索带来的昂贵开销。
步骤7，执行线搜索算法，在满足线搜索终止条件的情况下执行步骤8，否则再次执行步骤7。
步骤8，更新原子位置
Figure BDA0003412809360000141
并令k:＝k+1；执行步骤2。
基于SNLS，该执行线搜索算法的步骤，包括：
首先，将当前的原子位置替换为能量值连续上升之前的上一试探位置。该上一试探位置在经过上述连续上升次数的更新后可以得到当前的试探位置。接举前例，试探位置
Figure BDA0003412809360000142
体系能量
Figure BDA0003412809360000143
满足
Figure BDA0003412809360000144
则将当前的原子位置Rk替换为试探位置Rk-1，如前述试探位置Rk-1经过2次循环得到试探位置
Figure BDA0003412809360000145
适应性地，令Fk:＝Fk-1,Dk:＝Fk,Ek:＝Ek-1。
其次，计算试探步长及计算试探位置，以及判断是否满足线搜索终止条件。若是，则执行上述步骤8；若否，则再次执行步骤7。
基于MNLS或ANLS，该执行线搜索算法的步骤，包括：
首先，计算试探步长及计算试探位置，以及判断是否满足线搜索终止条件；若是，则执行步骤8；若否，则再次执行步骤7。
可选地，上述线搜索算法可以是Brent算法、插值方法、回溯方法及它们的变体。
本发明实施例提供的原子结构弛豫的非单调线搜索方法，使用BB步长作为试探步长，从而能够充分利用其优越于最速下降法的局部收敛速度，以及使用非单调线搜索策略来判定试探位置是否被接受，减少调用线搜索算法的次数，从而极大地减少线搜索带来的昂贵开销。相对于现有基于力的原子弛豫方法，本实施例提供的非单调线搜索方法属于一阶方法，但是基于BB步长提高了局部收敛速度；相对于预条件对CG具有潜在的加速效果，其不需要克服设计普适预条件子的困难；相对于二阶方法，其空间复杂度或计算复杂度也较低。
示例性地，为说明本实施例中基于SNLS的交替BB方法(称为NBB)的全局收敛性和相较于CG加速效果的普遍性，在55个算例(包括有机分子、金属合金、二维材料、半导体等，其中既有平衡位置附近的，也有远离平衡位置的)上进行了测试。在展示结果前，对于s∈{CG,NBB}，定义#KSe,s为算法s求解第e个算例时求解Kohn-Sham方程的次数，cpue,s为算法s求解第e个算例时所用的cpu运行时间，
Figure BDA0003412809360000151
图3示出了CG与NBB在55个算例上的性能分析示意图，左侧为#KSe,s的示意图，右侧为cpue,s的示意图。在图3中上方接近横直线的为NBB方法，下方倾斜上升线为CG方法。
其中
Figure BDA0003412809360000152
Figure BDA0003412809360000153
类似定义。
从图3可知，NBB在所有的算例上均收敛，这体现了NBB的全局收敛性。若以NBB与CG的#KS、cpu比值作为加速比，则NBB在约80％的算例上加速比为1.5，在约30％的算例上加速比为2.0，在约20％的算例上加速比为2.5。这表明NBB带来的加速效果是普遍的。
本发明实施例还提供了一种原子结构弛豫的非单调线搜索装置，图4示出了本发明实施例提供的一种原子结构弛豫的非单调线搜索装置的结构示意图，该装置包括：
初始输入模块401，用于输入初始原子位置R0，并令k:＝0；
收敛判断模块402，用于判断原子受力Fk是否满足收敛准则。若是，则终止；若否，则触发方向更新模块运行；
方向更新模块403，用于更新搜索方向Dk＝Fk；
步长计算模块404，用于计算Barzilai-Borwein步长；
试探位置计算模块405，用于根据所述Barzilai-Borwein步长计算试探位置
Figure BDA0003412809360000154
非单调线搜索判断模块406，用于根据非单调线搜索策略判断是否接受所述试探位置
Figure BDA0003412809360000155
若接受，则触发位置更新模块运行；若不接受，则触发线搜索执行模块运行；
线搜索执行模块407，用于执行线搜索算法，在满足线搜索终止条件的情况下触发位置更新模块运行；
位置更新模块408，用于更新原子位置
Figure BDA0003412809360000161
并令k:＝k+1；触发所述收敛判断模块运行。
本发明实施例提供的原子结构弛豫的非单调线搜索装置，使用BB步长作为试探步长，从而能够充分利用其优越于SD的局部收敛速度，以及使用非单调线搜索策略来判定试探位置是否被接受，减少调用线搜索算法的次数，从而极大地减少线搜索带来的昂贵开销。相对于现有基于力的原子弛豫方法，本发明实施例提供的非单调线搜索方法属于一阶方法，但是基于BB步长提高了局部收敛速度；相对于预条件对CG具有潜在的加速效果，其不需要克服设计普适预条件子的困难；相对于二阶方法，其空间复杂度或计算复杂度也较低。
可选地，基于以下公式计算得到两个用于加速梯度法的Barzilai-Borwein步长：
Figure BDA0003412809360000162
其中，S为前后原子位置之差，Y为后前原子受力之差；
或者，基于以下公式确定Barzilai-Borwein步长：
Figure BDA0003412809360000163
可选地，所述非单调线搜索判断模块，具体用于：判断试探位置能量值的连续上升次数是否大于或等于预设次数阈值；若否则确定接受所述试探位置
Figure BDA0003412809360000164
若是则确定不接受所述试探位置
Figure BDA0003412809360000165
或者，根据预先储存的最近m个试探位置的能量值
Figure BDA0003412809360000166
判断当前的所述试探位置
Figure BDA0003412809360000167
的能量值是否充分小于
Figure BDA0003412809360000168
若是则确定接受所述试探位置
Figure BDA0003412809360000169
若否则确定不接受所述试探位置
Figure BDA00034128093600001610
或者，根据预先储存的试探位置的能量值，判断当前的所述试探位置
Figure BDA00034128093600001611
的能量值是否小于预先存储的所有能量值的凸组合，若是则确定接受所述试探位置
Figure BDA0003412809360000171
若否则确定不接受所述试探位置
Figure BDA0003412809360000172
可选地，所述线搜索执行模块，具体用于：计算试探步长及计算试探位置，以及判断是否满足线搜索终止条件；若是，则触发位置更新模块运行；若否，则再次调用所述线搜索执行模块。
可选地，所述线搜索执行模块，还用于若非单调线搜索策略采用判断试探位置能量值的连续上升次数大于或等于所述预设次数阈值，则在调用线搜索执行模块之前，将当前的原子位置替换为能量值连续上升之前的上一试探位置。
本领域技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程度来指令控制装置来完成，所述的程序可存储于一计算机可读取的存储介质中，所述程序在执行时可包括如上述各方法实施例的流程，其中所述的存储介质可为存储器、磁盘、光盘等。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。
本说明书中所描述的以上内容仅仅是对本发明所作的举例说明。本发明不受上述实施例的限制，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (10)

1.一种原子结构弛豫的非单调线搜索方法，其特征在于，包括以下步骤：
步骤1，输入初始原子位置R0，并令k:＝0；
步骤2：判断原子受力Fk是否满足收敛准则。若是，则终止；若否，则执行步骤3；
步骤3，更新搜索方向Dk＝Fk；
步骤4，计算Barzilai-Borwein步长；
步骤5，根据所述Barzilai-Borwein步长计算试探位置
Figure FDA0003412809350000011
步骤6，根据非单调线搜索策略判断是否接受所述试探位置
Figure FDA0003412809350000012
若接受，则执行步骤8；若不接受，则执行步骤7；
步骤7，执行线搜索算法，在满足线搜索终止条件的情况下执行步骤8；
步骤8，更新原子位置
Figure FDA0003412809350000013
并令k:＝k+1；执行所述步骤2。
2.根据权利要求1所述的方法，其特征在于，基于以下公式计算得到两个用于加速梯度法的Barzilai-Borwein步长：
Figure FDA0003412809350000014
其中，S为前后原子位置之差，Y为后前原子受力之差；
或者，基于以下公式确定Barzilai-Borwein步长：
Figure FDA0003412809350000015
3.根据权利要求1或2所述的方法，其特征在于，所述根据非单调线搜索策略判断是否接受所述试探位置
Figure FDA0003412809350000016
包括：
判断试探位置能量值的连续上升次数是否大于或等于预设次数阈值；若否则确定接受所述试探位置
Figure FDA0003412809350000017
若是则确定不接受所述试探位置
Figure FDA0003412809350000018
或者，
根据预先储存的最近m个试探位置的能量值
Figure FDA0003412809350000021
判断当前的所述试探位置
Figure FDA0003412809350000022
的能量值是否充分小于
Figure FDA0003412809350000023
若是则确定接受所述试探位置
Figure FDA0003412809350000024
若否则确定不接受所述试探位置
Figure FDA0003412809350000025
或者，
根据预先储存的试探位置的能量值，判断当前的所述试探位置
Figure FDA0003412809350000026
的能量值是否小于预先存储的所有能量值的凸组合，若是则确定接受所述试探位置
Figure FDA0003412809350000027
若否则确定不接受所述试探位置
Figure FDA0003412809350000028
4.根据权利要求1所述的方法，其特征在于，所述执行线搜索算法，包括：
计算试探步长及计算试探位置，以及判断是否满足线搜索终止条件；
若是，则执行所述步骤8；若否，则再次执行所述步骤7。
5.根据权利要求3所述的方法，其特征在于，若非单调线搜索策略采用判断试探位置能量值的连续上升次数大于或等于所述预设次数阈值，则在执行所述执行步骤7之前，所述方法还包括：
将当前的原子位置替换为能量值连续上升之前的上一试探位置。
6.根据权利要求4或5所述的方法，其特征在于，所述线搜索算法包括：Brent算法、插值方法或回溯方法。
7.一种原子结构弛豫的非单调线搜索装置，其特征在于，包括：
初始输入模块，用于输入初始原子位置R0，并令k:＝0；
收敛判断模块，用于判断原子受力Fk是否满足收敛准则。若是，则终止；若否，则触发方向更新模块运行；
方向更新模块，用于更新搜索方向Dk＝Fk；
步长计算模块，用于计算Barzilai-Borwein步长；
试探位置计算模块，用于根据所述Barzilai-Borwein步长计算试探位置
Figure FDA0003412809350000029
非单调线搜索判断模块，用于根据非单调线搜索策略判断是否接受所述试探位置
Figure FDA0003412809350000031
若接受，则触发位置更新模块运行；若不接受，则触发线搜索执行模块运行；
线搜索执行模块，用于执行线搜索算法，在满足线搜索终止条件的情况下触发位置更新模块运行；
位置更新模块，用于更新原子位置
Figure FDA0003412809350000032
并令k:＝k+1；触发所述收敛判断模块运行。
8.根据权利要求7所述的装置，其特征在于，基于以下公式计算得到两个用于加速梯度法的Barzilai-Borwein步长：
Figure FDA0003412809350000033
其中，S为前后原子位置之差，Y为后前原子受力之差；
或者，基于以下公式确定Barzilai-Borwein步长：
Figure FDA0003412809350000034
9.根据权利要求7或8所述的装置，其特征在于，所述非单调线搜索判断模块，具体用于：
判断试探位置能量值的连续上升次数是否大于或等于预设次数阈值；若否则确定接受所述试探位置
Figure FDA0003412809350000035
若是则确定不接受所述试探位置
Figure FDA0003412809350000036
或者，
根据预先储存的最近m个试探位置的能量值
Figure FDA0003412809350000037
判断当前的所述试探位置
Figure FDA0003412809350000038
的能量值是否充分小于
Figure FDA0003412809350000039
若是则确定接受所述试探位置
Figure FDA00034128093500000310
若否则确定不接受所述试探位置
Figure FDA00034128093500000311
或者，
根据预先储存的试探位置的能量值，判断当前的所述试探位置
Figure FDA00034128093500000312
的能量值是否小于预先存储的所有能量值的凸组合，若是则确定接受所述试探位置
Figure FDA0003412809350000041
若否则确定不接受所述试探位置
Figure FDA0003412809350000042
10.根据权利要求9所述的装置，其特征在于，所述线搜索执行模块，具体用于：
计算试探步长及计算试探位置，以及判断是否满足线搜索终止条件；
若是，则触发所述位置更新模块运行；若否，则再次调用所述线搜索执行模块。
CN202111534901.5A 2021-12-15 2021-12-15 原子结构弛豫的非单调线搜索方法及装置 Active CN114220489B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111534901.5A CN114220489B (zh) 2021-12-15 2021-12-15 原子结构弛豫的非单调线搜索方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111534901.5A CN114220489B (zh) 2021-12-15 2021-12-15 原子结构弛豫的非单调线搜索方法及装置

Publications (2)

Publication Number Publication Date
CN114220489A true CN114220489A (zh) 2022-03-22
CN114220489B CN114220489B (zh) 2022-12-02

Family

ID=80702270

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111534901.5A Active CN114220489B (zh) 2021-12-15 2021-12-15 原子结构弛豫的非单调线搜索方法及装置

Country Status (1)

Country Link
CN (1) CN114220489B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115563447A (zh) * 2022-09-30 2023-01-03 北京应用物理与计算数学研究所 固定晶格体积晶体结构弛豫的计算方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012061475A2 (en) * 2010-11-02 2012-05-10 University Of Florida Research Foundation, Inc. Systems and methods for fast magnetic resonance image reconstruction
CN107713995A (zh) * 2017-11-15 2018-02-23 陕西师范大学 一种基于惩罚算法的生物发光断层成像光源重建方法
CN108830807A (zh) * 2018-06-01 2018-11-16 哈尔滨工业大学 一种基于mems陀螺辅助的星敏感器图像解运动模糊方法
CN112084655A (zh) * 2020-09-08 2020-12-15 南京众诚土地规划设计咨询有限公司 一种基于非单调线搜索的探地雷达参数反演方法
CN112968535A (zh) * 2021-03-30 2021-06-15 西安理工大学 一种用于全方向无线输电系统的负载线圈位置检测方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012061475A2 (en) * 2010-11-02 2012-05-10 University Of Florida Research Foundation, Inc. Systems and methods for fast magnetic resonance image reconstruction
CN107713995A (zh) * 2017-11-15 2018-02-23 陕西师范大学 一种基于惩罚算法的生物发光断层成像光源重建方法
CN108830807A (zh) * 2018-06-01 2018-11-16 哈尔滨工业大学 一种基于mems陀螺辅助的星敏感器图像解运动模糊方法
CN112084655A (zh) * 2020-09-08 2020-12-15 南京众诚土地规划设计咨询有限公司 一种基于非单调线搜索的探地雷达参数反演方法
CN112968535A (zh) * 2021-03-30 2021-06-15 西安理工大学 一种用于全方向无线输电系统的负载线圈位置检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
PENG-WEI HOU ET AL.: ""Influence of helium on the evolution of irradiation-induced defects in tungsten: An object kinetic Monte Carlo simulation"", 《CHIN. PHYS. B》 *
王丽芳: ""电阻和内耗法探索Cu-Zr-Al-Ag系非晶合金热稳定性和晶化/弛豫行为与玻璃形成能力"", 《中国博士学位论文全文数据库工程科技Ⅰ辑》 *
薛伟等: ""改进的OWL-QN方法解稀疏logistic回归问题"", 《中国科学》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115563447A (zh) * 2022-09-30 2023-01-03 北京应用物理与计算数学研究所 固定晶格体积晶体结构弛豫的计算方法及装置

Also Published As

Publication number Publication date
CN114220489B (zh) 2022-12-02

Similar Documents

Publication Publication Date Title
Larson et al. Derivative-free optimization methods
Denzel et al. Gaussian process regression for geometry optimization
Fritzen et al. On-the-fly adaptivity for nonlinear twoscale simulations using artificial neural networks and reduced order modeling
Griffith et al. An adaptive, formally second order accurate version of the immersed boundary method
Li et al. An accumulative error based adaptive design of experiments for offline metamodeling
Kulunchakov et al. Estimate sequences for stochastic composite optimization: Variance reduction, acceleration, and robustness to noise
Wang et al. A projection method for a system of nonlinear monotone equations with convex constraints
Yu et al. Dual coordinate descent methods for logistic regression and maximum entropy models
Lee et al. A parallel implementation of the simplex function minimization routine
US20240168934A1 (en) Similarity Analysis Using Enhanced MinHash
Liu Piecewise length scale control for topology optimization with an irregular design domain
CN114220489B (zh) 原子结构弛豫的非单调线搜索方法及装置
Branislav et al. A survey of gradient methods for solving nonlinear optimization
Zhang et al. Splash: User-friendly programming interface for parallelizing stochastic algorithms
Carstensen et al. A posteriori finite element error control for the p-Laplace problem
Jardow et al. A new hybrid conjugate gradient algorithm for unconstrained optimization with inexact line search
Chu et al. An alternating rank-k nonnegative least squares framework (ARkNLS) for nonnegative matrix factorization
Xu et al. On the Convergence of (Stochastic) Gradient Descent with Extrapolation for Non-Convex Minimization.
Zhang et al. Fast k-means clustering with Anderson acceleration
Long et al. A kernel approach for pde discovery and operator learning
KR20190064948A (ko) 준지도 학습에서의 꼭지점 중요도를 고려한 레이블 추론 방법 및 시스템
Heida et al. Stochastic homogenization of rate-dependent models of monotone type in plasticity
Gaviano et al. A local search method for continuous global optimization
Berglund et al. Zeroth-order randomized subspace Newton methods
Misev et al. Steepest descent optimisation of Runge–Kutta coefficients for second order implicit finite volume CFD codes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant