FiELD OF THE INVENTION

The present invention relates to methods and systems for using extrapolation analysis or techniques to express an approximate self-consistent solution or a change in a self- consistent solution based on a change in the value of one or more external parameters. The self-consistent solution may be used in a model of a system or nano-scale system having at least two probes or electrodes, and the model may be based on an electronic structure calculation comprising a self-consistent determination of an effective one- electron potential energy function and/or an effective one-electron Hamiltonian.

BACKGROUND OF THE INVENTION

Most common examples of methods within the field of atomic scale modelling, where the modelling is based on electronic structure calculations that require a self-consistent de¬ termination of an effective one-electron potential energy function are Density Functional Theory (DFT) and Hartree-Fock (HF) theory. Many applications of DFT are studies of how a system responds when an external parameter is varied. In such studies, it is necessary to perform a self-consistent calculation for each value of the external parameter, and this can be very time consuming. An important application is the calculation of the current- voltage (I-U) characteristics of a nano-scale device. An example of such a calculation is given in Stokbro, Computational Materials Science 27, 151 (2003), where the I-U charac- teristics of a Di-Thiol-Phenyl (DTP) molecule coupled with gold surfaces is calculated. The system is illustrated in Fig. 2, and the calculation follows the steps outlined in flowcharts 2 and 3 shown in Figs. 5 and 6. The calculation is very computationally demanding, due to the self-consistent loop for each voltage.

It is an objective of the present invention is to provide an efficient and reasonable accu¬ rate method for determining a change in a self-consistent solution caused by a variation of one or more external parameters. SUMMARY OF THE INVENTION

According to the present invention there is provided a method of using extrapolation analysis or technique to express an approximate self-consistent solution or a change in a self-consistent solution based on a change in the value of one or more external parame¬ ters, said self-consistent solution being used in a model of a system having at least two probes or electrodes, which model is based on an electronic structure calculation compris¬ ing a self-consistent determination of an effective one-electron potential energy function and/or an effective one-electron Hamiltonian, the method comprising: determining a first self-consistent solution to a selected function for a first value of a first external parameter by use of self-consistent loop calculation; determining a second self-consistent solution to the selected function for a second value of the first selected external parameter by use of self-consistent loop calculation, said second value of the first selected external parameter being different to the first value of the first selected external parameter; and expressing an approximate self-consistent solution or a change in the self-consistent solu¬ tion for the selected function for at least one selected value of the first selected external parameter by use of extrapolation based on at least the determined first and second self- consistent solutions and the first and second values of the first selected external parame- ter. Here, the approximate self-consistent solution or change in the self-consistent solution may be expressed by use of linear extrapolation.

According to an embodiment of the invention the method may further comprise that a third self-consistent solution to the selected function is determined for a third value of the first selected external parameter by use of self-consistent loop calculation, said third value of the first selected external parameter being different to the first and second values of the first selected external parameter. Here, the approximate self-consistent solution or change in the self-consistent solution for the selected function for at least one selected value of the first selected external parameter may be expressed by use of extrapolation based on at least the determined first, second and third self-consistent solutions and the first, sec¬ ond and third values of the first selected external parameter. Here, it is preferred that the approximate self-consistent solution or change in the self-consistent solution is expressed by use of second order extrapolation. It is preferred that the system being modelled is a nano-scale device or a system compris¬ ing a nano-scale device. It is also preferred that the modelling of the system comprises providing one or more of the external parameters as inputs to said probes or electrodes.

According to an embodiment of the invention the system being modelled is a two-probe system and the external parameter is a voltage bias, U, across said two probes or elec¬ trodes, said two-probe system being modelled as having two substantially semi-infinite probes or electrodes being coupled to each other via an interaction region.

It is also within an embodiment of the invention that the system being modelled is a three- probe system with three probes or electrodes and the external parameters are a first se¬ lected parameter and a second selected parameter being of the same type as the first selected parameter. Here, the system being modelled may be a three-probe system with three probes or electrodes and the external parameters are a first voltage bias, U1 , across a first and a second of said electrodes and a second voltage bias, U2, across a third and the first of said electrodes, said three-probe system being modelled as having three sub¬ stantially semi-infinite electrodes being coupled to each other via an interaction region.

When the system being modelled is a three-probe system, the method of the invention may further comprise: determining a fourth self-consistent solution to the selected function for a first value of the second selected external parameter by use of self-consistent loop calculation, determining a fifth self-consistent solution to the selected function for a second value of the second selected external parameter by use of self-consistent loop calculation, said second value of the second selected external parameter being different to the first value of the second selected external parameter; and wherein said expressing of the approximate self-consistent solution or change in the self- consistent solution for the selected function is expressed for the selected value of the first selected external parameter and a selected value of the second selected external para- meter by use of extrapolation based on at least the determined first and second self- consistent solutions together with the first and second values of the first selected external parameter, and further based on at least the determined fourth and fifth self-consistent solutions together with the first and second values of the second selected external para¬ meter. Here, the approximate self-consistent solution or change in the self-consistent so- lution may be expressed by use of linear extrapolation. The above described method of the invention provided for the three-probe system may further comprise that a sixth self-consistent solution to the selected function is determined for a third value of the second selected external parameter by use of self-consistent loop calculation, said third value of the second selected external parameter being different to the first and second values of the second selected external parameter; and that said expressing of the approximate self-consistent solution or change in the self- consistent solution for the selected function is expressed for the selected value of the first selected external parameter and the selected value of the second selected external pa- rameter by use of extrapolation based on at least the determined first, second and third self-consistent solutions together with the first, second and third values of the first se¬ lected external parameter, and further based on at least the determined fourth, fifth and sixth self-consistent solutions together with the first, second and third values of the sec¬ ond selected external parameter. Here, the approximate self-consistent solution or change in the self-consistent solution may be expressed by use of second order extrapo¬ lation.

For the methods of the invention provided for the three-probe system, the first value of the second selected external parameter may be equal to the first value of the first selected external parameter.

According to the present invention it is preferred that the selected function is selected from the functions represented by: the effective one-electron potential energy function, the effective one-electron Hamiltonian, and the electron density. Here, it is again preferred that the selected function is the effective one-electron potential energy function or the ef¬ fective one-electron Hamiltonian and the self-consistent loop calculation is based on the Density Functional Theory, DFT, or the Hartree-Fock Theory, HF.

According to an embodiment of the invention, the self-consistent loop calculation may be based on a loop calculation including the steps of: a) selecting a value of the electron density for a selected region of the model of the system, b) determining the effective one-electron potential energy function for the selected electron density and for a selected value of the external parameter, c) calculating a value for the electron density corresponding to the determined effec¬ tive one-electron potential energy function, d) comparing the selected value of the electron density with the calculated value of the electron density, and if the selected value and the calculated value of electron density are equal within a given numerical accuracy, then e) defining the solution to the effective one-electron potential energy function as the self-consistent solution to the effective one-electron potential energy function, and if not, then f) selecting a new value of the electron density and repeat steps b)- f) until the se- lected value and the calculated value of electron density are equal within said given nu¬ merical accuracy. Here, the self-consistent solution to the effective one-electron potential energy function may be determined for the probe or electrode regions of the system.

For embodiments where the self-consistent solution to the effective one-electron potential energy function is determined for the probe or electrode regions of the system, it is also preferred that Green's functions are constructed or determined for each of the probe or electrode regions based on the corresponding determined self-consistent solution to the effective one-electron potential energy function.

It is within an embodiment of the method of the invention that the selected function is the effective one-electron Hamiltonian for an interaction region of the system, and the deter¬ mination of a second self-consistent solution to the effective one-electron Hamiltonian of the interaction region of the system comprises the step of calculating a corresponding self-consistent solution to the effective one-electron potential energy function for the inter- action region at a given value of the first selected external parameter. Here, the determi¬ nation of a second self-consistent solution to the effective one-electron Hamiltonian may be based on a loop calculation including the steps of: aa) selecting a value of the electron density for the interaction region of the system, bb) determining the effective one-electron potential energy function for the selected electron density for a given value of the selected external parameter, cc) determining a solution to the effective one-electron Hamiltonian for the interaction region based on the in step bb) determined effective one-electron potential energy func¬ tion, dd) determining a solution to Green's function for the interaction region based on the in step cc) determined solution to the effective one-electron Hamiltonian, ee) calculating a value for the electron density corresponding to the determined Green's function for the interaction region, ff) comparing the selected value of the electron density with the calculated value of the electron density, and if the selected value and the calculated value of electron density are equal within a given numerical accuracy, then gg) defining the solution to the effective one-electron Hamiltonian as the self-consis¬ tent solution to the effective one-electron Hamiltonian, and if not, then hh) selecting a new value of the electron density and repeat steps bb)- hh) until the selected value and the calculated value of electron density are equal within said given numerical accuracy.

According to an embodiment of the invention the selected function may be the effective one-electron Hamiltonian being represented by a Hamiltonian matrix with each element of said matrix being a function having an approximate self-consistent solution or a change in the self-consistent solution being expressed by use of a corresponding extrapolation ex¬ pression,

The method of the present invention also covers an embodiment wherein the selected function is the effective one-electron Hamiltonian and the external parameter is a voltage bias across two probes of the system, an wherein a first and a second self-consistent so¬ lution is determined for the effective one-electron Hamiltonian for selected first and sec¬ ond values, respectively, of the external voltage bias, whereby an extrapolation expres¬ sion is obtained to an approximate self-consistent solution for the effective one-electron Hamiltonian when the external voltage bias is changed, said method further comprising: determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias using the obtained extrapolation expression, which expresses the approximate self-consistent solution or change in the self-consistent solution for the effective one-electron Hamiltonian. Here, the obtained extrapolation ex¬ pression may be a linear expression. The electrical current may be determined for a given range of the external voltage bias and for a given voltage step in the external voltage bias, and the electrical current may be determined using the following loop: aaa) determining the current for the lowest voltage within the given range of the exter¬ nal voltage bias, bbb) increasing the voltage bias by the given voltage step, ccc) determining the current for the new increased voltage bias, ddd) repeating steps bbb) and ccc) until the new increased voltage bias is larger than the highest voltage of the given range of the voltage bias.

It is also within an embodiment of the invention that the system being modelled is a two probe system and that the selected function is the effective one-electron Hamiltonian and the external parameter is a voltage bias across two probes of the system, said method comprising: dividing a determined voltage range for the external voltage bias in at least a first and a second voltage range, determining for the first and second voltage ranges a maximum and a minimum self-consistent solution to the effective one-electron Hamiltonian corresponding to the maximum and minimum values of said voltage ranges, obtaining a first extrapolation expression to the approximate self-consistent solu¬ tion for the effective one-electron Hamiltonian when the external voltage bias is changed, said first extrapolation expression being based on the determined maximum and minimum self-consistent solutions for the first voltage range and the maximum and minimum volt¬ age values of the first voltage range, obtaining a second extrapolation expression to the approximate self-consistent so¬ lution for the effective one-electron Hamiltonian when the external voltage bias is changed, said second extrapolation expression being based on the determined maximum and minimum self-consistent solutions for the second voltage range and the maximum and minimum voltage values of the second voltage range, determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias within the voltage range given by the minimum and maximum voltage of the first voltage range using the obtained first ex¬ trapolation expression, and determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias within the voltage range given by the minimum and maximum voltage of the second voltage range using the obtained second extrapola- tion expression. Here, the obtained first and second extrapolation expressions may be first and second linear expressions, respectively. It is also within an embodiment of the method of the invention that the determined voltage range is divided in at least three volt¬ age ranges, and that the method further comprises: determining for the third voltage range a maximum and a minimum self-consistent solution to the effective one-eiectron Hamiltonian corresponding to the maximum and minimum values of the third voltage range, obtaining a third extrapolation expression to the approximate self-consistent solu- tion for the effective one-electron Hamiltonian when the external voltage bias is changed, said third extrapolation expression being based on the determined maximum and mini¬ mum self-consistent solutions for the third voltage range and the maximum and minimum voltage values of the third voltage range, and determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias within the voltage range given by the minimum and maximum voltage of the third voltage range using the obtained third linear extrapola¬ tion. Also here, the obtained third extrapolation expression may be a third linear extrapola¬ tion expression.

The method of the present invention also covers an embodiment where the system being modelled is a two-probe system and wherein the selected function is the effective one- electron Hamiltonian and the external parameter is a voltage bias across two probes of the system, an wherein a first and a second self-consistent solution is determined for the effective one-electron Hamiltonian for selected first and second values, respectively, of the external voltage bias, with said second value being higher than the selected first value of the voltage bias, whereby a first extrapolation expression is obtained to an approximate self-consistent solution for the effective one-electron Hamiltonian when the external volt¬ age bias is changed, said method further comprising: aaaa) selecting a voltage range having a minimum value and a maximum value for the external voltage bias in order to determine the electrical current between the two probes of the system for a number of different values of the applied voltage bias within said range, bbbb) determining a maximum self-consistent solution to the effective one-electron Ham¬ iltonian for the selected maximum value of the external voltage bias by use of self-consis- tent loop calculation, cccc) determining the electrical current between the two probes of the system for the maximum value of the voltage bias based on the corresponding determined maximum self-consistent solution, dddd) determining the electrical current between the two probes of the system for the selected maximum value of the voltage bias based on the obtained first extrapolation ex¬ pression, eeee) comparing the current values determined in steps cccc) and dddd), and if they are equal within a given numerical accuracy, then ffff) determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias within the voltage range given by the selected first voltage value and the maximum voltage value using an extrapolation expression for an approximate self-consistent solution for the effective one-electron Ham- iltonian when the external voltage bias is changed. Here, the obtained first extrapolation expression may be a first linear extrapolation expression, and linear extrapolation may be used in step ffff) for expressing the approximate self-consistent solution for the effective one-electron Hamiltonian when the external voltage bias is changed. It is within a pre¬ ferred embodiment that a maximum extrapolation expression is obtained to the approxi- mate self-consistent solution for the effective one-electron Hamiltonian, said maximum extrapolation expression being based on the determined first and maximum self-consis¬ tent solutions and the first voltage bias and the maximum value of the voltage bias, and wherein said maximum extrapolation expression is used when determining the current in step ffff). The maximum extrapolation expression may be a maximum linear extrapolation expression. It is also preferred that when in step eeee) the current values determined in steps cccc) and dddd), are not equal within the given numerical accuracy, then the method further comprises: gggg) selecting a new maximum value of the external voltage bias between the first value and the previous maximum value, hhhh) repeating steps bbbb) to hhhh) until the in steps cccc) and dddd) determined cur¬ rent values are equal within said given numerical accuracy. According to an embodiment of the invention, the method may further comprise the steps:

iiii) determining a minimum self-consistent solution to the effective one-electron Hamil- tonian for the selected minimum value of the external voltage bias by use of self-consis¬ tent loop calculation, jjjj) determining the electrical current between the two probes of the system for the minimum value of the voltage bias based on the corresponding determined minimum self- consistent solution, kkkk) determining the electrical current between the two probes of the system for the selected minimum value of the voltage bias based on the obtained first extrapolation ex¬ pression, IHI) comparing the current values determined in steps jjjj) and kkkk), and if they are equal within a given numerical accuracy, then mmmm) determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias within the voltage range given by the selected first voltage value and the minimum voltage value using an extrapolation ex¬ pression for an approximate self-consistent solution for the effective one-electron Hamil- tonian when the external voltage bias is changed. Here, linear extrapolation may be used in step mmmm) for expressing the approximate self-consistent solution for the effective one-electron Hamiltonian when the external voltage bias is changed. Also here, it is within a preferred embodiment that a minimum extrapolation expression is obtained to the ap¬ proximate self-consistent solution for the effective one-electron Hamiltonian, where the minimum extrapolation expression is based on the determined first and minimum self- consistent solutions and the first voltage bias and the minimum value of the voltage bias, and wherein the minimum extrapolation expression is used when determining the current in step mmmm). Here, the minimum extrapolation expression may be a minimum linear extrapolation expression. Also here it is preferred that when in step IHI) the current values determined in steps jjjj) and kkkk), are not equal within the given numerical accuracy, then the method further comprises: nnnn) selecting a new minimum value of the external voltage bias between the first value and the previous minimum value, oooo) repeating steps iiii) to oooo) until the in steps jjjj) and kkkk) determined current values are equal within said given numerical accuracy.

According to the present invention there is also provided a computer system for using extrapolation analysis to express an approximate self-consistent solution or a change in a self-consistent solution based on a change in the value of one or more external parame- ters, said self-consistent solution being used in a model of a nano-scale system having at least two probes or electrodes, which model is based on an electronic structure calcula¬ tion comprising a self-consistent determination of an effective one-electron potential en¬ ergy function and/or an effective one-electron Hamiltonian, said computer system com¬ prising: means for determining a first self-consistent solution to a selected function for a first value of a first external parameter by use of self-consistent loop calculation; means for determining a second self-consistent solution to the selected function for a second value of the first selected external parameter by use of self-consistent loop calcu- lation, said second value of the first selected external parameter being different to the first value of the first selected external parameter; and means for expressing an approximate self-consistent solution or a change in the self- consistent solution for the selected function for at least one selected value of the first se¬ lected external parameter by use of extrapolation based on at least the determined first and second self-consistent solutions and the first and second values of the first selected external parameter. Here, the means for expressing the approximate self-consistent solu¬ tion or change in the self-consistent solution may be adapted for expressing such solution by use of linear extrapolation.

According to an embodiment of the invention the computer system may further comprise: means for determining a third self-consistent solution to the selected function for a third value of the first selected external parameter by use of self-consistent loop calculation, said third value of the first selected external parameter being different to the first and sec¬ ond values of the first selected external parameter. Here, the means for expressing the approximate self-consistent solution or change in the self-consistent solution for the se¬ lected function for at least one selected value of the first selected external parameter may be adapted for expressing such solution by use of extrapolation based on at least the de¬ termined first, second and third self-consistent solutions and the first, second and third values of the first selected external parameter. Here, it is preferred that the means for ex- pressing the approximate self-consistent solution or change in the self-consistent solution is adapted for expressing such solution by use of second order extrapolation.

For the computer system of the invention it is within an embodiment that the nano-scale system is a two-probe system and the external parameter is a voltage bias, U, across said two probes or electrodes, said two-probe system being modelled as having two substan¬ tially semi-infinite probes or electrodes being coupled to each other via an interaction re¬ gion.

The computer system of the invention also covers an embodiment wherein the nano-scale system is a three-probe system with three probes or electrodes and the external parame- ters are a first selected parameter and a second selected parameter being of the same type as the first selected parameter. Here it is preferred that the nano-scale system is a three-probe system with three probes or electrodes and the external parameters are a first voltage bias, U1 , across a first and a second of said electrodes and a second voltage bias, U2, across a third and the first of said electrodes, said three-probe system being modelled as having three substantially semi-infinite electrodes being coupled to each other via an interaction region.

Also here, when the nano-scale system being modelled is a three-probe system, the computer system of the invention may further comprise: means for determining a fourth self-consistent solution to the selected function for a first value of the second selected external parameter by use of self-consistent loop calculation; means for determining a fifth self-consistent solution to the selected function for a second value of the second selected external parameter by use of self-consistent loop calculation, said second value of the second selected external parameter being different to the first value of the second selected external parameter; and wherein said means for expressing of the approximate self-consistent solution or change in the self-consistent solution for the selected function is adapted to express the approxi¬ mate self-consistent solution for the selected value of the first selected external parameter and a selected value of the second selected external parameter by use of extrapolation based on the determined first and second self-consistent solutions together with the first and second values of the first selected external parameter, and further based on the de¬ termined fourth and fifth self-consistent solutions together with the first and second values of the second selected external parameter. Here, the means for expressing the approxi- mate self-consistent solution or change in the self-consistent solution may be adapted for expressing such solution by use of linear extrapolation.

The above described computer system for modelling a three-probe system may further comprise: means for determining a sixth self-consistent solution to the selected function for a third value of the second selected external parameter by use of self-consistent loop calculation, said third value of the second selected external parameter being different to the first and second values of the second selected external parameter. Here, the means for expressing the approximate self-consistent solution or change in the self-consistent solution for the selected function may be adapted to express the approximate self-consistent solution for the selected value of the first selected external parameter and the selected value of the second selected external parameter by use of extrapolation based on at least the deter¬ mined first, second and third self-consistent solutions together with the first, second and third values of the first selected external parameter, and further based on at least the de- termined fourth, fifth and sixth self-consistent solutions together with the first, second and third values of the second selected external parameter. Here, the means for expressing the approximate self-consistent solution or change in the self-consistent solution may be adapted for expressing such solution by use of second order extrapolation.

For the system of the invention provided for the three-probe system, the first value of the second selected external parameter may be equal to the first value of the first selected external parameter.

Also for the computer system of the present invention it is preferred that the selected func- tion is selected from the functions represented by: the effective one-electron potential en¬ ergy function, the effective one-electron Hamiltonian, and the electron density. Here, it is again preferred that the selected function is the effective one-electron potential energy function or the effective one-electron Hamiltonian and the self-consistent loop calculation is based on the Density Functional Theory, DFT, or the Hartree-Fock Theory, HF.

According to an embodiment of the invention, the computer may further comprise means for performing a self-consistent loop calculation based on a loop calculation including the steps of: a) selecting a value of the electron density for a selected region of the model of the nano-scale system, b) determining the effective one-electron potential energy function for the selected electron density and for a selected value of the external parameter, c) calculating a value for the electron density corresponding to the determined effec¬ tive one-electron potential energy function, d) comparing the selected value of the electron density with the calculated value of the electron density, and if the selected value and the calculated value of electron density are equal within a given numerical accuracy, then e) defining the solution to the effective one-electron potential energy function as the self-consistent solution to the effective one-electron potential energy function, and if not, then f) selecting a new value of the electron density and repeat steps b)- f) until the se¬ lected value and the calculated value of electron density are equal within said given nu¬ merical accuracy. Here, the means for performing the self-consistent loop calculation may be adapted to determine the self-consistent solution to the effective one-electron potential energy function for the probe or electrode regions of the system.

For embodiments wherein the means for performing the self-consistent loop calculation may be adapted to determine the self-consistent solution to the effective one-electron potential energy function for the probe or electrode regions of the system, it is also pre- ferred that the computer system further comprises means for determining Green's func¬ tions for each of the probe or electrode regions based on the corresponding determined self-consistent solution to the effective one-electron potential energy function.

For the computer system of the invention it is also within an embodiment that the selected function is the effective one-electron Hamiltonian for an interaction region of the system, and the means for determining a second self-consistent solution to the effective one- electron Hamiltonian of the interaction region of the system is adapted to perform said determination by including the step of calculating a corresponding self-consistent solution to the effective one-electron potential energy function for the interaction region at a given value of the first selected external parameter. Here, the means for determination of a sec¬ ond self-consistent solution to the effective one-electron Hamiltonian is adapted to per¬ form said determination based on a loop calculation including the steps of: aa) selecting a value of the electron density for the interaction region of the system, bb) determining the effective one-electron potential energy function for the selected electron density for a given value of the selected external parameter, cc) determining a solution to the effective one-electron Hamiltonian for the interaction region based on the in step b) determined effective one-electron potential energy function, dd) determining a solution to Green's function for the interaction region based on the in step c) determined solution to the effective one-electron Hamiltonian, ee) calculating a value for the electron density corresponding to the determined Green's function for the interaction region, ff) comparing the selected value of the electron density with the calculated value of the electron density, and if the selected value and the calculated value of electron density are equal within a given numerical accuracy, then gg) defining the solution to the effective one-electron Hamiltonian as the self-consis¬ tent solution to the effective one-electron Hamiltonian, and if not, then hh) selecting a new value of the electron density and repeat steps bb)- hh) until the selected value and the calculated value of electron density are equal within said given numerical accuracy.

Also the computer system of the invention covers an embodiment wherein the selected function is the effective one-electron Hamiltonian and the external parameter is a voltage bias across two probes of the system, wherein the means for determining a first and a second self-consistent solution is adapted to perform said determination for the effective one-electron Hamiltonian for selected first and second values, respectively, of the external voltage bias, and wherein the means for expressing an approximate self-consistent solu¬ tion by use of extrapolation analysis is adapted to obtain an extrapolation expression to an approximate self-consistent solution for the effective one-electron Hamiltonian when the external voltage bias is changed, said computer system further comprising: means for determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias using the obtained extrapolation expression, which expresses the approximate self-consistent solution or change in the self-consistent solution for the effective one-electron Hamiltonian. Here, the obtained extrapolation ex- pression may be a linear extrapolation expression. The means for determining the electri¬ cal current may be adapted to determine the electrical current for a given range of the external voltage bias and for a given voltage step in the external voltage bias, and the means for determining the electrical current may be adapted to perform said determina¬ tion using the following loop: aaa) determining the current for the lowest voltage within the given range of the exter¬ nal voltage bias, bbb) increasing the voltage bias by the given voltage step, ccc) determining the current for the new increased voltage bias, ddd) repeating steps bbb) and ccc) until the new increased voltage bias is larger than the highest voltage of the given range of the voltage bias.

It is also within an embodiment of the computer system of the invention that the system being modelled is a two-probe system and that the selected function is the effective one- electron Hamiltonian and the external parameter is a voltage bias across two probes of the system, and wherein the computer system further comprises: means for dividing a determined voltage range of the external voltage bias in at least a first and a second voltage range, means for determining for the first and second voltage ranges a maximum and a minimum self-consistent solution to the effective one-electron Hamiltonian corresponding to the maximum and minimum values of said voltage ranges, means for obtaining a first extrapolation expression to the approximate self-consis¬ tent solution for the effective one-electron Hamiltonian when the external voltage bias is changed, said first extrapolation expression being based on the determined maximum and minimum self-consistent solutions for the first voltage range and the maximum and mini- mum voltage values of the first voltage range, means for obtaining a second extrapolation expression to the approximate self- consistent solution for the effective one-electron Hamiltonian when the external voltage bias is changed, said second extrapolation expression being based on the determined maximum and minimum self-consistent solutions for the second voltage range and the maximum and minimum voltage values of the second voltage range, means for determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias within the voltage range given by the minimum and maximum voltage of the first voltage range using the obtained first extrapolation expression, and means for determining the electrical current between the two probes of the system for a number of different values of the applied voltage bias within the voltage range given by the minimum and maximum voltage of the second voltage range using the obtained second extrapolation expression. Here, the obtained first and second extrapolation ex¬ pressions may be first and second linear extrapolation expressions, respectively.

Other objects, features and advantages of the present invention will be more readily ap¬ parent from the detailed description of the preferred embodiments set forth below, taken in conjunction with the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

Fig. 1 is a flowchart (flowchart 1) illustrating the computational steps in a self-consistent loop of the Density Functional Theory. Fig. 2 illustrates a Benzene-Di-Thiol molecule coupled with two Gold (111 ) surfaces, where the gold surfaces are coupled to an external voltage source, and the electrodes have different chemical potentials μLax\ύμR .

Fig. 3 shows the self-consistent electron density of a carbon nano-tube coupled with a gold surface, where, when being outside the interaction region, the electron density is given by the bulk density of the electrodes.

Fig. 4a shows equivalent real axis (R) and complex contours (C) that can be used for the integral of Green's function G1 (z) .

Fig. 4b shows the variation of the spectral density (-ImG7(Z) ) along contour C (dashed) π of Fig. 4a and along the real axis R (solid) of Fig. 4a.

Fig. 5 is a flowchart (flowchart 2) showing steps required to calculate a self-consistent effective potential energy function of a two-probe system with applied voltage U using the Green's function approach, and where from the self-consistent effective one-electron po¬ tential energy function the electrical current / can be calculated.

Fig. 6 is a flowchart (flowchart 3) showing steps required for a self-consistent calculation of the I-U characteristics of a two-probe system.

Fig. 7a shows the self-consistent effective one-electron potential energy function of the system illustrated in Fig. 2 and calculated for different values of the applied voltage.

Fig. 7b shows the self-consistent effective one-electron potential energy function rescaled with the applied voltage.

Fig. 8 is a flowchart (flowchart 4) showing steps involved when using a linear extrapolation expression according to an embodiment of the invention to calculate the current-voltage characteristics, I-U.

Fig. 9 is a flowchart (flowchart 5) showing how an interpolation formula or linear extrapola¬ tion expression according to an embodiment of the invention can be used to calculate the I-U characteristics. Fig. 10 shows the result of a calculation of the / (U) characteristics of the system illus¬ trated in Fig. 2, with the line denoted "SCF" showing the result obtained with a full self- consistent calculation, while the line denoted "1. order" is showing the result obtained us- ing the scheme illustrated in Fig. 8, and the line denoted "2. order" is a second order ap¬ proximation.

Fig. 11 is a flowchart (flowchart 6) illustrating the use of an adaptive grid algorithm accord¬ ing to an embodiment of the invention for calculating the current voltage characteristics, /- U.

Fig. 12 is a flowchart (flowchart 7) being a recursive flowchart used by flowchart 6 of Fig. 11.

Fig. 13 is a flowchart (flowchart 8) being a recursive flowchart used by flowchart 6 of Fig. 11.

DETAILED DESCRIPTION OF THE INVENTION

Background theory

The purpose of atomic-scale modelling is to calculate the properties of molecules and materials from a description of the individual atoms in the systems. An atom consists of an ion core with charge Z, and an equal number of electrons that compensate this charge. We will use J^ , Zμ for the position and charge of the ions, where// = 1...N label the

ions, and N is the number of ions. The positions of the electrons are given by η, i = l...n , and n is the number of electrons.

Usually it is a good approximation to treat the ions as classical particles. The potential energy of the ions, V(R1,...,RN) , depends on the energy of the electronic system, E0 , through

where e is the electron charge. The electrons must be described as quantum particles, and the calculation of the electron energy requires that we solve the many-body Schrδdinger wave equation

H T(F15 ^ 5FJ = E0T(F1, ...,Fn) , Eq. 2

In Eq. 2, H is the many-body Ηamiltonian and Ψ the many-body wavefunction of the electrons. The "hat" over the many-body Ηamiltonian, H , symbolizes that the quantity is a quantum mechanical operator. The first term in Eq. 3 is the kinetic energy of the electrons, with ft = h/2π where h is Planck's constant, m the electron mass and V;the gradient

with respect to F1 . The second term is the electrostatic electron-ion attraction, and the last term is the electrostatic electron-electron repulsion.

The last term couples different electrons, and gives rise to a correlated motion between the electrons. Due to this complication an exact solution of the many-body Schrόdinger equation is only possible for systems with a single electron. Thus, approximations are required that can reduce the many-body Schrδdinger equation into a practical solvable model. A number of successful approaches have used an effective one-electron Ηamilto¬ nian to describe the electronic structure, and included the electron-electron interaction via an effective one-electron potential energy function in the one-electron Ηamiltonian.

Density functional method for electronic structure calculations

The invention can be used with electronic structure methods, which describe the electrons with an effective one-electron Ηamiltonian. DFT and HF theory are examples of such methods. In these methods the electrons are described as non-interacting particles mov¬ ing in an effective one-electron potential setup by the other electrons. The effective one- electron potential depends on the average position of the other electrons, and needs to be determined self consistently.

In Eq. 4 the term the effective one-

electron potential energy function andHlel is the one-electron Hamiltonian. The effective one-electron potential energy function depends on the electron density n. The kinetic en¬ ergy is given by a simple differential operator, and therefore independent of the density. This means that the effective one-electron potential energy function and the Hamiltonian has the same variation as function of the density, and when we are interested in determin¬ ing the self-consistent change of the effective one-electron potential energy function it is equivalent to specifying the self-consistent change of the Hamiltonian. Furthermore, for the self-consistent solution there is a one to one relation between the electron density and the effective one-electron potential energy function, thus specifying the self-consistent electron density, Hamiltonian or effective one-electron potential are equivalent.

In DFT the effective one-electron potential energy function is given by Veff [n] = Vim + Vxc[n]+ VH [n] . Eq . 5

The first term is the ion potential energy function which is given by the electrostatic poten¬ tial energy from the ion cores

and therefore independent of n. The second term is the exchange-correlation potential energy function

Vxc (F) = f(n(r), Vn(r ), Ψn(r )) , Eq . 7

which is a local function of the density and its gradients. The third term is the Hartree po¬ tential energy function, which is the electrostatic potential energy from the electron density and it can be calculated from the Poisson's equation

V2VH(r) = -4πen(r) . Eq. 8

Poisson's equation is a second-order differential equation and a boundary condition is required in order to fix the solution. For isolated systems the boundary condition is that the potential energy function asymptotically goes to zero, and in periodic systems the bound¬ ary condition is that the potential energy function is periodic. For such boundary condi- tions the solution of the Poisson's equation is straight-forward, and VH can be obtained from standard numerical software packages. For systems with an external voltage U, we solve the Hartree potential in separate parts of the system. This situation is discussed in more detail on page 26.

Thus from the density, we can obtain the effective one-electron potential energy function and thereby the Hamiltonian. The next step is to calculate the electron density from the Hamiltonian. It can be obtained by summing all occupied one-electron eigenstates.

Huι ψa (r) = εaΨa <J) , Eq. 9

For systems with a single chemical potential the occupied eigenstates are the states with an energy below the chemical potential. For systems with an applied external voltage U there are two chemical potentials and the situation more complicated. This situation is described on page 25.

The flowchart in Fig. 1 illustrates the self-consistent loop required to solve the equations. The system is defined by the position of the atoms Rμ (ionic coordinates), and external parameters like applied voltage U, temperature T, and pressure P, 102. Initially we make an arbitrary guess of the electron density of the system, 104. From the density we can construct the effective one-electron potential energy function using Eq. 5, 106. The effec¬ tive one-electron potential energy function defines the Hamiltonian through Eq. 4, 108. From the Hamiltonian we can calculate the electron density of the system by summing all occupied one-electron eigenstates as shown in Eq. 9, 10. If the new density is equal (within a specified numerical accuracy) to the density used to construct the effective one- electron potential energy function, 112, the self-consistent solution is obtained, 114, and we stop, 116. If the input and output electron densities are different, we make a new im- proved guess based on the previously calculated electron densities. In the simplest ver¬ sion the new guess is obtained from a linear mixing of the two electron densities, with a mixing parameter/? , 110. Application of DFT to closed and periodic systems

We will first show how Eq. 9 is most commonly solved for periodic and closed systems. A closed system is a system with a finite number of atoms. A periodic system is a system with an infinite number of atoms arranged in a periodic structure. For these systems, Eq. 9 is usually transformed into a matrix eigenvalue problem that can be solved with standard linear algebra packages. The transformation is obtained by writing the wave functions, ψα ,

as a linear combination of basis functions, ψα (F) = ∑α" φt (F) . Many different choices

exist for the basis functions, φt , some of the most common are plane-waves or atom- based functions with shapes resembling the atomic wave functions. Using the basis func¬ tions, Eq. 9 is transformed into

Sy = {φ, Ψj ) , Eq- 13

*(/=) = ∑ ∑(αJα)αy:{r)φJ{7) . Eq. 14 'J εα<μ

The symbol H denotes the Ηamiltonian matrix, and S the overlap matrix. The "bar" above the letters indicates that the quantities are matrixes.

For a molecular system the Ηamiltonian matrix is finite and it can be diagonalized with standard linear algebra packages. For a periodic structure it is only necessary to model the part of the system, which when repeated, generates the entire structure. Thus, again the Ηamiltonian matrix will be finite and the solution will be straight forward.

Application of DFT to open systems with an applied voltage

The application area of the invention is to systems where two (or more) semi-infinite elec¬ trodes are coupling with a nano-scale interaction region. We call such systems two-probe systems. The nano-scale interaction region can exchange particles with the electrodes and the two-probe systems are therefore open quantum mechanical systems. The left and right electrodes are electron reservoirs with definite chemical potentials, μL and μΛ . The difference between the chemical potentials, μL ~ μR = eU , Eq. 15

5 defines the voltage bias, U , applied to the system. For open systems the Hamiltonian matrix is infinite and the simple diagonalization technique in Eq. 11 for obtaining the one- electron eigenstates cannot be applied. Instead we will determine the electron density using the non-equilibrium Green's function formalism described in the following sections. Examples of two-probe systems are illustrated in Figs. 2 and 3. The system in Fig. 2 con- 10 sists of two semi-infinite gold electrodes coupling with a Phenyl Di-Thiol molecule. The interaction region 22 consists of the molecule and the first two layers of the electrodes. Regions 21 , 23 show the left and right electrode regions. Regions 24, 26 show the occu¬ pation of the one-electron levels within the electrodes; due to the applied voltage the chemical potential of the right electrode 26 is higher than for the left electrode 24. 15 Fig. 3 shows a semi-infinite carbon nano-tube coupling with a semi-infinite gold wire. The interaction region 32 is given by the nano-tube apex and the first layers of the gold wire. The left electrode 31 consists of a semi-infinite gold wire, and the right electrode 33 con¬ sists of a semi-infinite carbon nano-tube. The electron densities in the left electrode region 20 34 and in the right electrode region 36 are obtained from self-consistent bulk calculations. These densities seamlessly match the self-consistently calculated two probe density of the interaction region 35.

The screening approximation 25 The first step is to transform the open system into three subsystems that can be solved independently. Fig. 3a shows a carbon nano-tube coupled with a gold wire. The gold wire and the carbon nano-tube are metallic. Because of the metallic nature of the semi-infinite electrodes, the perturbation due to the interaction region only propagates a few Angstrøm 30 into the electrodes. This is illustrated in Fig. 3b, which shows the electron density. We see that when we move a few atomic distances away from the nano-tube to gold contact point, the electron density is periodic and resembles the bulk electron density. Thus, we can divide the electron density and the effective one-electron potential energy function into an interaction region and electrode regions, where the value in the electrode region is similar to the electrode bulk value. This is called the screening approximation.

Since the effective one-electron potential energy function is a local operator, the Hamilto- nian operator can also be separated into electrode and interaction region. Thus, if we ex¬ pand the Hamiltonian operator in a basis set with finite range, the Hamiltonian matrix can be separated into

where H1x , H11 and Hm denotes the Hamiltonian matrix of the left electrode, interaction

region, and right electrode, respectively, and H11 and H1R are the matrix elements in¬ volving the interaction region and the electrodes. Note that the size of the interaction re¬ gion is such that there are no couplings between the left and right electrode, i.e.

HLR = HRL = ° ■

Calculating the electron density using Green's functions

We will now show how the electron density is obtained within the Green's function formal¬ ism. For this purpose we introduce the spectral-density, D(ε) , and the electron density

operator N . The spectral density is the energy resolved electron density, and the total electron density is obtained by integrating the spectral density over all energies

D(ε) = δ(ε -H) , Eq. 17

In Eq. 17, the function δ(x) \s Dirac's delta function. The (retarded) Green's function is defined by

where δ+ is an infinitesimal positive number and / is the complex base. The Green's func¬ tion is related to the spectral density through D(ε) = —ϊmG(ε) , Eq. 20 π

where ImG is the imaginary part of G . Expanding the operators in basis functions, we

transform Eq. 19 into a matrix equation

From the Green's function we can obtain the spectral density matrix

£> O) = -ImG(^), Eq. 22 π

and thus the electron density

n(r) = ∑NiJφi(r)φj(r) , Eq. 23 u

The calculation of the electron density is now reduced to the matrix inversion in Eq. 21 , and the energy integral in Eq. 24. However, we have an open system and the matrix in Eq. 21 is therefore infinite. Due to the screening approximation we only need to calculate the electron density in the interaction region since in the electrode regions we can use the bulk electron density. From Eq. 24 we see that since our basis functions are localized, we only need to calculate the Green's function matrix of the interaction region and a few lay¬ ers of the electrodes.

Including the electrode region through a self energy term

In this section we will show how the Green's function matrix of the interaction region, Gn ,

can be calculated by inverting a matrix with the same size. To obtain this result we will

use perturbation theory in the coupling elements Hu(ε) = H11 - εSu

andHjyiε) = H10 - εSN . The unperturbed Green's functions, G0 , is calculated by setting

Hu = Hm = 0 and using that in this case Eq. 21 is block diagonal

GLQL{ε) = [{ε + iδ,)SLL -ΗlL]\ Eq. 25 Gf1 (S) = [(ε + iδ+ )S J1 -H11Y 1 Eq. 26

GfR(ε) = [(s + iδ+)SRR -HRRY . Eq. 27

Putting back the perturbation Hn and H ^ we find the Green's function from the Dyson's equation

G11 (S) = Gf1 (ε) + Gf1 (ε )\∑f, (ε) + If1 (ε)% (ε) , Eq . 28

?ni L (s) = H1L(s)GfL(s)Hn(s) , Eq. 29

where the terms ∑fj(ε) and ∑f7(£") are called the selfenergies of the electrodes. Rear¬ ranging the terms in the Dyson's equation, we arrive at

Gπ(ε) = [(S + IO+)S11 -H11 -If1(S) -If1(S)Y . Eq. 31

Calculation of the electrode Green's function

In order to determine the self energies we need to calculate the unperturbed Green's function, Gf1 , of the electrodes. Since, the Hamiltonian of the electrodes is semi-infinite, the Green's function cannot be obtained by simple matrix inversion. However, in cases where the electrode Hamiltonian is periodic, there exist very efficient algorithms for calcu¬ lating the electrodes Green's function. Below we will describe one of them. We will write the electrode Hamiltonian as periodic blocks, Hhh = H1^ = • •• , where the size of each block is such that only neighbouring blocks interact, i.e.

The Hamiltonian of each block, H1^ and the coupling matrix, HLχLι , are obtained from a bulk calculation of the electrode system. Using recursion, we build up a series of approxi¬ mations for the Green's function

: Eq. 36

In Eq. 33, 34, 35 the quantity G^(ε) is the n'order approximation to the Green's func¬

tion. The error, decreases as 1/n where n is the number of steps. Due to this poor convergence usually more than 1000 steps are required to obtain rea¬ sonable accuracy with this algorithm. The Green's function can be obtained in fewer steps by using a variant of the method described in Lopez-Sancho, J. Phys. F 14, 1205 (1984). With this variant of the algorithm only a few steps are needed to calculate the electrode Green's function, and the computational resources required for this part is usually negligi- ble compared to the resources required for the calculation of G77 .

Integration of the spectral density using a complex contour

We now have all the ingredients required in Eq. 31 to obtain G77 and thereby the electron density matrix,

The Green's function is a rapidly varying function along the real axis, and for realistic sys¬ tems often an accurate determination of the integral requires more than 5000 energy points along the real axis. To find a more efficient method we use that the Green's func¬ tion is an analytical function, and we can do the integral along a contour in the complex plane. In the complex plane the Green's function is very smooth. This is illustrated in Fig. 4. In Fig. 4a we show two equivalent lines of integrations, the contour C and the real axis line R. Fig. 4b shows the variation of the spectral density along C (dashed) and along R (solid). The function varies much more rapidly along R, and substantially more points are needed along R than along C to obtain the same accuracy. Typically, the use of contour integration reduces the number of integration points by a factor 100.

The electron density for a twoprobe system with external voltage bias

We have so far used that the system has a single chemical potential, i.e. μL = μR . How- ever, if we apply an external voltage, U , the two electrodes will have different chemical potentials linked through Eq. 15. Fig. 2 illustrates the system set up. The energy axis can be divided into two regions, the energy range below both chemical potentials we call the equilibrium region, and the energy range between the two chemical potentials we call the non-equilibrium region or voltage window. We will divide the electron density matrix into two parts,

Ny = N? + N™\ Eq. 38

where Njj? is the electron density matrix of the electrons with energies in the equilibrium

region, and N1"69 the electron density matrix of the electrons with energies in the non-

equilibrium region. We may say that N1"6" is the additional density due to the external voltage U.

Nygcan be calculated with the approach described in the previous sections, thus

- \ % — N? = - J Im Gy {ε)dε , Eq . 39

where we have assumed that//z < μR .

In the non-equilibrium region electrons are only injected from the right reservoir. Thus we need to divide the spectral density matrix into electron states originating from the left or right electrode, and only add the right electrode electron density. This division of the elec¬ tron density is accomplished in non-equilibrium Green's function theory, and we may write (ε)lm∑fI(ε) GjI(ε)dε. Eq. 40

The foundation of this equation can be found in Haug and A. P. Jauho, Quantum kinetics in transport and optics of semiconductors, (Springer-Verlag, Berlin, 1996) or Brandbyge Phys. Rev. B 65, 165401 (2002). Thus, we now have a description for how to calculate the electron density of the two probe system, including the situation with an external volt¬ age applied to the system.

Calculating the effective one-electron potential energy function in a two-probe system

In the previous sections we showed how to calculate the electron density from the Hamil- tonian using the Green's function approach. To complete the self-consistent cycle we need to calculate the Hamiltonian from the electron density, which means calculating the effective one-electron potential energy function, V® [n]. Within DFT the effective one- electron potential energy function is given by Eq. 5. For the two-probe system we need to solve Poisson's equation, Eq. 8, for the interaction region and the electrode regions sepa¬ rately. The Hartree potential energy function of the electrodes is obtained with the same approach as used for periodic systems, in this case the repeated structure is the electrode cell used to define H1^ in Eq. 32 and the corresponding cell for the right electrode HRΛ . These electrode Hartree potential energy functions now supply boundary conditions for the Hartree potential energy function of the interaction region. However, the electrodes are bulk systems and this means that we can add an arbitrary constant to their Hartree potential energy function and still obtain a valid solution. To fix this arbitrary constant we relate each electrode Hartree potential energy function to the chemical potential of the electrode, and use Eq. 15 to relate the left and right chemical potential. Thus, we now have fixed the Hartree potentials in the electrodes and they define the boundary condition of the Poisson's equation in the central region along the z direction. In the x and y direc¬ tion we will use periodic boundary conditions. With these boundary conditions the Hartree potential energy function of the interaction region can be obtained by a multigrid ap¬ proach, as described in Taylor, Phys. Rev. B 63, 245407 (2001). Electron transport coefficients and currents obtained from the Green's function

After finishing the self-consistent cycle we can calculate the transport properties of the system. The non-linear current through the contact, /, is obtained as

HU) = G0 Eq. 41

Q" where G0 = 2 — is the conduction quantum. The foundation of this equation is described h in H. Haug, Quantum kinetics in transport and optics of semiconductors, (Springer- Verlag, Berlin, 1996).

The self-consistent algorithm for the two-probe system

Fig. 5 shows required steps for a two-probe calculation of the electrical current from the left to the right electrode through a nano-scale device due to an applied voltage between the left and right electrode as described in Eq. 15. Initially we define the system by speci- fying the ionic positions, and the external parameters like the applied voltage and tem¬ perature, 202. Next we use the screening approximation to separate the system geometry into interaction and electrode regions, 204. The electron density and the effective one- electron potential energy function should approach their bulk value in the electrode region. Usually this will be the case around atoms in the third layer of a metallic surface, and it is therefore sufficient to include the first two layers of metallic surfaces within the interaction region. We calculate the self-consistent effective one-electron potential energy function for the isolated electrode regions using the flowchart in Fig. 1 , 206. From the self-consistent effective one-electron potential energy function we construct the electrode Greens func¬ tions, using Eq. 4, 12, 33-36, and the electrode selfenergies using Eq. 29, 30, 208. These initial calculations are now used as input to the two-probe calculation. Thus, we have cal¬ culated the self-consistent density of the electrode regions and only need to calculate the self-consistent density of the interaction region. Starting with an initial guess of the elec¬ tron density for the interaction region, 210, we perform a self-consistent loop similar to the flowchart in Fig. 1. First we calculate the effective one-electron potential energy function of the interaction region using Eq. 5-8, 212. From the effective one-electron potential en¬ ergy function we can obtain the Hamiltonian using Eq. 4, 12 and the Green's function through Eq. 31 , 214. From the Green's function we can calculate the electron density us- ing Eq. 23, 38, 39, 40, and thereby close the self-consistent cycle, 218. If the new electron density is different (within a specified numerical accuracy) from the electron density used to construct the effective one-electron potential energy function, 220, we make a new im¬ proved guess based on the previously calculated densities. In the simplest version the new guess is obtained from a linear mixing of the two densities, with a mixing parame¬ ter/? , 216. If the input and output densities are equal, we have obtained the self-consis¬ tent value of the electron density and thereby also the effective one-electron potential en¬ ergy function, Hamiltonian and Green's function, 222. From this Green's function we can calculate the current using Eq. 41 , 224. After the calculation of the current the algorithm stops, 226.

The procedure has been implemented in the TranSIESTA and McDCAL software. Further description of these softwares and the implementation details can be found in Brandbyge Phys. Rev. B 65, 165401 (2002), and Taylor Phys. Rev. B 63, 245407 (2001 ). To obtain the current-voltage characteristics, I-U curve, of a nano-scale device, we need to perform a self-consistent calculation for each voltage U. This is illustrated in the flowchart of Fig 6. Input system geometry and the voltage interval U0 , U1 and step size AU , 302. Set start¬

ing voltage to CT0 , 304. Follow the steps in flowchart 2 of Fig. 5 to perform a self-consistent calculation of the effective one-electron potential energy function at voltage U, and use the self-consistent potential energy function to calculate the current, 306. Increase the voltage with the step size, 308, if the new voltage is within the specified voltage interval, then perform a new self-consistent calculation, 310, else stop, 312. i

Example: Calculation of the I-U characteristics of DTP coupled with gold surfaces

We will now present results for the calculation of the I-U characteristics of the geometry illustrated in Fig. 2 using the TranSIESTA software. The calculation follows flowchart 3 of Fig. 6, and the points in Fig. 9 show the result of the calculation. A similar I-U characteris¬ tic was obtained in Stokbro Computational Materials Science 27, 151 (2003).

In Fig. 7 we show the change in the self-consistent effective one-electron potential energy function due to the applied voltage. The value of the effective one-electron potential en¬ ergy function is shown along a line starting in the left electrode, going through the center of the two sulphur atoms of the DTB molecule and ending in the right electrode. In the right electrode the effective one-electron potential energy function is shifted down due to the applied voltage. The main feature is that the effective one-electron potential energy function is flat in the electrode regions, and the main voltage drop is taken place within the molecular region.

The curves in Fig. 7a all have similar shapes. In Fig. 7b we have rescaled the curves with the applied voltage, and we observe that the rescaled effective one-electron potential en¬ ergy functions are nearly identical. This observation forms a basis for the invention as it shows that the self-consistent change in the effective one-electron potential energy func¬ tion has a simple variation with the applied voltage.

Linear interpolation using two voltage points

In one version of the algorithm, the effective one-electron potential energy function is cal¬ culated at zero voltage, U0 and for a small finite voltage, t/Δ . These data are now used to extrapolate to a general voltage. The effective one-electron potential energy function for the general voltage, U, is obtained by simple linear extrapolation

The Hamiltonian is related to the effective one-electron potential energy function by

H = -- V2 + Veff . Eq. 43 2m

This means that the same scaling relation applies to the Hamiltonian. Thus, the Hamilto¬ nian at a general voltage can be approximated by

where H5CF[t/0] and H5CF[t/Δ] are the self-consistent Hamiltonian at U0 and UL .

In most electronic structure methods the Hamiltonian is expanded in a basis settø. } , and represented by the matrix

HV = U H ΦJ) . Eq. 45 In this case the linear interpolation formula is applied to the Hamiltonian matrix elements.

From the Hamiltonian we can calculate all properties of the system, including the electrical current due to the applied voltage. The electrical current is obtained by first calculating the Green's function using Eq. 31 and from the Green's function calculate the current using Eq. 41. We may combine Eq. 44, 31 and 41 and write it as a mapping, M, that takes H5CF[t/0] , , U, and returns the current, /, at voltage U. We write the map¬ ping as

I(U) := M(U,HSCF [U0\HSCF[UA]), Eq. 46

The calculation of the I-U characteristics using the interpolation formula is summarized by flowchart 4 in Fig. 8. Input system geometry and the voltage interval Uλ , U2 , step size AU , and voltages U0 UA where we will calculate the self-consistent Hamiitonians that are used for the interpolation, 402. Use flowchart 2 of Fig. 5 to calculate the self-consis¬ tent effective one-electron potential energy function and Hamiltonian for voltage U0 , 404.

Self-consistent calculation for voltage UA , 406. Use flowchart 5 of Fig. 9 to calculate the /- U curve for the voltage interval U1JJ2 using Eq. 46 with the self-consistent results at U0

and UA io obtain an approximation for the current, 408. Stop, 410.

The calculation of the I-U curve follows flowchart 5 of Fig 9. Input voltage interval U1 , U2 , step size AU and the self-consistent Hamiltonian for two voltages, U0 and UA , 502. Set

starting voltage tot/ := U1 , 504. Use Eq. 46 with the self-consistent results at U0 and

U A to obtain an approximation for the current at U, 506. Increase the voltage with the step size, 508, if the new voltage is within the specified voltage interval, then continue calculat¬ ing the I-U curve, 510, else stop, 512.

Typical parameters for the calculation will be to select U0 = 0 Volt and UA = 0.4 Volt. It is most computationally efficient to choose a relative low value of the voltage, since the self- consistent calculation is more computationally demanding at a high voltage due to the calculation of the non-equilibrium density, Eq. 40, which involves an integral where the number of points is proportional to the size of the voltage.

Typical values for the range of the voltage in the I-U curve will be U1 = - 2.0 Volt and U2 =2.0 Volt. At higher voltages the electric field will be very high for a small nano-scale device, and such voltages are difficult to measure experimentally due to electrical break¬ down of the device.

In Fig. 10 we compare the result of calculating the current using the formula in Eq. 46 with the full self-consistent solution. The line denoted "1. order" shows the result obtained with Eq. 46, while the line denoted "SCF" shows the result obtained with the self-consistent calculation. We see that the results obtained with Eq. 46 are in excellent agreement with the full self-consistent calculation for V < 2.0 Volt, even though only calculations at V=O.0 Volt and V=O.4 Volt were used for the calculation.

Adaptive grid method for calculating the I-U characteristics.

In the previous section we used a two point interpolation formula to extrapolate the Hamil- tonian to a general voltage using the self-consistent Hamiltonian at two voltages t/oand C/Δ . We will now propose a systematic method to improve this scheme. The method is based on performing additional self-consistent calculations at selected voltage points, and using the self-consistent Hamiltonian at these voltage points to make im¬ proved interpolation formulas. With this method a series of I-U curves are produced that converges towards the self-consistently calculated I-U characteristics.

The target is to calculate the I-U characteristics in the interval [∑/15t/2]. Flowchart 6 in Fig. 11 shows the steps involved in the calculation. The initial steps are similar to flowchart 4 of Fig. 8; however, in this new algorithm we will improve the approximation by performing additional self-consistent calculations, where the new voltage points may be selected by the algorithm shown in flowcharts 7 and 8 of Figs. 12 and 13. Input system geometry and the voltage interval U1 , U2 , step size AU , and interpolation voltages U0 , UA , 602. Use flowchart 2 of Fig. 5 to calculate the self-consistent effective one-electron potential energy function and Hamiltonian for voltage U0 , 604. Self-consistent calculation for volt¬

age C/Δ , 606. Use flowchart 8 of Fig. 13 to calculate the I-U curve for the voltage interval Uλ,U0 using Eq. 46 with the self-consistent results at U0 and C/Δ to obtain an approxima- tion for the current, 608. Use flowchart 7 of Fig. 12 to calculate the I-U curve for the volt¬ age interval U0,U2 using Eq. 46 with the self-consistent results at U0 and t/Δ to obtain an approximation for the current, 610. Stop 612.

Flowcharts 7 and 8 of Figs. 12 and 13 show the algorithms for subdivision of the interval. The interval is subdivided until interpolated and self-consistent calculated currents agree within a specified accuracy, which we denoted . Flowchart 7 and 8 are similar except that flowchart 7 assumes the self-consistent Hamiltonian is known for the lowest voltage U A of the voltage interval where we request the I-U curve, while flowchart 8 assumes the self- consistent Hamiltonian is known for the highest voltage UB of the voltage interval. For flowchart 7, the input to the recursion step is the voltage interval UA , UB , and the self- consistent Hamiltonian at the endpoint U A and at an arbitrary voltage point Uc , 702.

Next we perform a self-consistent calculation at the highest voltage U B of the voltage interval, 704. We calculate the current from the interpolation formula, Eq. 46 and from the self-consistent Hamiltonian Eq. 31 , 41 , 706. If the interpolated current differs by more than δ from the self-consistent current, 708, we will further subdivide into intervals {UA,UM } and {UM ,UB } , where UM := (JJ Λ + UB)/2 , 714. The algorithm is recursively called with the interval ] , 716. For the interval {UM,UB } we know the Hamilto¬ nian at the last voltage point instead of for the first voltage point, and we use the slightly modified algorithm shown in flowchart 8, 718. The procedure is continued until the self- consistently calculated current for the new grid point agrees with the interpolated value within the prescribed accuracy δ . When the prescribed accuracy is obtained we can safely use Eq. 46 to calculate the I-U characteristics of the subinterval , 710. The recursive algorithm stops, 712.

The algorithm in flowchart 8 of Fig. 13 is a slight modification of the algorithm in flowchart 7 of Fig. 12, the only difference being that the input self-consistent Hamiltonian is calcu- lated at UB instead of U A . Here we just mention the differences in flowchart 8 when compared to flowchart 7. Input H5Ciτ [t/5] instead of HSCF ] , 802. Perform self-

consistent calculation at U A instead of at UB , 804. Calculate the current at U A , 806, compare currents calculated at U A , 808. The remainder of the algorithm is similar to the algorithm flowchart 7.

We note that in general this procedure will result in grid points unevenly distributed over the voltage window. The grid points will be most dense in the regions where the linear interpolation formula gives a poor description of the variation of the self-consistent poten- tial energy function. Thus the algorithm results in an adaptive formation of the grid points.

Using higher order approximations

For the methods described in the previous section the approximate solution was system- atically improved by performing additional self-consistent calculations. When more than two self-consistent calculations are performed it is possible to use higher order interpola¬ tion formulas. For instance, self-consistent calculations at U0 , Ux , and CZ2 ,can be com¬ bined to obtain a second order extrapolation formula

v£[u}:= VgF[Uo]+ (U-Uo)b + (U - U0)2 C Eq. 46b

for the effective potential, Similar second order extrapolation formulas can be used for the Hamiltonian, H[U] := H[U0]+ (U - U0)b + (U - U0)2 c Eq. 46e

b =H[Ux]l(Ux -U0) -c(Ux - U0) Eq. 46g

The line denoted "2. order" in figure 10 shows the result using a second order extrapola- tion formula obtained from self consistent calculations at 0.0 Volts, 0.4 Volts and 1.0 volts. The above can easily be generalized such that for n biases a (n-1 ) order extrapolation formula is used.

Generalization to multi-probe systems

The algorithm can be generalized to multi-probe systems, i.e. systems where there are more than two electrodes. Lets assume that we will include one additional electrode, then we can relate the chemical potential of this electrode, μ3 , to the chemical potential of the

left electrode through the applied voltage between the electrodes, UL3

μL ~ μ3 = eUL\ Eq. 47

We can now generalize Eq. 44 to a two-dimensional interpolation formula in the variables U13 and ULR , where the latter is the voltage difference between the left and the right electrode. It is convenient to choose U03 = U0R = U0 = 0 , since then we can use the

same self-consistent Hamiltonian for the U0 value in the interpolation formula. In this case

where U" , U^R are a small voltage increase in the left electrode-electrode 3 and left electrode-right electrode voltages, respectively. The self-consistent Hamiltonians HSCF [t/f ] are calculated for UL3 = JJf , ULR = 0 , and HSCF [t/f ] are calculated for

Q

Generalization to use electronic or ionic temperature

So far we have implicitly assumed that the electronic temperature is zero, since all inte¬ grals are written with fixed integration boundaries at the chemical potentials. To include a finite electronic temperature we must change the integrals in Eq. 18, 24, 37, 39, 40, 41 such that

where T is the temperature, /c the Boltzmanns constant, and f is the Fermi function

/W= -^. Eq.50

We can readily generalize this to use different electronic temperatures for the left and right electrode, by using different values of T in the Fermi function for the left and right elec¬ trode.

Those skilled in the art will appreciate that the invention is not limited by what has been particularly shown and described herein as numerous modifications and variations may be made to the preferred embodiment without departing from the spirit and scope of the in¬ vention.