CN101019122A - 在非平衡条件下用于分子的量子化学模拟的方法和计算机系统 - Google Patents

Abstract

本发明涉及根据一个或多个外部参数的值的变化来利用外推分析表达近似自洽解或自洽解的变化的方法和计算机系统，所述自洽解用在至少含有两个探针或电极的系统的模型中，该模型基于包含有效单电子势能函数和/或有效单电子哈密顿函数的自洽确定的电子结构计算。本发明的方法包含如下步骤：利用自洽循环计算为第一外部参数的第一值确定对所选函数的第一自洽解；利用自洽循环计算为第一所选外部参数的第二值确定对所选函数的第二自洽解，所述第一所选外部参数的第二值不同于第一所选外部参数的第一值；和根据至少所确定的第一和第二自洽解以及第一所选外部参数的第一和第二值，利用外推表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化。

Description

在非平衡条件下用于分子的量子化学模拟的方法和计算机系统

技术领域

本发明涉及根据一个或多个外部参数的值的变化来利用外推分析或技术表达近似自洽解(self-consistent solution)或自洽解的变化的方法和系统。该自洽解可以用在至少含有两个探针或电极的系统或纳米级系统模型中，并且该模型可能基于包含有效单电子势能函数和/或有效单电子哈密顿函数(Hamiltonian)的自洽确定的电子结构计算。

背景技术

原子级模拟(atomic scale modeling)领域内的方法的最常见例子是密度函数理论(density Functional Theory，DFT)和哈特里-福克(Hartree-Fock，HF)理论，其中，模拟基于需要自洽确定有效单电子势能函数的电子结构计算。DFT的许多应用是研究当外部参数发生变化时系统如何响应。在这样的研究中，必须针对外部参数的每个值进行自洽计算，这是非常耗时的。一种重要应用是计算纳米级器件的伏安(Current-Voltage，I-U)特性。在Stokbro的《计算材料科学)》(Stokbro，Computational Materials Science 27，151(2003))中给出了这种计算的例子，其中，计算了与金表面耦合的苯基二硫醇(Di-Thiol-Phenyl，DTP)分子的I-U特性。该系统例示在图2中，并且该计算遵循在如图5和6所示的流程图2和3中概述的步骤。由于针对每个电压的自洽循环，所以该计算需要大量计算。

本发明的目的是提供确定由一个或多个外部参数的变化引起的自洽解的变化的有效和合理的精确方法。

发明内容

根据本发明，提供了根据一个或多个外部参数的值的变化来利用外推分析或技术表达近似自洽解或自洽解的变化的方法，所述自洽解可以用在至少含有两个探针或电极的系统的模型中，该模型基于包含有效单电子势能函数和/或有效单电子哈密顿函数的自洽确定的电子结构计算，该方法包含：

利用自洽循环计算为第一外部参数的第一值确定对所选函数的第一自洽解；

利用自洽循环计算为第一所选外部参数的第二值确定对所选函数的第二自洽解，所述第一所选外部参数的第二值不同于第一所选外部参数的第一值；和

根据至少所确定的第一和第二自洽解以及第一所选外部参数的第一和第二值，利用外推表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化。这里，可以利用线性外推表达近似自洽解或自洽解的变化。

根据本发明的实施例，该方法可以进一步包含利用自洽循环计算为第一所选外部参数的第三值确定对所选函数的第三自洽解，所述第一所选外部参数的第三值不同于第一所选外部参数的第一和第二值。这里，可以根据至少所确定的第一、第二和第三自洽解以及第一所选外部参数的第一、第二和第三值，利用外推表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化。这里，最好，利用二阶外推表达近似自洽解或自洽解的变化。

最好，被模拟的系统是纳米级器件或包含纳米级器件的系统。此外，最好，系统的模拟包含将一个或多个外部参数作为输入提供给所述探针或电极。

根据本发明的实施例，被模拟的系统是双探针系统，而外部参数是跨越所述两个探针或电极的偏压U，所述双探针系统被模拟成具有通过相互作用区相互耦合的两个实质上半无限的探针或电极。

此外，在本发明的实施例中，被模拟的系统是带有三个探针或电极的三探针系统，而外部参数是第一所选参数和类型与第一所选参数相同的第二所选参数。这里，被模拟的系统可以是带有三个探针或电极的三探针系统，而外部参数是跨越所述电极中的第一和第二电极的第一偏压U1以及跨越所述电极中的第三和第一电极的第二偏压U2，所述三探针系统被模拟成具有通过相互作用区相互耦合的三个实质上半无限的探针或电极。

当被模拟的系统是三探针系统时，本发明的方法可以进一步包含：

利用自洽循环计算为第二所选外部参数的第一值确定对所选函数的第四自洽解；

利用自洽循环计算为第二所选外部参数的第二值确定对所选函数的第五自洽解，所述第二所选外部参数的第二值不同于第二所选外部参数的第一值；和

其中，根据至少所确定的第一和第二自洽解，以及第一所选外部参数的第一和第二值，并进一步根据至少所确定的第四和第五自洽解，以及第二所选外部参数的第一和第二值，利用外推为第一所选外部参数的所选值和第二所选外部参数的所选值表达所选函数的近似自洽解或自洽解的变化的所述表达。这里，可以利用线性外推表达近似自洽解或自洽解的变化。

为三探针系统提供的上述本发明方法可以进一步包含：利用自洽循环计算为第二所选外部参数的第三值确定对所选函数的第六自洽解，所述第二所选外部参数的第三值不同于第二所选外部参数的第一和第二值；和根据至少所确定的第一、第二和第三自洽解，以及第一所选外部参数的第一、第二和第三值，并进一步根据至少所确定的第四、第五和第六自洽解，以及第二所选外部参数的第一、第二和第三值，利用外推为第一所选外部参数的所选值和第二所选外部参数的所选值表达所选函数的近似自洽解或自洽解的变化的所述表达。这里，可以利用二阶外推表达近似自洽解或自洽解的变化。

对于为三探针系统提供的本发明方法，第二所选外部参数的第一值可以等于第一所选外部参数的第一值。

根据本发明，最好，从由如下表示的函数中选择所选函数：有效单电子势能函数、有效单电子哈密顿函数和电子密度。这里，最好，所选函数是有效单电子势能函数或有效单电子哈密顿函数，而自洽循环计算基于密度函数理论DFT或哈特里-福克理论HF。

根据本发明的实施例，自洽循环计算可以基于包括如下步骤的循环计算：

a)为系统模型的所选区选择电子密度的值；

b)为所选电子密度和外部参数的所选值确定有效单电子势能函数；

c)计算与确定的有效单电子势能函数相对应的电子密度值；

d)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

e)将有效单电子势能函数的解定义成有效单电子势能函数的自洽解；而如果不相等，那么

f)选择新的电子密度值，并重复步骤b)-f)，直到电子密度的所选值和计算值在所述给定数值精度内相等。这里，可以针对系统的探针或电极区确定有效单电子势能函数的自洽解。

对于针对系统的探针或电极区确定有效单电子势能函数的自洽解的实施例，最好，根据有效单电子势能函数的相应所确定的自洽解，为每个探针或电极区构建或确定格林函数(Green’s function)。

在本发明的方法实施例中，所选函数是系统相互作用区的有效单电子哈密顿函数，而系统相互作用区的有效单电子哈密顿函数的第二自洽解的确定包含在第一所选外部参数的给定值上计算相互作用区的有效单电子势能函数的相应自洽解的步骤。这里，有效单电子哈密顿函数的第二自洽解的确定可以基于包括如下步骤的循环计算：

aa)为系统的相互用作区选择电子密度的值；

bb)为所选外部参数的给定值确定所选电子密度的有效单电子势能函数；

cc)根据在步骤bb)中确定的有效单电子势能函数确定相互作用区的有效单电子哈密顿函数的解；

dd)根据在步骤cc)中确定的有效单电子哈密顿函数的解确定相互作用区的格林函数的解；

ee)计算与相互作用区的确定格林函数相对应的电子密度值；

ff)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

gg)将有效单电子哈密顿函数的解定义成有效单电子哈密顿函数的自洽解；而如果不相等，那么

hh)选择新的电子密度值，并重复步骤bb)-hh)，直到电子密度的所选值和计算值在所述给定数值精度内相等。

根据本发明的实施例，所选函数可以是用哈密顿矩阵表示的有效单电子哈密顿函数，所述矩阵的每个元素是具有利用相应外推表达式表达的近似自洽解或自洽解的变化的函数。

本发明的方法还涵盖了所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压的实施例，其中，分别针对外部偏压的所选第一和第二值，为有效单电子哈密顿函数确定第一和第二自洽解，从而，当外部偏压发生变化时，获取外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述方法进一步包含：利用表达有效单电子哈密顿函数的近似自洽解或自洽解的变化的所获得的外推表达式，针对施加偏压的许多不同值确定系统的两个探针之间的电流。这里，获得的外推表达式可以是线性表达式。可以针对外部偏压的给定范围和针对外部偏压的给定电压步长确定电流，和可以利用如下循环确定电流：

aaa)在外部偏压的给定范围内确定最低电压的电流；

bbb)按给定电压步长升高偏压；

ccc)确定新升偏压的电流；和

ddd)重复步骤bbb)和ccc)，直到新升偏压大于偏压的给定范围的最高电压。

此外，在本发明的实施例中，被模拟的系统是双探针系统，并且所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压，所述方法包含：

将外部偏压的所确定的电压范围划分成至少第一和第二电压范围；

为第一和第二电压范围确定与所述电压范围的最大和最小值相对应的有效单电子哈密顿函数的最大和最小自洽解；

当外部偏压发生变化时，获取第一外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述第一外推表达式基于第一电压范围的确定最大和最小自洽解以及第一电压范围的最大和最小电压值；

当外部偏压发生变化时，获取第二外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述第二外推表达式基于第二电压范围的确定最大和最小自洽解以及第二电压范围的最大和最小电压值；

利用获得的第一外推表达式，在由第一电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流；和

利用获得的第二外推表达式，在由第二电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。这里，获得的第一和第二外推表达式可以分别是第一和第二线性表达式。此外，在本发明的方法实施例中，将确定的电压范围划分成至少三个电压范围，并且该方法进一步包含：

为第三电压范围确定与第三电压范围的最大和最小值相对应的有效单电子哈密顿函数的最大和最小自洽解；

当外部偏压发生变化时，获取第三外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述第三外推表达式基于第三电压范围的确定最大和最小自洽解以及第三电压范围的最大和最小电压值；和

利用获得的第三线性外推表达式，在由第三电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。此外，这里，获得的第三推表达式可以是第三线性外推表达式。

本发明的方法还涵盖了被模拟的系统是双探针系统，并且所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压的实施例，其中，分别针对外部偏压的所选第一和第二值，为有效单电子哈密顿函数确定第一和第二自洽解，所述第二值高于偏压的所选第一值，从而，当外部偏压发生变化时，获取第一外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述方法进一步包含：

aaaa)为外部偏压选择具有最大值和最小值的电压范围，以便在所述范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流；

bbbb)利用自洽循环计算，针对外部偏压的所选最大值确定有效单电子哈密顿函数的最大自洽解；

cccc)根据相应所确定的最大自洽解，针对偏压的最大值确定系统的两个探针之间的电流；

dddd)根据获得的第一外推表达式，针对偏压的所选最大值确定系统的两个探针之间的电流；

eeee)比较在步骤cccc)和dddd)中确定的电流，如果它们在给定数值精度内相等，那么

ffff)当外部偏压发生变化时，将外推表达式用作有效单电子哈密顿函数的近似自洽解，在由所选第一电压值和最大电压值给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。这里，获得的第一外推表达式可以是第一线性外推表达式，和当外部偏压发生变化时，可以将线性外推表达式用在步骤ffff)中表达有效单电子哈密顿函数的近似自洽解。在优选实施例中，获取最大外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述最大外推表达式基于所确定的第一和最大自洽解以及第一偏压和该偏压的最大值，并且其中当在步骤ffff)中确定电流时，使用所述最大外推表达式。最大外推表达式可以是最大线性外推表达式。此外，最好，当在步骤eeee)中确定在步骤cccc)和dddd)中确定的电流值在给定数值精度内不相等时，该方法进一步包含：

gggg)在第一值和前最大值之间选择外部偏压的新最大值；

hhhh)重复步骤bbbb)到hhhh)，直到在步骤cccc)和dddd)中确定的电流值在所述给定数值精度内相等。根据本发明的实施例，该方法进一步包含：

iiii)利用自洽循环计算，针对外部偏压的所选最小值计算有效单电子哈密顿函数的最小自洽解；

jjjj)根据相应所确定的最小自洽解，针对偏压的最小值确定系统的两个探针之间的电流；

kkkk)根据获得的第一外推表达式，针对偏压的所选最小值确定系统的两个探针之间的电流；

llll)比较在步骤jjjj)和kkkk)中确定的电流，如果它们在给定数值精度内相等，那么

mmmm)当外部偏压发生变化时，将外推表达式用作有效单电子哈密顿函数的近似自洽解，在由所选第一电压值和最小电压值给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。这里，当外部偏压发生变化时，可以将线性外推表达式用在步骤mmmm)中表达有效单电子哈密顿函数的近似自洽解。此外，这里，在优选实施例中，获取最小外推表达式，作为有效单电子哈密顿函数的近似自洽解，其中，最小外推表达式基于所确定的第一和最小自洽解以及第一偏压和该偏压的最小值，并且其中当在步骤mmmm)中确定电流时，使用最小外推表达式。这里，最小外推表达式可以是最小线性外推表达式。此外，这里，最好，当在步骤llll)中在步骤jjjj)和kkkk)中确定的电流值在给定数值精度内不相等时，该方法进一步包含：

nnnn)在第一值和前最小值之间选择外部偏压的新最小值；和

oooo)重复步骤iiii)到oooo)，直到在步骤jjjj)和kkkk)中确定的电流值在所述给定数值精度内相等。

根据本发明，还提供了根据一个或多个外部参数的值的变化来利用外推分析表达近似自洽解或自洽解的变化的计算机系统，所述自洽解可以用在至少含有两个探针或电极的纳米级系统的模型中，该模型基于包含有效单电子势能函数和/或有效单电子哈密顿函数的自洽确定的电子结构计算，所述计算机系统包含：

利用自洽循环计算为第一外部参数的第一值确定对所选函数的第一自洽解的装置；

利用自洽循环计算为第一所选外部参数的第二值确定对所选函数的第二自洽解的装置，所述第一所选外部参数的第二值不同于第一所选外部参数的第一值；和

根据至少所确定的第一和第二自洽解以及第一所选外部参数的第一和第二值，利用外推表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化的装置。这里，表达近似自洽解或自洽解的变化的装置适用于利用线性外推表达这样的解。

根据本发明的实施例，该计算机系统可以进一步包含：利用自洽循环计算为第一所选外部参数的第三值确定对所选函数的第三自洽解的装置，所述第一所选外部参数的第三值不同于第一所选外部参数的第一和第二值。这里，表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化的装置适用于根据至少所确定的第一、第二和第三自洽解以及第一所选外部参数的第一、第二和第三值，利用外推表达这样的解。这里，最好，表达近似自洽解或自洽解的变化的装置适用于利用二阶外推表达这样的解。

对于本发明的计算机系统，在实施例中，纳米级系统是双探针系统和外部参数是跨越所述两个探针或电极的偏压U，所述双探针系统被模拟成具有通过相互作用区相互耦合的两个实质上半无限的探针或电极。

本发明的计算机系统还涵盖了纳米级系统是带有三个探针或电极的三探针系统和外部参数是第一所选参数和类型与第一所选参数相同的第二所选参数的实施例。这里，最好，纳米级系统是带有三个探针或电极的三探针系统，而外部参数是跨越所述电极中的第一和第二电极的第一偏压U1以及跨越所述电极中的第三和第一电极的第二偏压U2，所述三探针系统被模拟成具有通过相互作用区相互耦合的三个实质上半无限的电极。

此外，这里，当被模拟的纳米级系统是三探针系统时，本发明的计算机系统可以进一步包含：

利用自洽循环计算为第二所选外部参数的第一值确定对所选函数的第四自洽解的装置；

利用自洽循环计算为第二所选外部参数的第二值确定对所选函数的第五自洽解的装置，所述第二所选外部参数的第二值不同于第二所选外部参数的第一值；和

其中，表达所选函数的近似自洽解或自洽解的变化的所述装置适用于根据所确定的第一和第二自洽解，以及第一所选外部参数的第一和第二值，并进一步根据所确定的第四和第五自洽解，以及第二所选外部参数的第一和第二值，利用外推表达第一所选外部参数的所选值和第二所选外部参数的所选值的近似自洽解。这里，表达近似自洽解或自洽解的变化的装置适用于利用线性外推表达这样的解。

将三探针系统模拟的上述计算机系统可以进一步包含：

利用自洽循环计算为第二所选外部参数的第三值确定对所选函数的第六自洽解的装置，所述第二所选外部参数的第三值不同于第二所选外部参数的第一和第二值。这里，表达所选函数的近似自洽解或自洽解的变化的装置适用于根据至少所确定的第一、第二和第三自洽解，以及第一所选外部参数的第一、第二和第三值，并进一步根据至少所确定的第四、第五和第六自洽解，以及第二所选外部参数的第一、第二和第三值，利用外推表达第一所选外部参数的所选值和第二所选外部参数的所选值的近似自洽解。这里，表达近似自洽解或自洽解的变化的装置适用于利用二阶外推表达这样的解。

对于为三探针系统提供的本发明系统，第二所选外部参数的第一值可以等于第一所选外部参数的第一值。

此外，对于本发明的计算机系统，最好，从由如下表示的函数中选择所选函数：有效单电子势能函数、有效单电子哈密顿函数和电子密度。这里，最好，所选函数是有效单电子势能函数或有效单电子哈密顿函数，而自洽循环计算基于密度函数理论DFT或哈特里-福克理论HF。

根据本发明的实施例，该计算机系统可以进一步包含根据包括如下步骤的循环计算进行自洽循环计算的装置：

a)为纳米级系统模型的所选区选择电子密度的值；

b)为所选电子密度和外部参数的所选值确定有效单电子势能函数；

c)计算与确定的有效单电子势能函数相对应的电子密度；

d)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

e)将有效单电子势能函数的解定义成有效单电子势能函数的自洽解；而如果不相等，那么

f)选择新的电子密度值，并重复步骤b)-f)，直到电子密度的所选值和计算值在所述给定数值精度内相等。这里，进行自洽循环计算的装置适用于针对系统的探针或电极区确定有效单电子势能函数的自洽解。

此外，对于进行自洽循环计算的装置适用于针对系统的探针或电极区确定有效单电子势能函数的自洽解的实施例，最好，该计算机系统包含根据有效单电子势能函数的相应所确定的自洽解，为每个探针或电极区确定格林函数的装置。

此外，对于本发明的计算机系统，在实施例中，所选函数是系统相互作用区的有效单电子哈密顿函数，和确定系统相互作用区的有效单电子哈密顿函数的第二自洽解的装置适用于通过包括在第一所选外部参数的给定值上计算相互作用区的有效单电子势能函数的相应自洽解的步骤进行所述确定。这里，确定有效单电子哈密顿函数的第二自洽解的装置适用于根据包括如下步骤的循环计算进行所述确定：

aa)为系统的相互用作区选择电子密度的值；

bb)为所选外部参数的给定值确定所选电子密度的有效单电子势能函数；

cc)根据在步骤bb)中确定的有效单电子势能函数确定相互作用区的有效单电子哈密顿函数的解；

dd)根据在步骤cc)中确定的有效单电子哈密顿函数的解确定相互作用区的格林函数的解；

ee)计算与相互作用区的确定格林函数相对应的电子密度值；

ff)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

gg)将有效单电子哈密顿函数的解定义成有效单电子哈密顿函数的自洽解；而如果不相等，那么

hh)选择新的电子密度值，并重复步骤bb)-hh)，直到电子密度的所选值和计算值在所述给定数值精度内相等。

此外，本发明的计算机系统还涵盖了所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压的实施例，其中，确定第一和第二自洽解的装置适用于分别针对外部偏压的所选第一和第二值，为有效单电子哈密顿函数进行所述确定，并且其中利用外推分析表达近似自洽解的装置适用于当外部偏压发生变化时，获取外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述计算机系统进一步包含：利用表达有效单电子哈密顿函数的近似自洽解或自洽解的变化的所获得的外推表达式，针对施加偏压的许多不同值确定系统的两个探针之间的电流的装置。这里，获得的外推表达式可以是线性外推表达式。确定电流的装置适用于针对外部偏压的给定范围和针对外部偏压的给定电压步长确定电流，和确定电流的装置适用于利用如下循环进行所述确定：

aaa)在外部偏压的给定范围内确定最低电压的电流；

bbb)按给定电压步长升高偏压；

ccc)确定新升偏压的电流；和

ddd)重复步骤bbb)和ccc)，直到新升偏压大于偏压的给定范围的最高电压。

此外，在本发明的计算机系统的实施例中，被模拟的系统是双探针系统，并且所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压，并且其中该计算机系统进一步包含：

将外部偏压的所确定的电压范围划分成至少第一和第二电压范围的装置；

为第一和第二电压范围确定与所述电压范围的最大和最小值相对应的有效单电子哈密顿函数的最大和最小自洽解的装置；

当外部偏压发生变化时、获取第一外推表达式作为有效单电子哈密顿函数的近似自洽解的装置，所述第一外推表达式基于第一电压范围的确定最大和最小自洽解以及第一电压范围的最大和最小电压值；

当外部偏压发生变化时、获取第二外推表达式作为有效单电子哈密顿函数的近似自洽解的装置，所述第二外推表达式基于第二电压范围的确定最大和最小自洽解以及第二电压范围的最大和最小电压值；

利用获得的第一外推表达式、在由第一电压范围的最小和最大电压给出的电压范围内针对施加偏压的许多不同值确定系统的两个探针之间的电流的装置；和

利用获得的第二外推表达式、在由第二电压范围的最小和最大电压给出的电压范围内针对施加偏压的许多不同值确定系统的两个探针之间的电流的装置。这里，获得的第一和第二外推表达式可以分别是第一和第二线性表达式。

通过结合附图对本发明优选实施例进行详细描述，本发明的其它目的、特征和优点将更加清楚。

附图说明

图1是例示密度函数理论(DFT)的自洽循环中的计算步骤的流程图(流程图1)；

图2例示了与两个金(111)表面耦合的苯二硫醇分子，其中，金表面与外部电压源耦合，并且电极具有不同的化学势μ_L和μ_R。

图3示出了与金表面耦合的碳纳米管的自洽电子密度，其中，当在相互作用区以外时，电子密度由电极的体积密度(bulk density)给出；

图4a示出了可以用于积分格林函数G_I(z)的等效实轴(R)和复轮廓线(complex contour)(C)；

图4b示出了谱密度(

ImG_I(z))沿着图4a的轮廓线C(虚线)和沿着图4a的实轴R(实线)的变化；

图5是示出利用格林函数方法计算具有所施加的电压U的双探针系统的自洽有效势能函数所需的步骤的流程图(流程图2)，其中可以从自洽有效单电子势能函数中计算电流I；

图6是示出自洽计算双探针系统的I-U特性所需的步骤的流程图(流程图3)；

图7a示出了针对施加电压的不同值计算的例示在图2中的系统的自洽有效单电子势能函数；

图7b示出了用施加电压重新调整的自洽有效单电子势能函数；

图8是示出利用根据本发明实施例的线性外推表达式计算伏安特性I-U时所涉及的步骤的流程图(流程图4)；

图9是示出如何利用根据本发明实施例的内插公式或线性外推表达式计算I-U特性的流程图(流程图5)；

图10示出了例示在图2中的系统的I(U)特性的计算结果，用“SCF”表示的线段示出了利用自洽计算获得的结果，而用“1阶”表示的线段示出了利用例示在图8中的方案获得的结果，和用“2阶”表示的线段是二阶近似；

图11是例示将根据本发明实施例的自适应网格算法用于计算伏安特性I-U的流程图(流程图6)；

图12是作为图11的流程图6使用的递归流程图的流程图(流程图7)；和

图13是作为图11的流程图6使用的递归流程图的流程图(流程图8)。

具体实施方式

背景理论

原子级模拟的目的是根据系统中的各个原子的描述计算分子和材料的特性。原子由电荷Z的离子实(ion core)和抵消这个电荷的等量电子组成。我们将 和Z_μ用于离子的位置和电荷，其中，μ＝1，...，N标记离子，N是离子数(number)。电子的位置由

i＝1，...，n给出，n是电子数。

通常，将离子当作经典粒子是不错的近似。离子的势能 通过如下方程依赖于电子系统的能量E₀：

$V=({\stackrel{\→}{R}}_1,...,{\stackrel{\→}{R}}_N)=E_0+\frac{1}{2}\underset{\mathrm{\μ},{\mathrm{\μ}}^{\′}=1}{\overset{N}{\mathrm{\Σ}}}\frac{Z_{\mathrm{\μ}}{Z}_{{\mathrm{\μ}}^{\′}}{e}^2}{|R_{\mathrm{\μ}}-{\stackrel{\→}{R}}_{{\mathrm{\μ}}^{\′}}|},$ 方程(1)

其中，e是电子电荷。电子必须被描述成量子粒子，和电子能量计算需要我们求解多体薛定谔(Schrdinger)波动方程：

$\hat{H}\mathrm{\Ψ}({\stackrel{\→}{r}}_1,...,{\stackrel{\→}{r}}_n)=E_0\mathrm{\Ψ}({\stackrel{\→}{r}}_1,...,{\stackrel{\→}{r}}_n),$ 方程(2)

方程(3)

在方程(2)中，

是多体哈密顿函数，Ψ是电子的多体波函数。多体哈密顿函数

上面的“^”表示该量是量子力学算符。方程(3)中的第一项是电子的动能，

，其中，h是普朗克(Planck)常数，m是电子质量，

是关于 的梯度。第二项是静电电子-离子吸引力，最后一项是静电电子-电子排斥力。

最后一项耦合不同的电子，并引起电子之间的关联运动。由于这种复杂性，多体薛定谔方程的精确解只有针对单电子的系统才有可能。因此，需要可以将多体薛定谔方程简化成实际可解模型的近似。许多成功的方法使用有效单电子哈密顿函数来描述电子结构，并通过单电子哈密顿函数中的有效单电子势能函数将电子-电子相互作用包括进来。

计算电子结构的密度函数方法

本发明可以与用有效单电子哈密顿函数描述电子的电子结构方法一起使用。DFT和HF理论是这样方法的例子。在这些方法中，电子被描述成在其它电子建立的有效单电子势中运动的非相互作用粒子。有效单电子势依赖于其它的电子平均位置，并需要自洽地确定。

方程(4)

在方程(4)中，

项描述动能，

是有效单电子势能函数，而

是单电子哈密顿函数。有效单电子势能函数依赖于电子密度n。动能由简单微分算符给出，因此，与密度无关。这意味着，有效单电子势能函数和哈密顿函数具有与密度函数相同的变量(variant)，当我们对确定有效单电子势能函数的自洽变化感兴趣时，这等效于指定哈密顿函数的自洽变化。此外，对于自洽解，在电子密度和有效单电子势能函数之间存在一一对应关系，因此，指定自洽电子密度、哈密顿函数或有效单电子势是等效的。

在DFT中，有效单电子势能函数由下式给出：

V^eff[n]＝V^ion+V^xc[n]+V^H[n]。                         方程(5)

第一项是通过来自离子实的静电势能给出的离子势能函数：

$V^{\mathrm{ion}}\left(\stackrel{\→}{r}\right)=\underset{\mathrm{\μ}=1}{\overset{N}{\mathrm{\Σ}}}\frac{Z_{\mathrm{\μ}}{e}^2}{|\stackrel{\→}{r}-{\stackrel{\→}{R}}_{\mathrm{\μ}}|},$ 方程(6)

因此，与n无关。第二项是交换相关势能函数(exchange-correlation potentialenergy function)：

$V^{\mathrm{xc}}\left(\stackrel{\→}{r}\right)=f(n\left(\stackrel{\→}{r}\right),\stackrel{\→}{\▿}n\left(\stackrel{\→}{r}\right),{\stackrel{\→}{\▿}}^{2}n\left(\stackrel{\→}{r}\right)),$ 方程(7)

它是密度和它的梯度的局部函数(local function)。第三项是哈特里势能函数，它是源自电子密度的静电势能，可以根据泊松(Poisson)方程计算出来：

${\stackrel{\→}{\▿}}^2{V}^H\left(\stackrel{\→}{r}\right)=-4\mathrm{\πen}\left(\stackrel{\→}{r}\right).$ 方程(8)

泊松方程是二阶微分方程，需要边界条件来固定解。对于孤立系统，边界条件是势能函数渐近地变成零；而在周期性系统中，边界条件是势能函数是周期性的。对于这样的边界条件，泊松方程的解是直截了当的，可以从标准数值计算软件包中获得V^H。对于存在外部电压U的系统，我们求解系统的分立部分中的哈特里势。这种情况将在下面作更详细讨论。

因此，我们可以从密度中获得有效单电子势能函数，从而获得哈密顿函数。下一步是从哈密顿函数中计算电子密度。这可以通过求和所有所占据的单电子本征态(occupied one-electron eigenstate)获得：

${\hat{H}}_{1\mathrm{el}}{\mathrm{\ψ}}_{\mathrm{\α}}\left(\stackrel{\→}{r}\right)={\mathrm{\ϵ}}_{\mathrm{\α}}{\mathrm{\ψ}}_{\mathrm{\α}}\left(\stackrel{\→}{r}\right),$ 方程(9)

$n\left(\stackrel{\→}{r}\right)=\underset{\mathrm{\α}\∈\mathrm{occ}}{\mathrm{\Σ}}{\left|{\mathrm{\ψ}}_{\mathrm{\α}}\left(\stackrel{\→}{r}\right)\right|}^2.$ 方程(10)

对于具有单化学势的系统，所占据的本征态是具有低于化学势的能量的态。对于具有所施加的外部电压U的系统，存在两个化学势，情况更复杂。这种情况将在下面加以描述。

图1中的流程图例示了求解方程所需的自洽循环。系统通过原子的位置R_μ(离子坐标)以及像施加电压U、温度T和压强P那样的外部参数来定义(102)。首先，我们任意猜测系统的电子密度(104)。根据该密度，我们利用方程(5)构建有效单电子势能函数(106)。有效单电子势能函数通过方程(4)定义哈密顿函数(108)。如方程(9)和(10)所示，根据哈密顿函数，我们可以通过求和所有所占据的单电子本征态计算系统的电子密度。如果新密度等于(在指定数值精度内)用于构建有效单电子势能函数的密度(112)，我们就获取自洽解(114)，然后停止(116)。如果输入和输出电子密度不同，我们根据前面计算的电子密度作出新的改进猜测。在最简单的形式中，以混合参数β线性混合两个电子密度获取新猜测(110)。

将DFT应用于封闭周期性系统

我们将首先示出如何针对周期性封闭系统最一般地求解方程(9)。封闭系统是原子数有限的系统。周期性系统是无限个原子以周期性结构排列的系统。对于这些系统，通常将方程(9)转换成可以用标准线性代数软件包求解的矩阵本征值问题。这种转换是通过将波函数ψ_α写成基函数的线性组合

获得的。对于基函数_i，存在许多不同选择，一些最常用的是平面波或形状像原子波函数的基于原子的函数。利用基函数，将方程(9)转换成：

$\underset{j}{\mathrm{\Σ}}{\stackrel{\‾}{H}}_{\mathrm{ij}}{a}_j^{\mathrm{\α}}={\mathrm{\ϵ}}_{\mathrm{\α}}\underset{j}{\mathrm{\Σ}}{\stackrel{\‾}{S}}_{\mathrm{ij}}{a}_j^{\mathrm{\α}},$ 方程(11)

方程(12)

方程(13)

方程(14)

符号

表示哈密顿矩阵，

表示重叠矩阵。字母上面的“-”表示这些量是矩阵。

对于分子系统，哈密顿矩阵是有限的，可以用标准线性代数软件包对角化。对于周期性结果，只需模拟当重复时生成整个结构的部分系统。因此，哈密顿矩阵也是有限的，解也是直截了当的。

将DFT应用于存在施加电压的开放系统

本发明的应用领域是两个(或更多个)半无限电极与纳米级相互作用区耦合的系统。我们将这样的系统称为双探针系统。纳米级相互作用区可以与电极交换粒子，因此，双探针系统是开放量子力学系统。左右电极是具有确定化学势μ_L和μ_R的电子库。化学势之差：

μ_L-μ_R＝eU                                      方程(15)

定义施加给系统的偏压U。对于开放系统，哈密顿矩阵是无限的，和不能应用获取单电子本征态的方程(11)中的简单对角化技术。取而代之，我们利用在如下部分中描述的非平衡格林函数形式确定电子密度。双探针系统的例子例示在图2和3中。图2中的系统由与苯基二硫醇分子耦合的两个半无限金电极组成。相互作用层22由分子和电极的前两层组成。区域21和23示出了左右电极区。区域24和26示出了电极内单电子能级的占据(occupation ofthe one-electron level)；由于施加的电压，右电极26的化学势高于左电极24的化学势。

图3示出了与半无限金线耦合的半无限碳纳米管。相互作用层32由纳米管顶端和金线的前几层给出。左电极31由半无限金线组成，而右电极33由半无限碳纳米管组成。左电极区34中和右电极区36中的电子密度通过自洽体积计算(bulk calculation)获得。这些密度与自洽计算的相互作用区35的双探针密度无缝匹配。

筛选近似(screening approximation)

第一步是将开放系统转换成可以独立求解的三个子系统。图3a示出了与金线耦合的碳纳米管。金线和碳纳米管都是金属的。由于半无线电极的金属性质，由相互作用区引起的摄动只进入电极几个埃()。这例示在示出电子密度的图3b中。可以看到，当我们离开纳米管向金触点移动几个原子距离时，电子密度是周期性的，像体积电子密度(bulk electron density)。因此，我们可以将电子密度和有效单电子势能函数划分成相互作用区和电极区，其中，电极区中的值类似于电极体积值。这被称为筛选近似。

由于有效单电子势能函数是局部算符，所以也可以将哈密顿算符分离成电极区和相互作用区。因此，如果我们用有限范围的基集展开哈密顿算符，可以将哈密顿矩阵分离成：

$\stackrel{\‾}{H}=\left(\begin{array}{ccc}{\stackrel{\‾}{H}}_{\mathrm{LL}}& {\stackrel{\‾}{H}}_{\mathrm{LI}}& 0\\ {\stackrel{\‾}{H}}_{\mathrm{IL}}& {\stackrel{\‾}{H}}_{\mathrm{II}}& {\stackrel{\‾}{H}}_{\mathrm{IR}}\\ 0& {\stackrel{\‾}{H}}_{\mathrm{RI}}& {\stackrel{\‾}{H}}_{\mathrm{RR}}\end{array}\right),$ 方程(16)

其中，

和 分别表示左电极、相互作用区和右电极的哈密顿矩阵，而

和

是与相互作用区和电极有关的矩阵元。注意，相互作用区的大小使得左电极与右电极之间没有耦合，即，H_LR＝H_RL＝0。

利用格林函数计算电子密度

现在，我们将示出如何在格林函数形式下获取电子密度。为此，我们引入谱密度

和电子密度算符

谱密度是单位能量电子密度，总电子密度通过对所有能量积分谱密度获得：

$\hat{D}\left(\mathrm{\ϵ}\right)=\mathrm{\δ}(\mathrm{\ϵ}-\hat{H}),$ 方程(17)

$\hat{N}=\underset{-\∞}{\overset{\mathrm{\μ}}{\∫}}\hat{D}\left(\mathrm{\ϵ}\right)\mathrm{d\ϵ}.$ 方程(18)

在方程(17)中，δ(x)是狄拉克(Dirac)δ函数。(推迟)格林函数由下式定义：

$\hat{G}\left(\mathrm{\ϵ}\right)={[\mathrm{\ϵ}-\hat{H}+i{\mathrm{\δ}}_{+}]}^{-1},$ 方程(19)

其中，是δ₊无穷小正数，i是复基数(complex base)。格林函数通过下式与谱密度相联系：

$\hat{D}\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\π}}\mathrm{Im}\hat{G}\left(\mathrm{\ϵ}\right),$ 方程(20)

其中，

是

的虚部。当用基函数展开算符时，我们将方程(19)转换成矩阵方程：

$\stackrel{\‾}{G}\left(\mathrm{\ϵ}\right){=[(\mathrm{\ϵ}+{\mathrm{i\δ}}_{+})\stackrel{\‾}{S}-\stackrel{\‾}{H}]}^{-1},$ 方程(21)

我们可以从格林函数中获得谱密度矩阵：

$\stackrel{\‾}{D}\left(\mathrm{\ϵ}\right)=\frac{1}{\mathrm{\π}}\mathrm{Im}\stackrel{\‾}{G}\left(\mathrm{\ϵ}\right),$ 方程(22)

因此，获得电子密度：

方程(23)

$\stackrel{\‾}{N}=\underset{-\∞}{\overset{\mathrm{\μ}}{\∫}}\stackrel{\‾}{D}\left(\mathrm{\ϵ}\right)\mathrm{d\ϵ}.$ 方程(24)

现在，电子密度的计算被简化成方程(21)中的矩阵求逆和方程(24)中的能量积分。但是，我们拥有开放系统，因此，方程(21)中的矩阵是无限的。由于筛选近似，我们只需计算相互作用区中的电子密度，因为在电极区中，我们可以使用体积电子密度。从方程(24)中可以看出，由于我们的基函数是局部化的，我们只需计算相互作用区和电极的几个层的格林函数矩阵。

通过自能(self energy)项将电极区包括进来

在这个部分中，我们将示出如何通过求逆大小相同的矩阵计算相互作用区的格林函数矩阵

为了获得这个结果，我们将摄动理论应用在耦合元 ${\stackrel{\‾}{\tilde{H}}}_{\mathrm{LI}}\left(\mathrm{\ϵ}\right)={\stackrel{\‾}{H}}_{\mathrm{LI}}-\mathrm{\ϵ}{\stackrel{\‾}{S}}_{\mathrm{LI}}$ 和 ${\stackrel{\‾}{\tilde{H}}}_{\mathrm{RI}}\left(\mathrm{\ϵ}\right)={\stackrel{\‾}{H}}_{\mathrm{RI}}-\mathrm{\ϵ}{\stackrel{\‾}{S}}_{\mathrm{RI}}$ 中。通过设置 ${\stackrel{\‾}{\tilde{H}}}_{\mathrm{LI}}={\stackrel{\‾}{\tilde{H}}}_{\mathrm{RI}}=0$ 和利用在这种情况下方程(21)是块对角的特点计算无摄动格林函数G⁰：

${\stackrel{\‾}{G}}_{\mathrm{LL}}^0\left(\mathrm{\ϵ}\right)={[(\mathrm{\ϵ}+i{\mathrm{\δ}}_{+}){\stackrel{\‾}{S}}_{\mathrm{LL}}-{\stackrel{\‾}{H}}_{\mathrm{LL}}]}^{-1},$ 方程(25)

${\stackrel{\‾}{G}}_{\mathrm{II}}^0\left(\mathrm{\ϵ}\right)={[(\mathrm{\ϵ}+i{\mathrm{\δ}}_{+}){\stackrel{\‾}{S}}_{\mathrm{II}}-{\stackrel{\‾}{H}}_{\mathrm{II}}]}^{-1},$ 方程(26)

${\stackrel{\‾}{G}}_{\mathrm{RR}}^0\left(\mathrm{\ϵ}\right)={[(\mathrm{\ϵ}+i{\mathrm{\δ}}_{+}){\stackrel{\‾}{S}}_{\mathrm{RR}}-{\stackrel{\‾}{H}}_{\mathrm{RR}}]}^{-1},$ 方程(27)

将摄动

和

放回去，我们从代森(Dyson)方程中求出格林函数：

${\stackrel{\‾}{G}}_{\mathrm{II}}\left(\mathrm{\ϵ}\right)={\stackrel{\‾}{G}}_{\mathrm{II}}^0\left(\mathrm{\ϵ}\right)+{\stackrel{\‾}{G}}_{\mathrm{II}}^0\left(\mathrm{\ϵ}\right)[{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{II}}^L\left(\mathrm{\ϵ}\right)+{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{II}}^R\left(\mathrm{\ϵ}\right)]{\stackrel{\‾}{G}}_{\mathrm{II}}\left(\mathrm{\ϵ}\right),$ 方程(28)

${\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{II}}^L\left(\mathrm{\ϵ}\right)={\stackrel{\‾}{\tilde{H}}}_{\mathrm{IL}}\left(\mathrm{\ϵ}\right){\stackrel{\‾}{G}}_{\mathrm{LL}}^0\left(\mathrm{\ϵ}\right){\stackrel{\‾}{\tilde{H}}}_{\mathrm{LI}}\left(\mathrm{\ϵ}\right),$ 方程(29)

${\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{II}}^R\left(\mathrm{\ϵ}\right)={\stackrel{\‾}{\tilde{H}}}_{\mathrm{IR}}\left(\mathrm{\ϵ}\right){\stackrel{\‾}{G}}_{\mathrm{RR}}^0\left(\mathrm{\ϵ}\right){\stackrel{\‾}{\tilde{H}}}_{\mathrm{RI}}\left(\mathrm{\ϵ}\right),$ 方程(30)

其中，项

和

被称为电极的自能。重新排列代森方程中的项，我们得出：

${\stackrel{\‾}{G}}_{\mathrm{II}}\left(\mathrm{\ϵ}\right)={[(\mathrm{\ϵ}+{\mathrm{i\δ}}_{+}){\stackrel{\‾}{S}}_{\mathrm{II}}-{\stackrel{\‾}{H}}_{\mathrm{II}}-{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{II}}^L\left(\mathrm{\ϵ}\right)-{\stackrel{\‾}{\mathrm{\Σ}}}_{\mathrm{II}}^R\left(\mathrm{\ϵ}\right)]}^{-1}.$ 方程(31)

电极格林函数的计算

为了确定自能，我们需要计算电极的无摄动格林函数

由于电极的哈密顿函数是半无限的，不能通过简单矩阵求逆获取格林函数。但是，在电极哈密顿函数周期性的情况下，存在非常有效的计算电极格林函数的算法。下面我们将描述这些算法之一。我们将电极哈密顿函数写成周期性块 ${\stackrel{\‾}{H}}_{L_1{L}_1}={\stackrel{\‾}{H}}_{L_2{L}_2}=\·\·\·,$ 其中，每个块的大小使得只有相邻块才有相互作用，即：

方程(32)

每个块的哈密顿函数

和耦合矩阵

通过电极系统的体积计算获得。利用递归，我们为格林函数构建出一系列近似：

${\stackrel{\‾}{G}}_{L_1{L}_1}^{0\left[0\right]}\left(\mathrm{\ϵ}\right)={[(\mathrm{\ϵ}+i{\mathrm{\δ}}_{+}){\stackrel{\‾}{S}}_{L_1{L}_1}-{\stackrel{\‾}{H}}_{L_1{L}_1}]}^{-1},$ 方程(33)

${\stackrel{\‾}{G}}_{L_1{L}_1}^{0\left[1\right]}\left(\mathrm{\ϵ}\right)={[(\mathrm{\ϵ}+i{\mathrm{\δ}}_{+}){\stackrel{\‾}{S}}_{L_1{L}_1}-{\stackrel{\‾}{H}}_{L_1{L}_1}-{\stackrel{\‾}{H}}_{L_1{L}_2}{\stackrel{\‾}{G}}_{L_2{L}_2}^{0\left[0\right]}{\stackrel{\‾}{H}}_{L_2{L}_1}]}^{-1},$ 方程(34)

${\stackrel{\‾}{G}}_{L_1{L}_1}^{0\left[2\right]}\left(\mathrm{\ϵ}\right)={[(\mathrm{\ϵ}+i{\mathrm{\δ}}_{+}){\stackrel{\‾}{S}}_{L_1{L}_1}-{\stackrel{\‾}{H}}_{L_1{L}_1}-{\stackrel{\‾}{H}}_{L_1{L}_2}{\stackrel{\‾}{G}}_{L_2{L}_2}^{0\left[1\right]}{\stackrel{\‾}{H}}_{L_2{L}_1}]}^{-1},$ 方程(35)

                                      方程(36)

在方程(33)、(34)和(35)中，量

是格林函数的n阶近似。误差

随1/n而减小，其中，n是步数。由于这个收敛很慢，通常需要1000多步才能获得这种算法的合理精度。利用Lopez-Sancho在J.Phys.F14，1205(1984)中描述的方法的一种变体，可以在几步内获得格林函数。借助于这种算法的变体，只需要几步就可以计算出电极的格林函数，与计算

所需的资源相比，这个部分所需的计算资源通常可忽略不计。

利用复轮廓线积分谱密度

现在，我们拥有在方程(31)中获取

从而获得以下电子密度矩阵所需的所有成分：

${\stackrel{\‾}{N}}_{\mathrm{ij}}=\frac{1}{\mathrm{\π}}\underset{-\∞}{\overset{\mathrm{\μ}}{\∫}}\mathrm{Im}{\stackrel{\‾}{G}}_{\mathrm{ij}}\left(\mathrm{\ϵ}\right)\mathrm{d\ϵ}.$ 方程(37)

格林函数是沿着实轴迅速变化的函数，对于真实的系统，沿着实轴精确确定积分往往需要5000多个能量点。为了找出更有效的方法，我们利用格林函数是解析函数的特点，我们可以沿着复平面内的轮廓线进行积分。在复平面内，格林函数是非常平滑的。这例示在图4中。在图4中，我们示出了两条等效积分线，即，轮廓线C和实轴R。图4b示出了谱密度沿着C(虚线)和沿着R(实线)的变化。函数沿着R变化非常迅速，为了获得相同的精度，沿着R比沿着C需要多得多的点。通常，使用轮廓线积分使积分点数减少到百分之一。

具有外部偏压的双探针系统的电子密度

到目前为止，我们利用了系统具有单化学势，即，μ_L＝μ_R的特点。但是，当我们施加外部电压U时，两个电极将具有通过方程(15)联系的不同化学势。图2例示了系统的建立。能量轴可以划分成两个区域，即，我们称为平衡区的在两个化学势以下的能量范围以及我们称为非平衡区或电压窗的在两个化学势之间的能量范围。我们将电子密度矩阵划分成两个部分：

${\stackrel{\‾}{N}}_{\mathrm{ij}}={\stackrel{\‾}{N}}_{\mathrm{ij}}^{\mathrm{eq}}+{\stackrel{\‾}{N}}_{\mathrm{ij}}^{\mathrm{neq}},$ 方程(38)

其中，

是能量在平衡区内的电子的电子密度矩阵，而

是能量在非平衡区内的电子的电子密度矩阵。我们可以认为

是由外部电压U引起的附加密度。

可以利用在前面部分中描述的方法计算，因此，

${\stackrel{\‾}{N}}_{\mathrm{ij}}^{\mathrm{eq}}=\frac{1}{\mathrm{\π}}\underset{-\∞}{\overset{{\mathrm{\μ}}_L}{\∫}}\mathrm{Im}{\stackrel{\‾}{G}}_{\mathrm{ij}}\left(\mathrm{\ϵ}\right)\mathrm{d\ϵ},$ 方程(39)

其中，我们假设μ_L＜μ_R。

在非平衡区内，电子只从右库注入。因此，我们需要将谱密度划分成源自左或右电极的电子状态，并且只需加上右电极电子密度。这种电子密度划分是在非平衡格林函数理论下完成的，我们可以写成：

方程(40)

这个方程的理论基础可以在Haug和A.P.Jauho的《半导体的量子输运动力学和光学》(Haug and A.P.Jauho，Quantum kinetics in transport and optics ofsemiconductors，Springer-Verlag，Berlin，1996)或Brandbyge Phys.Rev.B65，165401(2002)中找到。因此，我们现在对如何计算双探针系统的电子密度，包括将外部电压施加在系统上的状况加以描述。

计算双探针系统中的有效单电子势能函数

在前面的部分中，我们示出了如何利用格林函数方法从哈密顿函数中计算电子密度。为了完成自洽循环，我们需要从电子密度中计算哈密顿函数，这意味着计算有效单电子势能函数V^eff[n]。在DFT中，有效单电子势能函数由方程(5)给出，对于双探针系统，我们需要针对相互作用区和电极区分别求解泊松方程，即，方程(8)。电极的哈特里势能函数利用与用于周期性系统相同的方法获得，在这种情况下，重复结构是用于定义方程(32)中的H_L1L1的电极晶胞(electrode cell)和用于定义H_R1R1的右电极的相应晶胞。这些电极哈特里势能函数现在提供相互作用区的哈特里势能函数的边界条件。但是，电极是体积系统(bulk system)，这意味着我们可以将任何常数加入它们的哈特里势能函数中，仍然获得有效解。为了固定这个任意常数，我们将每个电极哈特里势能函数与电极的化学势相联系，并利用方程(15)与左右化学势相联系。因此，我们现在固定了电极中的哈特里势，它们限定了中心区中的泊松方程沿着z方向的边界条件。在x和y方向，我们使用周期性边界条件。借助于这些边界条件，可以通过在Taylor Phys.Rev.B63，245407(2001)中所述的多网格方法获得相互作用区的哈特里势能函数。

从格林函数中获得的电子输运系数和电流

在完成了自洽循环之后，我们可以计算系统的传递特性。通过触点的非线性电流I像如下那样获得：

方程(41)

其中， $G_0=2\frac{e^2}{h}$ 是传导量子。这个方程的理论基础描述在H.Haug的《半导体的量子输运动力学和光学》(H.Haug，Quantum kinetics in transport and opticsof semiconductors，Springer-Verlag，Berlin，1996)中。

双探针系统的自洽算法

图5示出了双探针计算由于如方程(15)所述的左右电极之间的施加电压，通过纳米级器件从左电极流到右电极的电流所需的步骤。首先，我们通过指定离子位置和像施加电压和温度那样的外部参数定义系统(202)。接着，我们利用筛选近似将系统几何结构分解成相互作用区和电极区(204)。在电极区中，电子密度和有效单电子势能函数应该接近它们的体积值(bulk value)。通常，在金属表面的第三层中的原子附近就是这种情况，因此，使金属表面的前两层包括在相互作用区内就足够了。我们利用图1中的流程图计算分离电极区的自洽有效单电子势能函数(206)。我们利用方程(4)、(12)和(33)-(36)从自洽有效单电子势能函数中构建电极格林函数，并利用方程(29)和(30)构建电极自能(208)。现在，将这些初始计算值用作双探针计算的输入。因此，我们已经计算出电极区的自洽密度，和只需计算相互作用区的自洽密度。从首先猜测相互作用区的电子密度开始(210)，我们进行与图1中的流程图相似的自洽循环。首先，我们利用方程(5)-(8)计算相互作用区的有效单电子势能函数(212)。我们可以利用方程(4)和(12)从有效单电子势能函数中获取哈密顿函数，并通过方程(31)获取格林函数(214)。我们可以利用方程(23)、(38)、(39)和(40)从格林函数中计算电子密度，从而结束自洽循环(218)。如果新电子密度不同于(在指定数值精度内)用于构建有效单电子势能函数的电子密度(220)，我们就根据前面计算的密度作出新的改进猜测。在最简单的形式中，以混合参数β线性混合两个密度获取新猜测(216)。如果输入和输出密度相等，我们就获得了电子密度的自洽值，从而也获得了有效单电子势能函数、哈密顿函数和格林函数(222)。我们可以利用方程41从这个格林函数中计算电流。在计算出电流之后，结束该算法(226)。

该进程已经在TranSIESTA和McDCAL软件中得到实现。并且，对这些软件和实现细节的描述可以在Brandbyge Phys.Rev.B65，165401(2002)和Taylor Phys.Rev.B63，245407(2001)中找到。为了获得纳米级器件的伏安特性I-U曲线，我们需要对每个电压U进行自洽计算。这例示在图6的流程图中。输入系统几何条件和电压区间U₀-U₁和步长ΔU(302)。将开始电压设置成U₀(304)。在图5的流程图2中的步骤之后，在电压U上进行有效单电子势能函数的自洽计算，并将自洽势能函数用于计算电流(306)。按步长升高电压(308)，如果新电压在指定电压区间内，那么进行新自洽计算(310)，否则，停止(312)。

例子：与金表面耦合的DTP的I-U特性的计算

现在，我们将给出利用TranSIESTA软件计算例示在图2中的几何结构的I-U特性的结果。该计算遵循图6的流程图3，图9中的点示出了计算结果。在Stokbro的《计算材料科学》(Stokbro，Computational Materials Science 27，151(2003))中获得了相似的I-U特性。

在图7中，我们示出了由施加电压引起的自洽有效单电子势能函数的变化。有效单电子势能函数的值是沿着从左电极开始，经过DTB分子的两个硫原子的中心，终止在右电极上的线段示出的。在右电极中，由于施加的电压，有效单电子势能函数下移了(shift down)。主要特征是，在电极区内，有效单电子势能函数是平坦的，而在分子区内，发生主电压降。

图7a中的曲线都具有相似的形状。在图7b中，我们用施加电压重新调整了曲线，可以观察到，重新调整的有效单电子势能函数几乎相同。这种观察形成本发明的基础，因为它表明，有效单电子势能函数的自洽变化与施加的电压存在简单的变化关系。

利用二电压点的线性内插

在该算法的一种形式中，在零电压U₀上和针对小有限电压U_Δ计算有效单电子势能函数。现在将这些数据用于外推一般电压。一般电压U的有效单电子势能函数通过简单线性外推获得：

$V_{\mathrm{int}}^{\mathrm{eff}}\left[U\right]:=V_{\mathrm{SCF}}^{\mathrm{eff}}\left[U_0\right]+\frac{U-U_0}{U_{\mathrm{\Δ}}-U_0}\left(V_{\mathrm{SCF}}^{\mathrm{eff}}\right[U_{\mathrm{\Δ}}]-V_{\mathrm{SCF}}^{\mathrm{eff}}[U_0\left]\right).$ 方程(42)

哈密顿函数通过下式与有效单电子势能函数联系：

方程(43)

这意味着相同的换算关系应用于哈密顿函数。因此，一般电压上的哈密顿函数可以通过下式近似：

${\hat{H}}_{\mathrm{int}}\left[U\right]:={\hat{H}}_{\mathrm{SCF}}\left[U_0\right]+\frac{U-U_0}{U_{\mathrm{\Δ}}-U_0}\left({\hat{H}}_{\mathrm{SCF}}\right[U_{\mathrm{\Δ}}]-{\hat{H}}_{\mathrm{SCF}}[U_0\left]\right),$ 方程(44)

其中， 和

是U₀和U_Δ上的自洽哈密顿函数。

在大多数电子结构方法中，哈密顿函数用基组{φ_i}展开，并通过如下矩阵表示：

${\stackrel{\&OverBar;}{H}}_{\mathrm{ij}}=<{\mathrm{\&phi;}}_i\left|\hat{H}\right|{\mathrm{\&phi;}}_j>.$ 方程(45)

在这种情况下，将线性内插公式应用于哈密顿矩阵元。

我们可以从哈密顿函数中计算系统的所有特性，包括由施加电压引起的电流。该电流通过首先利用方程(31)计算格林函数，然后利用方程(41)从格林函数中计算电流获得。我们可以将方程(44)、(31)和(41)组合在一起，并将其写成采用

U，和在电压U上返回I的映射M。

我们将该映射写成：

$I\left(U\right):M(U,{\stackrel{\&OverBar;}{H}}_{\mathrm{SCF}}[U_0],{\stackrel{\&OverBar;}{H}}_{\mathrm{SCF}}[U_{\mathrm{\&Delta;}}\left]\right),$ 方程(46)

利用内插公式计算I-U特性总结在图8的流程图4中。输入系统几何条件和电压区间U₁-U₂、步长ΔU以及我们将计算用于内插的自洽哈密顿函数的电压U₀、U_Δ(402)。利用图5的流程图2计算自洽有效单电子势能函数和电压U₀的哈密顿函数(404)。对电压U_Δ进行自洽计算(406)。利用图9的流程图5计算电压区间U₁-U₂的I-U曲线，以便利用方程(46)，以及U₀和U_Δ上的自洽结果获取电流的近似值(408)。然后停止(410)。

I-U曲线的计算遵循图9的流程图5。输入电压区间U₁-U₂、步长ΔU和两个电压U₀和U_Δ的自洽哈密顿函数(502)。将开始电压设置成U：＝U₁(504)。利用方程(46)以及U₀和U_Δ上的自洽结果获取U上的电流的近似值(506)。按步长(step size)升高电压(508)，如果新电压在指定电压区间内，那么继续计算I-U曲线(510)，否则，停止(512)。

用于计算的典型参数是选择U₀＝0伏特和U_Δ＝0.4伏特。由于非平衡密度，即，方程(40)的计算涉及到点数与电压大小成正比的积分，自洽计算在高压上需要更多的计算，因此，选择相对低的电压值是最有效的。

I-U曲线中电压范围的典型值是U₁＝-2.0伏特和U₂＝2.0伏特。对于纳米级小器件，在较高电压上，电场强度非常高，而且由于器件的电击穿，这样的电压在实验上难以测量。

在图10中，我们将利用方程(46)中的公式计算电流的结果与完全自洽解相比较。用“1阶”表示的线段示出了利用方程(46)获得的结果，而用“SCF”表示的线段示出了利用自洽计算获得的结果。可以看出，尽管只有V＝0.0伏特和V＝0.4伏特处的计算用于计算，但对于V＜2.0伏特，利用方程(46)获得的结果与完全自洽计算极其一致。

计算I-U特性的自适应网格方法

在前面的部分中，我们使用了利用两个电压U₀和U_Δ上的自洽哈密顿函数，将哈密顿函数外推到一般电压的二点内插公式。现在，我们将提出改进这种方案的系统方法。该方法基于在所选电压点上进行附加自洽计算，并利用这些电压点上的自洽哈密顿函数形成改进内插公式。借助于这种方法，生成收敛于自洽计算I-U特性的一系列I-U曲线。

我们的目的是计算区间[U₁，U₂]内的I-U特性。图11中的流程图6示出了计算所涉及的步骤。最初步骤类似于图8的流程图4；但是，在这种新算法中，我们将通过进行附加自洽计算改善近似，其中，新电压点可以通过如图12和13的流程图7和8所示的算法选择。输入系统几何条件和电压区间U₁-U₂、步长ΔU和内插电压U₀、U_Δ(602)。利用图5的流程图2计算自洽有效单电子势能函数和电压U₀的哈密顿函数(604)。对电压U_Δ进行自洽计算(606)。利用图13的流程图8计算电压区间U₁-U₀的I-U曲线，以便利用方程(46)以及U₀和U_Δ上的自洽结果获取电流的近似值(608)。利用图12的流程图7计算电压区间U₀-U₂的I-U曲线，以便利用方程(46)以及U₀和U_Δ上的自洽结果获取电流的近似值(610)。然后停止(612)。

图12和图13的流程图7和8示出了细分区间的算法。细分区间，直到内插电流和自洽计算电流在我们用δ表示的指定精度内一致。除了流程图7假设自洽哈密顿函数在我们需要I-U曲线的电压区间的最低电压U_A上已知，而流程图8假设自洽哈密顿函数在电压区间的最高电压U_B上已知之外，流程图7和8是相似的。对于流程图7，到递归步骤的输入是电压区间U_A-U_B和端点U_A和任意电压点U_C处的自洽哈密顿函数(702)。接着，我们在电压区间的最高电压U_B处进行自洽计算(704)。我们从内插公式，即，方程(46)和自洽哈密顿函数，即，方程(31)和(41)中计算电流(706)。如果内插电流与自洽电流相差超过δ(708)，我们将区间进一步细分成区间{U_A，U_M}和{U_M，U_B}，其中，U_M：＝(U_A+U_B)/2(714)。对于{U_A，U_M}，递归地调用算法(710)。对于{U_M，U_B}，我们知道最后电压点上的哈密顿函数，而不是知道第一电压点上的哈密顿函数，所以我们使用如图8所示的轻微修正算法(718)。该过程一直持续到新网格点的自洽计算电流在规定精度δ内与内插值一致。当获得规定的精度时，我们可以安全地使用方程(46)来计算子区间{U_A，U_B}的I-U特性(710)。然后，停止递归算法(712)。

图13的流程图8中的算法是图12的流程图7中的算法的轻微修正，唯一不同之处在于，取代U_A，在U_B上计算输入的自洽哈密顿函数。这里，我们仅仅提到流程图8与流程图7相比的不同之处。取代H_SCF[U_A]，输入H_SCF[U_B](802)。取代U_B，在U_A上进行自洽计算(804)。在U_A上计算电流(806)，比较在U_A上计算的电流(808)。该算法的其余部分与流程图7的算法相似。

我们注意到，一般说来，这个过程将导致网格点不均匀地分布在电压窗上。网格点主要集中在线性内插公式不能很好地描述自洽势能函数的变化的区域中。因此，该算法导致网格点自适应地形成。

利用高阶近似

对于在前面部分中所述的方法，通过进行附加自洽计算系统性地改善近似解。当进行不止两次的自洽计算时，可以使用高阶内插公式。例如，对于有效势V_int ^eff[U]，可以组合U₀、U₁和U₂上的自洽计算，以获得二阶外推公式：

$V_{\mathrm{int}}^{\mathrm{eff}}\left[U\right]:=V_{\mathrm{SCF}}^{\mathrm{eff}}\left[U_0\right]+\left({U-U}_0\right)b+{(U-U_0)}^{2}c,$ 方程(46b)

$c=\left(V_{\mathrm{SCF}}^{\mathrm{eff}}\right[U_1]-\frac{U_1-U_0}{U_2-U_0}{V}_{\mathrm{SCF}}^{\mathrm{eff}}[U_2\left]\right)/(U_2{U}_2-U_1{U}_1),$ 方程(46c)

$b=V_{\mathrm{SCF}}^{\mathrm{eff}}\left[U_1\right]/(U_1-U_0)-c(U_1-U_0).$ 方程(46d)

类似的二阶外推公式可以用于哈密顿函数：

H[U]：＝H[U₀]+(U-U₀)b+(U-U₀)²c，                    方程(46e)

$c=\left(H\right[U_1]-\frac{U_1-U_0}{U_2-U_0}H[U_2\left]\right)/(U_2{U}_2-U_1{U}_1),$ 方程(46f)

b＝H[U₁]/(U₁-U₀)-c(U₁-U₀)。                         方程(46g)

图10中用“2阶”表示的线段示出了利用从0.0伏特、0.4伏特和1.0伏特处的自洽计算获得的二阶外推公式所得的结果。

上面的描述可以容易地推广到，对于n个偏压，使用(n-1)阶外推公式。

推广到多探针系统

可以将算法推广到多探针系统，即，存在不止两个探针的系统。让我们假设系统还包括一个附加电极，那么我们可以通过电极之间的施加电压U^L3将这个电极的化学势μ₃与左电极的化学势相联系：

μ_L-μ₃＝eU^L3。                                      方程(47)

现在，我们可以将方程(44)推广成变量为U^L3和U^LR的二维内插公式，

其中，后者是左电极与右电极之间的电压差。选择 $U_0^{L3}=U_0^{\mathrm{LR}}=U_0=0$ 是方便的，因为，我们以后可以将相同的自洽哈密顿函数用于内插公式中的U₀值。在这种情况下，

${\hat{H}}_{\mathrm{int}}[U^{L3},U^{\mathrm{LR}}]:{=\hat{H}}_{\mathrm{SCF}}\left[U_0\right]+\frac{U^{L3}{-U}_0}{U_{\mathrm{\&Delta;}}^{L3}-U_0}\left({\hat{H}}_{\mathrm{SCF}}\right[U_{\mathrm{\&Delta;}}^{L3}]-{\hat{H}}_{\mathrm{SCF}}[U_0\left]\right),$

$+\frac{U^{\mathrm{LR}}-U_0}{U_{\mathrm{\&Delta;}}^{\mathrm{LR}}-U_0}\left({\hat{H}}_{\mathrm{SCF}}\right[U_{\mathrm{\&Delta;}}^{\mathrm{LR}}]-{\hat{H}}_{\mathrm{SCF}}[U_0\left]\right)$ 方程(48)

其中，U_Δ ^L3和U_Δ ^LR分别是左电极-电极3电压和左电极-右电极电压的小电压增量。针对 $U^{L3}=U_{\mathrm{\&Delta;}}^{L3}$ 和U^LR＝0计算自洽哈密顿函数

并针对 $U^{\mathrm{LR}}=U_{\mathrm{\&Delta;}}^{\mathrm{LR}}$ 和U^L3＝0计算

推广到使用电子或离子温度

到目前为止，我们暗中假设电子温度是零，因为所有积分都被写成积分边界固定在化学势上。为了包括有限电子温度，我们必须改变方程(18)、(24)、(37)、(39)、(40)、(41)中的积分，以使得

$\stackrel{\mathrm{\&mu;}}{\&Integral;}\&RightArrow;\stackrel{\&infin;}{\&Integral;}f[(\mathrm{\&epsiv;}-\mathrm{\&mu;})/\mathrm{kT}],$ 方程(49)

其中，T是温度，k是玻尔兹曼(Boltzmann)常数，f是费米(Fermi)函数

$f\left[x\right]=\frac{1}{e^x+1}.$ 方程(50)

通过在费米函数中对左右电极使用不同的T值，我们可以容易地将这种情况推广到对左右电极使用不同的电子温度。

本领域的普通技术人员应该认识到，本发明不局限于本文具体所示和所述的内容，可以对本发明的优选实施例作出许多修正和改变，而不偏离本发明的精神和范围。

Claims (56)

1.一种根据一个或多个外部参数的值的变化来利用外推分析表达近似自洽解或自洽解的变化的方法，所述自洽解用在至少含有两个探针或电极的系统的模型中，该模型基于包含有效单电子势能函数和/或有效单电子哈密顿函数的自洽确定的电子结构计算，该方法包含：

利用自洽循环计算为第一外部参数的第一值确定对所选函数的第一自洽解；

利用自洽循环计算为第一所选外部参数的第二值确定对所选函数的第二自洽解，所述第一所选外部参数的第二值不同于第一所选外部参数的第一值；和

根据至少所确定的第一和第二自洽解以及第一所选外部参数的第一和第二值，利用外推表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化。

2.根据权利要求1所述的方法，其中，利用线性外推表达近似自洽解或自洽解的变化。

3.根据权利要求1所述的方法，其中，

利用自洽循环计算为第一所选外部参数的第三值确定对所选函数的第三自洽解，所述第一所选外部参数的第三值不同于第一所选外部参数的第一和第二值，并且其中

根据至少所确定的第一、第二和第三自洽解以及第一所选外部参数的第一、第二和第三值，利用外推表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化。

4.根据权利要求3所述的方法，其中，利用二阶外推表达近似自洽解或自洽解的变化。

5.根据权利要求1-4的任何一项所述的方法，其中，被模拟的系统是纳米级器件或包含纳米级器件的系统。

6.根据权利要求1-5的任何一项所述的方法，其中，系统的模拟包含将一个或多个外部参数作为输入提供给所述探针或电极。

7.根据权利要求1-6的任何一项所述的方法，其中，该系统是双探针系统，而外部参数是跨越所述两个探针或电极的偏压U，所述双探针系统被模拟成具有通过相互作用区相互耦合的两个实质上半无限的探针或电极。

8.根据权利要求1-6的任何一项所述的方法，其中，该系统是带有三个探针或电极的三探针系统，而外部参数是第一所选参数和类型与第一所选参数相同的第二所选参数。

9.根据权利要求8所述的方法，其中，该系统是带有三个探针或电极的三探针系统，而外部参数是跨越所述电极中的第一和第二电极的第一偏压U1以及跨越所述电极中的第三和第一电极的第二偏压U2，所述三探针系统被模拟成具有通过相互作用区相互耦合的三个实质上半无限的探针或电极。

10.根据权利要求8或9所述的方法，所述方法进一步包含：

利用自洽循环计算为第二所选外部参数的第一值确定对所选函数的第四自洽解；

利用自洽循环计算为第二所选外部参数的第二值确定对所选函数的第五自洽解，所述第二所选外部参数的第二值不同于第二所选外部参数的第一值；和

其中，根据至少所确定的第一和第二自洽解，以及第一所选外部参数的第一和第二值，并进一步根据至少所确定的第四和第五自洽解，以及第二所选外部参数的第一和第二值，利用外推为第一所选外部参数的所选值和第二所选外部参数的所选值表达所选函数的近似自洽解或自洽解的变化的所述表达。

11.根据权利要求10所述的方法，其中，利用线性外推表达近似自洽解或自洽解的变化。

12.根据权利要求3和10所述的方法，其中，

利用自洽循环计算为第二所选外部参数的第三值确定对所选函数的第六自洽解，所述第二所选外部参数的第三值不同于第二所选外部参数的第一和第二值；并且其中

根据至少所确定的第一、第二和第三自洽解，以及第一所选外部参数的第一、第二和第三值，并进一步根据至少所确定的第四、第五和第六自洽解，以及第二所选外部参数的第一、第二和第三值，利用外推为第一所选外部参数的所选值和第二所选外部参数的所选值表达所选函数的近似自洽解或自洽解的变化的所述表达。

13.根据权利要求12所述的方法，其中，利用二阶外推表达近似自洽解或自洽解的变化。

14.根据权利要求10-13的任何一项所述的方法，其中，第二所选外部参数的第一值等于第一所选外部参数的第一值。

15.根据权利要求1-14的任何一项所述的方法，其中，从由如下表示的函数中选择所选函数：有效单电子势能函数、有效单电子哈密顿函数和电子密度。

16.根据权利要求15所述的方法，其中，所选函数是有效单电子势能函数或有效单电子哈密顿函数，而自洽循环计算基于密度函数理论DFT或哈特里-福克理论HF。

17.根据权利要求1-16的任何一项所述的方法，其中，自洽循环计算基于包括如下步骤的循环计算：

a)为系统模型的所选区选择电子密度的值；

b)为所选电子密度和外部参数的所选值确定有效单电子势能函数；

c)计算与确定的有效单电子势能函数相对应的电子密度值；

d)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

e)将有效单电子势能函数的解定义成有效单电子势能函数的自洽解；而如果不相等，那么

f)选择新的电子密度值，并重复步骤b)-f)，直到电子密度的所选值和计算值在所述给定数值精度内相等。

18.根据权利要求17所述的方法，其中，针对系统的探针或电极区确定有效单电子势能函数的自洽解。

19.根据权利要求16-18的任何一项所述的方法，其中，所选函数是系统相互作用区的有效单电子哈密顿函数，而系统相互作用区的有效单电子哈密顿函数的第二自洽解的确定包含在第一所选外部参数的给定值上计算相互作用区的有效单电子势能函数的相应自洽解的步骤。

20.根据权利要求18和19所述的方法，其中，根据有效单电子势能函数的相应所确定的自洽解，为每个探针或电极区构建或确定格林函数。

21.根据权利要求19或20所述的方法，其中，有效单电子哈密顿函数的第二自洽解的确定基于包括如下步骤的循环计算：

aa)为系统的相互用作区选择电子密度的值；

bb)为所选外部参数的给定值确定所选电子密度的有效单电子势能函数；

cc)根据在步骤bb)中确定的有效单电子势能函数确定相互作用区的有效单电子哈密顿函数的解；

dd)根据在步骤cc)中确定的有效单电子哈密顿函数的解确定相互作用区的格林函数的解；

ee)计算与相互作用区的确定格林函数相对应的电子密度值；

ff)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

gg)将有效单电子哈密顿函数的解定义成有效单电子哈密顿函数的自洽解；而如果不相等，那么

hh)选择新的电子密度值，并重复步骤bb)-hh)，直到电子密度的所选值和计算值在所述给定数值精度内相等。

22.根据权利要求1-21的任何一项所述的方法，其中，所选函数是用哈密顿矩阵表示的有效单电子哈密顿函数，所述矩阵的每个元素是具有利用相应外推表达式表达的近似自洽解或自洽解的变化的函数。

23.根据权利要求7-22的任何一项所述的方法，其中，所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压，其中，分别针对外部偏压的所选第一和第二值，为有效单电子哈密顿函数确定第一和第二自洽解，从而，当外部偏压发生变化时，获取外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述方法进一步包含：

利用表达有效单电子哈密顿函数的近似自洽解或自洽解的变化的所获得的外推表达式，针对施加偏压的许多不同值确定系统的两个探针之间的电流。

24.根据权利要求23所述的方法，其中，针对外部偏压的给定范围和针对外部偏压的给定电压步长确定电流。

25.根据权利要求24所述的方法，其中，利用如下循环确定电流：

aaa)在外部偏压的给定范围内确定最低电压的电流；

bbb)按给定电压步长升高偏压；

ccc)确定新升偏压的电流；和

ddd)重复步骤bbb)和ccc)，直到新升偏压大于偏压的给定范围的最高电压。

26.根据权利要求7-22的任何一项所述的方法，其中，所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压，所述方法包含：

将外部偏压的所确定的电压范围划分成至少第一和第二电压范围；

为第一和第二电压范围确定与所述电压范围的最大和最小值相对应的有效单电子哈密顿函数的最大和最小自洽解；

当外部偏压发生变化时，获取第一外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述第一外推表达式基于第一电压范围的确定最大和最小自洽解以及第一电压范围的最大和最小电压值；

当外部偏压发生变化时，获取第二外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述第二外推表达式基于第二电压范围的确定最大和最小自洽解以及第二电压范围的最大和最小电压值；

利用获得的第一外推表达式，在由第一电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流；和

利用获得的第二外推表达式，在由第二电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。

27.根据权利要求26所述的方法，其中，将确定的电压范围划分成至少三个电压范围，所述方法进一步包含：

为第三电压范围确定与第三电压范围的最大和最小值相对应的有效单电子哈密顿函数的最大和最小自洽解；

当外部偏压发生变化时，获取第三外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述第三外推表达式基于第三电压范围的确定最大和最小自洽解以及第三电压范围的最大和最小电压值；和

利用获得的第三线性外推表达式，在由第三电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。

28.根据权利要求7-22的任何一项所述的方法，其中，所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压，其中，分别针对外部偏压的所选第一和第二值，为有效单电子哈密顿函数确定第一和第二自洽解，所述第二值高于偏压的所选第一值，从而，当外部偏压发生变化时，获取第一外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述方法进一步包含：

aaaa)为外部偏压选择具有最大值和最小值的电压范围，以便在所述范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流；

bbbb)利用自洽循环计算，针对外部偏压的所选最大值确定有效单电子哈密顿函数的最大自洽解；

cccc)根据相应所确定的最大自洽解，针对偏压的最大值确定系统的两个探针之间的电流；

dddd)根据获得的第一外推表达式，针对偏压的所选最大值确定系统的两个探针之间的电流；

eeee)比较在步骤cccc)和dddd)中确定的电流，如果它们在给定数值精度内相等，那么

ffff)当外部偏压发生变化时，将外推表达式用作有效单电子哈密顿函数的近似自洽解，在由所选第一电压值和最大电压值给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。

29.根据权利要求23所述的方法，其中，获取最大外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述最大外推表达式基于所确定的第一和最大自洽解以及第一偏压和该偏压的最大值，并且其中当在步骤ffff)中确定电流时，使用所述最大外推表达式。

30.根据权利要求28和29所述的方法，其中，当在步骤eeee)中在步骤cccc)和dddd)中确定的电流值在给定数值精度内不相等时，

gggg)在第一值和前最大值之间选择外部偏压的新最大值；和

hhhh)重复步骤bbbb)到hhhh)，直到在步骤cccc)和dddd)中确定的电流值在所述给定数值精度内相等。

31.根据权利要求28-30的任何一项所述的方法，所述方法进一步包含：

iiii)利用自洽循环计算，针对外部偏压的所选最小值计算有效单电子哈密顿函数的最小自洽解；

jjjj)根据相应所确定的最小自洽解，针对偏压的最小值确定系统的两个探针之间的电流；

kkkk)根据获得的第一外推表达式，针对偏压的所选最小值确定系统的两个探针之间的电流；

llll)比较在步骤jjjj)和kkkk)中确定的电流，如果它们在给定数值精度内相等，那么

mmmm)当外部偏压发生变化时，将外推表达式用作有效单电子哈密顿函数的近似自洽解，在由所选第一电压值和最小电压值给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流。

32.根据权利要求31所述的方法，其中，获取最小外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述最小外推表达式基于所确定的第一和最小自洽解以及第一偏压和该偏压的最小值，并且其中当在步骤mmmm)中确定电流时，使用最小外推表达式。

33.根据权利要求31和32所述的方法，其中，当在步骤llll)中在步骤jjjj)和kkkk)中确定的电流值在给定数值精度内不相等时，

nnnn)在第一值和前最小值之间选择外部偏压的新最小值；和

oooo)重复步骤iiii)到oooo)，直到在步骤jjjj)和kkkk)中确定的电流值在所述给定数值精度内相等。

34.一种根据一个或多个外部参数的值的变化来利用外推分析表达近似自洽解或自洽解的变化的计算机系统，所述自洽解用在至少含有两个探针或电极的纳米级系统的模型中，该模型基于包含有效单电子势能函数和/或有效单电子哈密顿函数的自洽确定的电子结构计算，所述计算机系统包含：

利用自洽循环计算为第一外部参数的第一值确定对所选函数的第一自洽解的装置；

利用自洽循环计算为第一所选外部参数的第二值确定对所选函数的第二自洽解的装置，所述第一所选外部参数的第二值不同于第一所选外部参数的第一值；和

根据至少所确定的第一和第二自洽解以及第一所选外部参数的第一和第二值，利用外推表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化的装置。

35.根据权利要求34所述的计算机系统，其中，表达近似自洽解或自洽解的变化的装置适用于利用线性外推表达这样的解。

36.根据权利要求34所述的计算机系统，所述计算机系统进一步包含：

利用自洽循环计算为第一所选外部参数的第三值确定对所选函数的第三自洽解的装置，所述第一所选外部参数的第三值不同于第一所选外部参数的第一和第二值，并且其中

表达第一所选外部参数的至少一个所选值的所选函数的近似自洽解或自洽解的变化的装置适用于根据至少所确定的第一、第二和第三自洽解以及第一所选外部参数的第一、第二和第三值，利用外推表达这样的解。

37.根据权利要求36所述的计算机系统，其中，表达近似自洽解或自洽解的变化的装置适用于利用二阶外推表达这样的解。

38.根据权利要求34-37的任何一项所述的计算机系统，其中，纳米级系统是双探针系统和外部参数是跨越所述两个探针或电极的偏压U，所述双探针系统被模拟成具有通过相互作用区相互耦合的两个实质上半无限的探针或电极。

39.根据权利要求34-37的任何一项所述的计算机系统，其中，纳米级系统是带有三个探针或电极的三探针系统和外部参数是第一所选参数和类型与第一所选参数相同的第二所选参数。

40.根据权利要求39所述的计算机系统，其中，纳米级系统是带有三个探针或电极的三探针系统，而外部参数是跨越所述电极中的第一和第二电极的第一偏压U1以及跨越所述电极中的第三和第一电极的第二偏压U2，所述三探针系统被模拟成具有通过相互作用区相互耦合的三个实质上半无限的电极。

41.根据权利要求39或40所述的计算机系统，其中，所述计算机系统进一步包含：

利用自洽循环计算为第二所选外部参数的第一值确定对所选函数的第四自洽解的装置；

利用自洽循环计算为第二所选外部参数的第二值确定对所选函数的第五自洽解的装置，所述第二所选外部参数的第二值不同于第二所选外部参数的第一值；和

其中，表达所选函数的近似自洽解或自洽解的变化的所述装置适用于根据所确定的第一和第二自洽解，以及第一所选外部参数的第一和第二值，并进一步根据所确定的第四和第五自洽解，以及第二所选外部参数的第一和第二值，利用外推表达第一所选外部参数的所选值和第二所选外部参数的所选值的近似自洽解。

42.根据权利要求41所述的计算机系统，其中，表达近似自洽解或自洽解的变化的装置适用于利用线性外推表达这样的解。

43.根据权利要求36和41所述的计算机系统，所述计算机系统进一步包含：

利用自洽循环计算为第二所选外部参数的第三值确定对所选函数的第六自洽解的装置，所述第二所选外部参数的第三值不同于第二所选外部参数的第一和第二值，并且其中

表达所选函数的近似自洽解或自洽解的变化的装置适用于根据至少所确定的第一、第二和第三自洽解，以及第一所选外部参数的第一、第二和第三值，并进一步根据至少所确定的第四、第五和第六自洽解，以及第二所选外部参数的第一、第二和第三值，利用外推表达第一所选外部参数的所选值和第二所选外部参数的所选值的近似自洽解。

44.根据权利要求43所述的计算机系统，其中，表达近似自洽解或自洽解的变化的装置适用于利用二阶外推表达这样的解。

45.根据权利要求41-44的任何一项所述的计算机系统，其中，第二所选外部参数的第一值等于第一所选外部参数的第一值。

46.根据权利要求34-45的任何一项所述的计算机系统，其中，从由如下表示的函数中选择所选函数：有效单电子势能函数、有效单电子哈密顿函数和电子密度。

47.根据权利要求46所述的计算机系统，其中，所选函数是有效单电子势能函数或有效单电子哈密顿函数，而自洽循环计算基于密度函数理论DFT或哈特里一福克理论HF。

48.根据权利要求34-47的任何一项所述的计算机系统，进一步包含根据包括如下步骤的循环计算进行自洽循环计算的装置：

a)为纳米级系统模型的所选区选择电子密度的值；

b)为所选电子密度和外部参数的所选值确定有效单电子势能函数；

c)计算与确定的有效单电子势能函数相对应的电子密度；

d)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

e)将有效单电子势能函数的解定义成有效单电子势能函数的自洽解；而如果不相等，那么

f)选择新的电子密度值，并重复步骤b)-f)，直到电子密度的所选值和计算值在所述给定数值精度内相等。

49.根据权利要求48所述的计算机系统，其中，进行自洽循环计算的装置适用于针对系统的探针或电极区确定有效单电子势能函数的自洽解。

50.根据权利要求47-49的任何一项所述的计算机系统，其中，所选函数是系统相互作用区的有效单电子哈密顿函数，和确定系统相互作用区的有效单电子哈密顿函数的第二自洽解的装置适用于通过包括在第一所选外部参数的给定值上计算相互作用区的有效单电子势能函数的相应自洽解的步骤进行所述确定。

51.根据权利要求49和50所述的计算机系统，进一步包含根据有效单电子势能函数的相应所确定的自洽解，为每个探针或电极区确定格林函数的装置。

52.根据权利要求50或51所述的计算机系统，其中，确定有效单电子哈密顿函数的第二自洽解的装置适用于根据包括如下步骤的循环计算进行所述确定：

aa)为系统的相互用作区选择电子密度的值；

bb)为所选外部参数的给定值确定所选电子密度的有效单电子势能函数；

cc)根据在步骤bb)中确定的有效单电子势能函数确定相互作用区的有效单电子哈密顿函数的解；

dd)根据在步骤cc)中确定的有效单电子哈密顿函数的解确定相互作用区的格林函数的解；

ee)计算与相互作用区的确定格林函数相对应的电子密度值；

ff)将电子密度的所选值与电子密度的计算值相比较，如果电子密度的所选值和计算值在给定数值精度内相等，那么

gg)将有效单电子哈密顿函数的解定义成有效单电子哈密顿函数的自洽解；而如果不相等，那么

hh)选择新的电子密度值，并重复步骤bb)-hh)，直到电子密度的所选值和计算值在所述给定数值精度内相等。

53.根据权利要求38-52的任何一项所述的计算机系统，其中，所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压，其中，确定第一和第二自洽解的装置适用于分别针对外部偏压的所选第一和第二值，为有效单电子哈密顿函数进行所述确定，并且其中利用外推分析表达近似自洽解的装置适用于当外部偏压发生变化时，获取外推表达式，作为有效单电子哈密顿函数的近似自洽解，所述计算机系统进一步包含：

利用表达有效单电子哈密顿函数的近似自洽解或自洽解的变化的所获得的外推表达式，针对施加偏压的许多不同值确定系统的两个探针之间的电流的装置。

54.根据权利要求53所述的计算机系统，其中，确定电流的装置适用于针对外部偏压的给定范围和针对外部偏压的给定电压步长确定电流。

55.根据权利要求54所述的计算机系统，其中，确定电流的装置适用于利用如下循环进行所述确定：

aaa)在外部偏压的给定范围内确定最低电压的电流；

bbb)按给定电压步长升高偏压；

ccc)确定新升偏压的电流；和

ddd)重复步骤bbb)和ccc)，直到新升偏压大于偏压的给定范围的最高电压。

56.根据权利要求38-52的任何一项所述的计算机系统，其中，所选函数是有效单电子哈密顿函数，而外部参数是跨越系统的两个探针的偏压，所述计算机系统进一步包含：

将外部偏压的所确定的电压范围划分成至少第一和第二电压范围的装置；

为第一和第二电压范围确定与所述电压范围的最大和最小值相对应的有效单电子哈密顿函数的最大和最小自洽解的装置；

当外部偏压发生变化时，获取第一外推表达式，作为有效单电子哈密顿函数的近似自洽解的装置，所述第一外推表达式基于第一电压范围的确定最大和最小自洽解以及第一电压范围的最大和最小电压值；

当外部偏压发生变化时，获取第二外推表达式，作为有效单电子哈密顿函数的近似自洽解的装置，所述第二外推表达式基于第二电压范围的确定最大和最小自洽解以及第二电压范围的最大和最小电压值；

利用获得的第一外推表达式，在由第一电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流的装置；和

利用获得的第二外推表达式，在由第二电压范围的最小和最大电压给出的电压范围内，针对施加偏压的许多不同值确定系统的两个探针之间的电流的装置。

Images