Numerical Methods for Analysis of Charged Vacancy 

Diffusion in Dielectric Solids 

by John D. Clayton, Peter W. Chung, Michael A. Greenfield, 

and William D. Nothwang 

ARL-TR-4002 December 2006 

Approved for public release; distribution is unlimited. 
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents. 

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the 

use thereof. 

Destroy this report when it is no longer needed. Do not return it to the originator. 
Army Research Laboratory 

Aberdeen Proving Ground, MD 21005-5066 

ARL-TR-4002 December 2006 

Numerical Methods for Analysis of Charged Vacancy 

Diffusion in Dielectric Solids 

John D. Clayton, Michael A. Greenfield, and William D. Nothwang 

Weapons and Materials Research Directorate, ARL 

Peter W. Chung

Computational and Information Sciences Directorate, ARL 

Approved for public release; distribution is unlimited. 
ii

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

December 2006 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

October 2005–October 2006 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

Numerical Methods for Analysis of Charged Vacancy Diffusion in Dielectric 

Solids 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

FY06WMR15 

5e. TASK NUMBER 

6. AUTHOR(S) 

John D. Clayton, Peter W. Chung, Michael A. Greenfield, 

and William D. Nothwang 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Research Laboratory 

ATTN: AMSRD-ARL-WM-TA 

Aberdeen Proving Ground, MD 21005-5066 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-4002 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

A theory for charged vacancy diffusion in elastic dielectric materials is formulated and implemented numerically in a finite 

difference code. The governing equations consist of Maxwell’s equations of electrostatics coupled with kinetic equations for 

vacancy diffusion, with the chemical potential accounting for both mixing energy of vacancies and electrostatically-driven 

charge migration. A second-order accurate implicit scheme is used to solve Maxwell’s parabolic equations, while an explicit 

method is used to integrate the elliptic evolution equations for transient vacancy concentration. In addition to the theoretical 

background and numerical methodology, user documentation is included for the computer implementation, presently limited to 

one-dimensional analysis. Provided here are descriptions of the code structure, user instructions, and a representative 

application of the software for analysis of barium strontium titanate thin films containing charged oxygen vacancies. The 

source code is included in the appendix. 

15. SUBJECT TERMS

thin films, dielectrics, electromechanics, diffusion 

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

John Clayton 

a. REPORT 

UNCLASSIFIED 

b. ABSTRACT 

UNCLASSIFIED

c. THIS PAGE 

UNCLASSIFIED

17. LIMITATION 

OF ABSTRACT 

UL 

18. NUMBER 

OF PAGES 

40 19b. TELEPHONE NUMBER (Include area code) 

410-306-0975 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18
iii

Contents 

List of Figures iv

List of Tables iv

1. Introduction 1

2. Theory 2

2.1 Model Framework ...........................................................................................................2

2.2 One-Dimensional Model .................................................................................................5 

3. Numerical Methods 7

3.1 Electrostatics ...................................................................................................................7 

3.2 Transient Diffusion..........................................................................................................8 

4. Software Manual 10

4.1 Code Structure...............................................................................................................10 

4.2 Input...............................................................................................................................11 

4.3 Output............................................................................................................................13 

5. Sample Problem-BST Film 14

6. Conclusions 17

7. References 18

Appendix. Source Code 19

Distribution List 33 
iv

List of Figures 

Figure 1. One-dimensional problem domain. .................................................................................5 

Figure 2. Spatial discretization. ......................................................................................................7 

Figure 3. Flowchart for code execution. ......................................................................................11 

Figure 4. Vacancy concentration and electric field: 0 c = 1 ppm , V = 0 .........................................16 

List of Tables 

Table 1. Properties of BST film at 298 K. ....................................................................................14 
1

1. Introduction 

Dielectric materials with charged defects exhibit a variety of physical phenomena whose origins 

are not fully understood. Early continuum theories (Devonshire, 1954; Toupin, 1956) of the 

electromechanical behavior of dielectric media have been set forth, though these do not consider 

defects explicitly. However, defects such as vacancies have been the focus of considerable study 

(Lifshitz, 1963), particularly with regards to crystalline ceramics of current interest to the U.S. 

Army Research Laboratory. 

The unique aspect of the present study is consideration of electrically charged, as opposed to 

neutral, vacancies in dielectric solids. Defect concentrations and excess charges in ferroic 

ceramics may be adjusted during processing via heat treatments and/or addition of doping 

chemicals (Cole et al., 2003). Charged point defects have been identified as a major factor 

affecting the reliability of ferroelectric devices (Damjanovic, 1998), including gate dielectric 

semiconductors, particularly those of thin film geometry (Buchanan, 1999). 

Here, a general modeling framework is constructed for elastic dielectric semiconductors with 

mobile charged point vacancies. This framework combines the physics of continuum elasticity, 

electrostatics, mass diffusion, and charged defect kinetics. Changes in surface morphology due 

to the boundary flux of charged vacancies are captured, extending a previous theory of one of the 

co-authors on neutral vacancy kinetics (Grinfeld and Hazzledine, 1997). 

The theory is implemented numerically in a finite difference code (Hoffman, 1992) enabling 

simultaneous solution of the elliptic equations of electrostatics of dielectrics and the transient 

parabolic equations of charged diffusion. The analysis is limited to a single spatial dimension. 

The time duration of the problem is decomposed into a sequence of steps. A second-order 

accurate fully implicit scheme is invoked to solve Maxwell's equations in each step, while a fully 

explicit scheme is used to integrate the transient vacancy concentration. The spatial domain and 

grid spacing are updated when the surface flux of concentration is nonzero, as vacancies exiting 

the domain influence its instantaneous dimensions. 

Documentation is presented for the computer implementation. Included here are descriptions of 

the source code structure, user instructions, and representative input files for the software, the 

latter specifically for analysis of barium strontium titanate (Ba1-xSrxTiO3) (BST) thin films 

containing charged oxygen vacancies. The source code is given in the appendix. 

In the notation that follows, the Einstein summation convention is used on repeated lower-case 

indices, unless indicated otherwise. Cartesian spatial coordinate indices span three dimensions 

and are written in Roman font, while curvilinear surface coordinate indices span two dimensions 
2

and are written in Greek font. Subscripted commas denote covariant differentiation. Capitalized 

subscripts are often used in physical constants and are not summed. In the description of 

numerical methods, subscripts are often used for node numbers, and superscripts for time 

increments. 

2. Theory 

The governing equations, thermodynamic framework, and constitutive relations are first 

presented, with limited derivations, in three-dimensional (3-D) form in section 2.1. The onedimensional (1-D) equations that are solved numerically follow in section 2.2. 

2.1 Model Framework 

The local electrostatic behavior of dielectric continua is dictated by Maxwell's equations 

(Stratton, 1941): 

, ˆ i Di = ρ , (1) 

Ei i = −φ, , (2) 

and 

0

i ii D EP = ε + , (3) 

where D is the electric displacement, ρˆ is the charge density, E is the electric field, φ is the 

electrostatic potential, and 0 ε is the permittivity of free space. The polarization vector P is 

defined only within the material and vanishes in free space. Local mechanical equilibrium and 

mass conservation are ensured by 

, 0 ij σ j = (4) 

and 

, 0 i

i ρ ρ  + u = , (5) 

where σ is the symmetric Cauchy stress, ρ is the mass density of the deformed solid with 

vacancies, and u is the displacement, which itself encompasses expansion or contraction due to 

vacancies within the bulk material. Small displacements are assumed henceforth in the present 

theory. The balance of energy and the dissipation inequality may be written in global rate form 

as 

Φ Ω =W − ≥ 0  , (6) 

with Φ the dissipation, W the rate of external work, and Ω the system energy. For a dielectric 

solid containing mobile charged vacancies, the external work and energy of the system are 
3

ˆ ˆi i

i i

d W ds t u ds Q n ds dt

= +− σφ μ ∫∫∫  (7) 

and 

( ) ( ) 2

0 , 1/2 ˆ i

d d U dv dv ds

dt dt Ω ρ θη ε φ σφ =−+ + ∫∫ ∫    , (8) 

where s and v are the surface and volume of the system with unit normal n, σˆ is the surface 

charge density, t

ˆ is the traction vector, μ is the chemical potential for vacancy diffusion, Q is 

the vacancy flux, U is the internal energy per unit mass of the substance, and γ is the surface 

energy density. Heat conduction is not considered explicitly here, as is evident from equations 

6–8. The Helmholtz free energy ψ , specific entropy η , and temperature θ are related by the 

usual thermodynamic relationship 

ψ =U −ηθ . (9) 

The following constitutive assumptions are made regarding the free energy and charge density: 

( , , ,, ) i

i j ψ =ψ θξ u P , (10) 

and 

ρˆ = ezξ , (11) 

where ξ is the number of vacancies per unit volume. Assumption (equation 10) suggests a free 

energy dependence on mechanical strain, polarization, temperature, and vacancy concentration. 

Equation 11 denotes that the charge density is proportionate to the vacancy concentration, with e

and z the charge of an electron and the valence contribution of each defect, respectively. 

Substituting equations 9–11 into equation 6 and making use of the balance laws in equations 1–5 

and the divergence theorem, and restricting 0 φ = on s, the following thermodynamically 

admissible bulk constitutive relations may be derived: 

i E P i = ρ ψ∂ ∂ , (12) 

,

ij

i j σ ρψ = ∂ ∂u , (13) 

η = −∂ ∂ ψ θ , (14) 

,

i Qi ξ = −  , (15) 

,

i ij Q d = − μ j , (16) 

and 

μ = ρψ ξ φ ∂ ∂+ ez , (17) 
4

where the diffusivity tensor d is symmetric and positive definite. Let χ =αξ denote the volume 

fraction of vacancies, with α is a scalar conversion factor. Mass conservation requires that 

(1 ) ( , ) i

i a Qn qα − =− + χ α α  , (18) 

where a is the surface velocity of the substance with unit normal n and q is the flux tangential to 

the surface of charged vacancies. Conservation of mechanical work leads to 

, ˆi

ij j t un =∂ ∂ ρ ψ . (19) 

By requiring that the volumetric energy of the system Ω remain constant or decrease with time 

in the absence of mechanical working due to external stresses, the following relations for the 

in-plane and normal surface fluxes are formulated: 

(( ) ( )), q A ˆ / 1 α αβ

β

= +− ρψ ρφχ χ , (20) 

and 

() () ( ( )) 1 1 1 ˆ i χ Q n ez i β ρψ ρφ χ ρ ψ χ α φ − − = + + − ∂ ∂+ , (21) 

where A is the symmetric positive definite surface diffusivity tensor and β is a positive scalar 

characterizing the resistance of the surface to penetration by vacancies. Implicit in equations 4, 

13, and 19 is the vanishing of Maxwell's stress tensor (Toupin, 1956), meaning that contributions 

from terms of second order in the electric field and polarization are neglected in the mechanical 

equilibrium equations. 

The continuum theory is applied here to describe linear elastic dielectric solids. A specific form 

of the free energy density is thus postulated as 

( ) ( )( ) ( ) ( ) , , 1/2 1/2 , ijkl i j

ik jl ij ρψ λ φ ξθ = ++ ^ u u PP , (22) 

where ^ and λ are linear-elastic moduli and inverse dielectric susceptibility, respectively, and 

symmetrized indices are in parentheses. Presently, we address only the response of the material 

in its paraelectric state, at temperatures above the Curie point. Thus, phase transformations, 

piezoelectricity, pyroelectricity, and spontaneous polarization are not considered. From equation 

22, bulk thermodynamic relations in equations 12, 13, and 3 reduce to 

( ) ,

ij ijkl

k l σ = ^ u , (23) 

j E P i ij = λ , (24) 

and 

0

i ij D ER j = ε ε , (25) 
5

with 1 1

0

ij ij ij

R ε δ ελ − − = + the relative permittivity (i.e., real dielectric constant). The vacancy and 

temperature-dependent contribution φ is assumed here to follow the universal relation for the 

chemical potential of an ideal mixture (Fried et al., 1977), most relevant for noninteracting 

species and small vacancy concentrations: 

φ θ θν N N G Nk v A T AB = + ( 0 ( ) ln ) . (26) 

In equation 26, NA is Avagadro’s number, B k is Boltzmann’s constant, 1/ NT = α is the number 

of atomic sites per unit volume, v is the mole fraction of vacancies, and G0 is the bulk Gibbs free 

energy of unstressed, defect-free substance. As implemented here, we assume simply that 

G c 0 = − θ ⎡ˆ ˆ (θ ηθ ) ( )⎤ ⎣ ⎦, (27) 

with specific heat cˆ and specific entropy ηˆ . Upon application of equation 26, the kinetic 

equation for bulk diffusion then follows directly from equation 16 as 

( , ) i ij Q d k ezE =− − Bj j θχ χ . (28) 

2.2 One-Dimensional Model 

The physical system analyzed is a substance of length (i.e., thickness) T with applied bias 

voltage V, as shown in figure 1. Alternatively to voltage boundary conditions, electric fields may 

be applied at the boundaries (not shown in figure 1). No mechanical tractions are applied, and 

isotropic material properties are assumed ( ij ij

R R ε = ε δ , ij ij d d = δ ), thus rendering the analysis one 

dimensional in thickness direction x. Correspondingly in this analysis, surface flux 0 qα = and 

surface tension τ = 0 . A constant temperature θ is assumed. 

Figure 1. One-dimensional problem domain. 

Under the preceding conditions, the governing electrostatic equations 1–3 reduce to 

( ) 2 2

0 0 / ˆ R R dD dx dE dx d dx = =− = ε ε εε φ ρ . (29) 

The vacancy evolution equations 28 and 15 become 
6

Q d k d dx ezE =− − ( Bθχ χ ( ) ) (30) 

and 

χ α  = − (dQ dx) , (31) 

with the domain size evolving, from equation 18, as 

a QT T Q  = −− − αχ χ ( ) () () /1 0/1 0 ( ) () () ( ) . (32) 

Spatial boundary conditions are applied as follows: 

( ) ( ) t t t Lt ˆ ˆ 0, , 0 = = , (33) 

( ) 0 φ 0,t V= or ( ) 0 E 0,t E = , (34) 

( ) , φ Tt V= T or ETt E ( , ) = T , (35) 

( ) 0 χ 0,t = χ or ( ) 0 QtQ 0, = , (36) 

and 

( ) , χ χ T t = T or QTt Q ( , ) = T . (37) 

The boundary fluxes Q0 and QT may be prescribed as constant values, or determined 

instantaneously from the model thermodynamics (equation 21): 

( ) ( ( )( )) 1 1 1 χ Q ez β ρψ ρφ χ ρ ψ χ α φ ˆ − − =− + + − ∂ ∂ + . (38) 

Initial conditions for the vacancy concentration and electrostatic charge density are applied as 

χ χ ( x,0) = ( x) (39) 

and 

ρ ρ ˆ ( x,0) = ( x), (40) 

where ρ = ezχ α/ should be imposed for consistency with equation 11. 

Note from equations 30 and 31 that the electric field may act as a local source or sink for 

vacancies, implying that ions may be impinged within or released from the lattice due to 

electrostatic forces. In some cases, it is convenient and realistic to impose an additional 

constraint that the total number of vacancies in the system remains constant, i.e., 

0 0

T

χdx T = χ ∫ , (41) 

where χ 0 is the initial average concentration. 
7

3. Numerical Methods 

The time duration of the problem is decomposed into a series of steps. For each time step, the 

electrostatic problem is solved, with the corresponding solution for the electric field used in the 

transient equations (equation 30) for updating the vacancy concentration. In the following time 

step, this updated vacancy concentration is used to compute the local electric charge entering 

Maxwell’s equations (equation 29). The analysis thus marches forward in time, with the 

numerical solutions of the equations of electrostatics and diffusion coupled in this manner. 

The 1-D spatial domain is discretized into n −1 increments of equal length ∆x , with n the total 

number of nodes, as shown in figure 2. 

Figure 2. Spatial discretization. 

3.1 Electrostatics 

A second-order accurate implicit scheme is used to solve Maxwell's static equations in each step. 

Within the domain 0 < <x T , equation 29 is written, to accuracy of order ( )2

∆x , using a finite 

difference approximation (Hoffman, 1992) for the second spatial derivative of φ : 

( )2

11 0 2 / ˆ i ii i R φ φ φ ρ εε x + − − + =− ∆ , (42) 

where subscripts denote node numbers. The solution of equation 42 for all nodes i is then 

written in matrix form as 

[ ] [ ] [ ] −1

φ = Α ρ , (43) 

where [φ ] is the n-dimensional solution vector of nodal electric potentials, [Α] is the n n ×

coefficient matrix that is generally tridiagonal and sparse, and [ρ ] is an n-dimensional vector 

that corresponds to the right hand side of equation 42. In the solution procedure, the dimensions 

of [φ ], [Α], and [ρ ] are reduced by one for each electrostatic potential boundary condition 

enforced a priori. Once the electric potential is known, the electric field is then computed 

numerically as 
8

1 1 2 ii i − ∆= − xE φ + φ − . (44) 

The following second-order finite difference approximations are used at the boundaries (i = 1 

and i = n): 

1 1 23 −2 34 ∆ =− + − xE φ φ φ , (45) 

and 

2 1 2 43 nn n n −∆ = − + xE φ − − φ φ . (46) 

3.2 Transient Diffusion 

A fully explicit scheme is used to integrate the transient vacancy concentration, with all 

necessary spatial derivatives obtained using a second-order accurate finite difference approach. 

The vacancy flux (equation 30) is written as follows, where subscripts denote node numbers: 

Q d d dx ezE i ii =− ∂ ∂ − ( ( ρψ χ )/ ) , (47) 

where for nodes in the domain 0 < <x T , 

ln iB i ρ ∂ψχ θ χ ∂ = k (48) 

and 

d dx k x ( ) ρψ χ θ χ χ ∂∂ = − ∆ i Bi i / ln ln / 2 ( + − 1 1 ) . (49) 

Note that in equation 47, the electric field Ei is obtained from the electrostatic solution given by 

equations 44–46. At the boundaries, the following approximations were used, analogous to 

equations 45 and 46: 

d dx k x ( ) ρψ χ θ χ χ χ ∂∂ = − + − ∆ 1 3 21 / ln 4ln 3ln / 2 B ( ) , (50) 

and 

d dx k x ( ) ρψ χ θ χ χ χ ∂∂ = − + ∆ n Bn n n / ln 4ln 3ln / 2 ( − − 2 1 ) . (51) 

The rate equation for concentration is approximated as follows for 0 < x < T : 

( 1 1 )/ 2 χ α i ii QQ x = −− ∆ + −  , (52) 

and at the boundaries as 

χ α 1 3 21 =− − + − ∆ ( QQQ x 4 3 /2 ) (53) 

and 

χ α  n n nn =− − + ∆ (Q QQ x − − 2 1 4 3 /2 ) . (54) 

Application of concentration boundary conditions simply entails 0 χi = , while application of 

constant flux conditions at a boundary eliminates the need to solve (equation 47) at that 
9

boundary. When transient flux conditions (equation 38) are applied, the following 

approximation suffices: 

( ) ( ( ) ) 1 1 1 ln ˆ χi i i ii i B i i Q k ez β ρψ ρ φ χ θ χ α φ − − =− + + − + , (55) 

where the free energy per unit volume is found from 

( ( )( )) 2

0 0 / 2 / ln / / ρψ ε ε θ χ α χ α i R i T AB i T i T A = ++ E N G Nk N N N . (56) 

Finally, the concentration is updated explicitly as 

tt t

i ii χ χχ t +∆ = + ∆  , (57) 

where ∆t is the time increment. 

The logarithmic form of the chemical potential, and flux equations 48–51, require that 0 χi > . 

This is achieved in practice by enforcing 

min 0 t t

i i χ i χ χ → ∨∈ < +∆  , (58) 

where t t χ i

+∆ is the projected updated concentration from equation 57 and 21

min χ 10− = is a nearnegligible, default minimum concentration. Optional global conservation condition (equation 

41) is imposed by rescaling equation 57: 

( ) ( )

1

0

1

1 n tt t t

i ii ii

i

n T tx t

n

χ χ χχ χχ

−

+∆

=

⎛ ⎞ − = +∆ ∆ +∆ ⎜ ⎟ ⎝ ⎠

  ∑ . (59) 

The domain size is updated from equation 32 as 

aQ Q  =− − − − αχ χ ( n n /1 /1 ( ) 1 1 ( )) (60) 

and 

( )/ 1 ( ) t t x T at n +∆ ∆ = +∆ −  . (61) 

For the solution of equation 42 in the next time step, the updated charge concentration is found 

from equation 11: 

ˆ / tt tt

i i ρ ezχ α +∆ +∆ = . (62) 

A fixed time increment is used throughout the analysis, i.e., 

max ∆tt m = / , (63) 

where max t and m are the total simulation time and number of time increments, respectively. 

Convergence and stability of the numerical solution dictate the practical choice of ∆t
10

(Chapra and Canale, 1998). Let / D dk = B R θ ξ , where R ξ is a reference concentration per unit 

volume, on the order of the minimum local concentration in the domain. Then one may write 

( )2

∆≤ ∆ t xD Λ / , (64) 

where Λ ≤1/2 guarantees convergence and stability. Setting Λ ≤1/4 ensures that the solution 

will not oscillate, while the choice Λ =1/6 has been shown, for the pure diffusion problem, to 

minimize truncation error (Carnahan et al., 1969). A description of what effect the unique 

electric field term i ezE in equation 47 may have on stability and convergence is not presently 

available. 

4. Software Manual 

In what follows, the source code structure, input files, and output files are documented. Note 

that for reference, the complete code is contained in the appendix. 

4.1 Code Structure 

The source code is written in the standard FORTRAN90 language. It has been compiled with 

Portland Group's FORTRAN compiler (i.e., pgf90) (http://www.pgroup.com/) and executed on a 

64-bit Linux workstation. Floating point values are represented in scientific double precision. 

This is a serial, as opposed to parallel, code. The program permits solution of the elliptic 

differential equations of electrostatics and/or the parabolic differential equations for transient 

diffusion. Standard SI units are used throughout. The software enables all of the features 

described in section 3, as well as a few additional options offering additional user flexibility. 

The code consists of the following routines: 

• the main routine, 'semiconductor_1D', which controls the grid spacing, memory allocation, 

and time incrementation. 

• subroutine 'elliptic', which solves Maxwell's equations. 

• subroutine 'parabolic', which solves the transient diffusion equations. 

• subroutine 'chem_potential', which computes the local free energy and its gradient. 

• subroutine 'LU_decomp', which decomposes the [Α] matrix in (43) to LU-form. 

• subroutine 'LU_backsub', which solves equation 43 using a back substitution technique. 

Figure 3 is a flowchart demonstrating execution of the code. 
11

Figure 3. Flowchart for code execution. 

4.2 Input 

Up to three input files are required. The file read by the main routine, 'parameters.inp', is always 

required. Its format is as follows: 

'parameters.inp'

Line 1: problem type. Specify 'maxwell' to solve only the electrostatic equations, 'fick' to solve 

only the diffusion equations, or 'mixed' to solve the fully coupled electrostatic-diffusion 

problem 

Line 2: number of nodes n 

Line 3: solution end time max t [s] 

Line 4: fixed time increment ∆t [s] 

Line 5: number of time steps in between each write to the output files, t n

Line 6: domain thickness T [m] 

Line 7: temperature θ [K] 

Line 8: relative dielectric permittivity R,0 ε

Line 9: factor κ enabling permittivity to depend on local electric field (Johnson, 1962): 

( ) 1/3 2

,0 1 εε κ R R E −

= +

Line 10: diffusivity d 

Line 11: diffusivity ˆ

d for electrostatic contribution (set ˆd = d for consistency with equation 

47): 

( ( ) ) ˆ / Q d d dx dezE i ii =− ∂ ∂ + ρψ χ

Line 12: factor for conversion of concentration from [ppm] to [m–3], equal to 6 10 /α − [m–3] 

Line 13: factor for conversion of concentration from [ppm] to volume fraction, nominally 6 10−

Line 14: valence contribution per vacancy, z

Line 15: specific entropy ηˆ [J/mol K] 

Line 16: specific heat capacity cˆ [J/mol K] 
12

The file 'electrostatic.inp' is read only if the problem type is 'maxwell' or 'mixed'. Its format is 

specified as follows: 

'electrostatic.inp'

Line 1: boundary condition at x = 0 . Use 'potential' to supply a voltage φ , or 'flux' to supply 

an electric field E 

Line 2: value of φ ( ) x = 0 [V] for a 'potential' condition, or value of E x( = 0) [V/m] for a 

'flux' condition 

Line 3: boundary condition at x = T . Use 'potential' to supply a voltage φ , or 'flux' to supply 

an electric field E 

Line 4: value of φ ( ) x = T [V] for 'potential' condition, or value of E ( x T = ) [V/m] for 'flux' 

condition 

Line 5: text heading (not used) 

Line 6: initial charge density at node 1, 0

1 ˆt ρi

=

= [C/m3

] 

Line 7: initial charge density at node 2, 0

2 ˆt ρi

=

= [C/m3

] 

.

.

.

Line n + 5 : initial charge density at node n, 0 ˆt ρi n

=

= [C/m3

] 

The file 'diffusion.inp' is read only if the problem type is 'fick' or 'mixed'. Its format is given as 

follows: 

'diffusion.inp'

Line 1: boundary condition at x = 0 . Use 'potential' to supply a constant vacancy 

concentration c , 'flux' to supply a constant vacancy flux Q, or 'special' to apply equation 55. 

Line 2: value of c x( ) = 0 [ppm] for a 'potential' condition, value of Q x( = 0)[ppm m/s] for a 

'flux' condition, or value of β in (55) for 'special' condition 

Line 3: boundary condition at x = T . Use 'potential' to supply a constant vacancy 

concentration c , 'flux' to supply a constant vacancy flux Q, or 'special' to apply equation 55. 

Line 4: value of cx T ( ) = [ppm] for a 'potential' condition, value of Qx T ( = ) [ppm m/s] for a 

'flux' condition, or value of β in (55) for 'special' condition 

Line 5: flag to apply constraint (59): 'yes' or 'no' 

Line 6: initial vacancy concentration at node 1, 0

1

t

i c =

= [ppm] 

Line 7: initial vacancy concentration at node 2, 0

2

t

i c =

= [ppm] 

.

.

.

Line n + 5 : initial vacancy concentration at node n, t 0

i n c =

= [ppm] 
13

4.3 Output 

Up to four output files are generated, each containing computational results. The file 'outelec.txt' 

is written for problems of type 'maxwell' or 'mixed', and provides the following: 

'outelec.txt'

At output request timest t t nt t n t n t nt t == = =− = 0, / , 2 / ,..., 1 / , max max max max t tt t ( ) : 

Column 1: node i 

Column 2: time t [s] 

Column 3: position x [m] 

Column 4: electric potential φi [V] 

Column 5: electric field Ei [V/m] 

Column 6: electric charge ˆρi [C/m3] 

The file 'outdiff.txt' is written for problems of type 'fick' or 'mixed', and provides the following: 

'outdiff.txt'

At output request timest t t nt t n t n t nt t == = =− = 0, / , 2 / ,..., 1 / , max max max max t tt t ( ) : 

Column 1: node i 

Column 2: time t [s] 

Column 3: position x [m] 

Column 4: vacancy concentration c [ppm] 

Column 5: concentration rate c [ppm/s] 

Column 6: vacancy flux Q [ppm m/s] 

The file 'earray.txt' contains the electric field data in tabular form for easy generation of 3-D 

mesh plots. Zeros are returned if the 'elliptic' routine is not called. The format is as follows: 

'earray.txt'

Rows correspond to node numbers i, while columns correspond to time increments. 

Column 1: x coordinate of each node i, at max t t =

Column 2: electric field Ei [kV/cm] at each node i, at t = 0

Column 3: electric field Ei [kV/cm] at each node i, at max / t tt n =

.

.

.

Column 2 t n + : electric field Ei [kV/cm] at each node i, at max t t =

The file 'carray.txt' contains the vacancy concentration data in tabular form for easy generation of 

3-D mesh plots. Zeros are returned if the 'parabolic' routine is not called. The format is listed as 

follows: 
14

'carray.txt'

Rows correspond to node numbers i, while columns correspond to time increments. 

Column 1: x coordinate of each node i, at max t t =

Column 2: concentration i c [ppm] at each node i, at t = 0

Column 3: concentration i c [ppm] at each node i, at max / t tt n =

.

.

.

Column 2 t n + : concentration i c [ppm] at each node i, at max t t =

5. Sample Problem-BST Film 

An example problem is discussed here to demonstrate the software capabilities and input file 

format. The physical system analyzed here is a thin film of uniform thickness T =100 nm

(figure 1). Isothermal conditions are assumed with θ = 298 K , a temperature at which undoped 

BST remains cubic in phase for molar concentrations of Sr greater than 0.3 (Tinte et al., 2004). 

Doping with Mg further lowers the Curie point (Cole et al., 2003). 

Requisite material properties are listed in table 1, deemed representative of particular 

composition Ba0.6Sr0.4TiO3. The dielectric constant R ε is chosen as representative of the BST 

film systems of interest (Cole et al., 2003), while thermodynamic properties cˆ and ηˆ of the bulk 

substance are used (Todd and Lorenson, 1952). The vacancy diffusivity d is computed from 

/ / D k ez dk = = μ B BR θ θξ , (65) 

where D is a thermally-activated diffusivity associated with Fick's first law, μ is the drift 

mobility of ionic charges, and R ξ is a reference defect concentration. We assume here, from 

drift mobility data on 100 nm-thick BST films (Zafar et al., 1998), a constant value of 

8 –1 –1 –1 2(10) V m s μξ R− = , leading to the value of d reported in table 1. Parameter α is estimated 

from consideration of the crystal structure and lattice parameters, while the value of z indicates 

that each oxygen vacancy contributes a charge of magnitude of two free electrons. 

Table 1. Properties of BST film at 298 K. 

Parameter Value Parameter Value 

R ε 500 d 10 6.24(10) J–1 m–1 s–1

cˆ 100 J mol–1 K–1 α 29 3 10 m −

ηˆ 115 J mol–1 K–1 z 2.0 
15

The particular initial-boundary value problem includes the following boundary conditions, 

corresponding to shorted or ground electrodes at the film boundaries, which themselves are 

constrained to remain impenetrable to vacancy flow: 

φ ( ) 0, 0 t = , φ (T t, 0 ) = , (66) 

 Q t ( ) 0, 0 = , QTt ( , 0 ) = . (67) 

Initial conditions are applied as 

( ) 6 χ x,0 10− = (68) 

and 

( ) 3 ρˆ x,0 32043 C/m = . (69) 

The spatially-constant value of ( ) 6 χ x,0 10− = corresponds to an initial concentration of 

0 c =1 ppm . Additionally, vacancy conservation constraint (59) is activated throughout the 

analysis. The spatial domain is decomposed into 300 segments of equal length ∆ = x 1/ 3 nm , 

with a total of 301 nodes. The time domain is chosen as 0 1 s ≤ t ≤ μ . Input files used for the 

computation are as follows: 

'parameters.inp' 

mixed 

301 

1.0d-6 

1.0d-12 

100000 

1.0d-7 

298 

500 

0.0 

6.24d10 

6.24d10 

1.0d23 

1.0d-6 

2.0 

115.0 

100.0 

'electrostatic.inp' 

potential 

0.0 

potential 

0.0 

initial_charge_distribution 

32043.0 

32043.0 

.

.

.

32043.0 
16

'diffusion.inp' 

flux 

0.0 

flux 

0.0 

yes 

1.0 

1.0 

.

.

.

1.0 

1.0 

The solution is shown in figure 4, with concentration and electric field mesh plots generated 

from output files 'carray.txt' and 'earray.txt', respectively. Although a steady-state concentration 

is not reached in the short duration of this transient simulation, a buildup of vacancies at the 

boundaries x = 0 and x T = =100 nm occurs, as the charges migrate from their initially uniform 

distribution towards the boundaries where the voltage 0 φ = . Vacancy migration has a negligible 

effect on the local electric field E, however, which is initially linear due to the constant initial 

charge distribution and exhibits an average value of zero. 

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

0

20

40

60

80

100

0.0 0.2 0.4 0.6 0.8

c [ppm]

x [nm]

t [μs]

0.8 

0.9 

1.0 

1.1 

1.2 

1.3 

1.4 

-4

-3

-2

-1

0

1

2

3

4

0

20

40

60

80

100

0.0 0.2 0.4 0.6 0.8 1.0

E [kV/cm]

x [nm]

t [μs]

-4 

-3 

-2 

-1 

0 

1 

2 

3 

4 

Figure 4. Vacancy concentration and electric field: 0 c = 1 ppm , V = 0 . 

Vacancy concentration distributions can be explained upon further examination of chemical 

potential μ of equation 17. The term ρ ∂ψ ξ∂ establishes an energetic penalty for large 

concentrations and forces the vacancies to diffuse towards a uniform state. The term ezφ

accounts for electrostatic interactions due to charges of the vacancies, and leads to their 

migration towards locations where the potential φ is small. The final distribution of vacancies is 

thus established from a balance of the contributions of these two terms in the chemical potential. 
17

6. Conclusions 

A continuum model for elastic dielectric semiconductors with mobile charged point vacancies 

has been presented. The complete 3-D theoretical framework includes the physics of continuum 

elasticity, electrostatics, mass diffusion, and charged defect kinetics. 

A restricted version of the theory has been implemented in a finite difference code allowing 

solution of the elliptic equations of electrostatics of dielectrics coupled with the transient 

parabolic equation of charged diffusion. The numerical analysis is at present limited to one 

spatial dimension. 

User documentation has been provided for the computer implementation. Given in the report are 

descriptions of the source code structure, user instructions, and representative input files for 

analysis of BST thin films containing charged oxygen vacancies. The source code has been 

included in the appendix. 
18

7. References 

Buchanan, D. A. IBM J. Res. Develop. 1999, 43, 245. 

Carnahan, B.; Luther, H. A.; Wilkes, J. O. Applied Numerical Methods; Wiley: New York, 

1969. 

Chapra, S. C.; Canale, R. P. Numerical Methods for Engineers with Programming and Software 

Applications; 3rd ed.; McGraw-Hill: Boston, 1998. 

Cole, M. W.; Nothwang, W. D.; Hubbard, C.; Ngo, E.; Ervin, M. J. Appl. Phys. 2003, 93, 9218. 

Damjanovic, D. Rep. Prog. Phys. 1998, 61, 1267. 

Devonshire, A. F. Phil. Mag. 1954, 3, 85. 

Fried, V.; Hameka, H. F.; Blukis, U. Physical Chemistry; Macmillan: New York, 1977. 

Grinfeld, M. A.; Hazzledine, P. M. Europhys. Lett. 1997, 37, 409. 

Hoffman, J. D. Numerical Methods for Engineers and Scientists; McGraw-Hill: Boston, 1992. 

Johnson, K. M. J. Appl. Phys. 1962, 33, 2826. 

Lifshitz, I. M. Sov. Phys. JETP, 1963, 17, 909. 

Stratton, J. A. Electromagnetic Theory; McGraw-Hill: New York, 1941. 

Tinte, S.; Stachiotti, M. G.; Phillpot, S. R.; Sepliarsky, M.; Wolf, D.; Migoni, R. L. J. Phys. 

Condens. Matter 2004, 16, 3495. 

Todd, S. S.; Lorenson, R. E. J. Amer. Chem. Soc. 1952, 74, 2043. 

Toupin, R. A. J. Rat. Mech. Anal. 1956, 5, 849. 

Zafar, S.; Jones, R. E.; Jiang, B.; White, B.; Chu, P.; Taylor, D.; Gillespie, S. Appl. Phys. Lett.

1998, 73, 175. 
19

Appendix. Source Code 

 

This appendix appears in its original form, without editorial change. 
20

INTENTIONALLY LEFT BLANK.
21

 program semiconductor_1D 

c 

 implicit none 

c 

 character(len=50) problem_type 

c 

 integer i,j,k,l,m,n 

 integer nmax !max # of nodes in grid 

 integer nsteps !max * of timesteps 

 integer nout,kount,count !# steps between output requests 

 integer ke,kc !output flags 

c 

 real*8 t,tmax,dt !time,max time, delta-time 

 real*8 p1,p2 !permittivity constants 

 real*8 d1,d2 !bulk diffusivity constants 

 real*8 z !valence per vacancy 

 real*8 temp !temperature (K) 

 real*8 alpha !convert to per m^3 from ppm 

 real*8 beta !convert to vol. fraction from ppm 

 real*8 L0,LC !initial and current domain length (m) 

 real*8 Ldot !velocity of domain (m/s) 

 real*8 dx !grid spacing (m) 

 real*8 x !coordinate (nm) 

 real*8 pfree !permittivity of free space (C/Vm) 

 real*8 kb !Boltzmann's constant (J/K) 

 real*8 e !charge of an electron (C) 

 real*8 s !specific entropy (J/molK) 

 real*8 cp !specific heat (J/molK) 

 real*8 cave0 !initial avg concentration (ppm) 

 real*8,pointer,dimension(:)::dconc !cumulative change in conc (ppm) 

 real*8,pointer,dimension(:)::dconst !relative dielectric constant 

 real*8,pointer,dimension(:)::efield !electric field (V/m) 

 real*8,pointer,dimension(:)::potential !electric potential (V) 

 real*8,pointer,dimension(:)::charge !electric charge (C) 

 real*8,pointer,dimension(:,:)::earray !electric field output matrix 

 real*8,pointer,dimension(:,:)::carray !concentration output matrix 

c 

 500 format(12e14.5) 

c 

 pfree = 8.854187817d-12 

 kb = 1.3806505d-23 

 e = 1.60217653d-19 

c 

 open(unit=1,file='parameters.inp',status='unknown') 

 read(1,*)problem_type 

 read(1,*)nmax 

 read(1,*)tmax 

 read(1,*)dt 

 read(1,*)nout 

 read(1,*)L0 

 read(1,*)temp 

 read(1,*)p1 

 read(1,*)p2 

 read(1,*)d1 

 read(1,*)d2 

 read(1,*)alpha 

 read(1,*)beta 
22

 read(1,*)z 

 read(1,*)s 

 read(1,*)cp 

 close(1) 

c 

 print*,'Problem ',problem_type 

 print*,'Domain initial length ',L0 

 print*,'Number of nodes ',nmax 

 print*,'Max time and delta_t ',tmax,dt 

 print*,'Bulk diffusion constants ',d1,d2 

 print*,'Bulk dielectric constants ',p1,p2 

 print*,'Valence per vacancy ',z 

 print*,'Temperature ',temp 

 print*,'Entropy ',s 

 print*,'Specific heat',cp 

 print*,'Conversion ppm to per cubic meter ',alpha 

 print*,'Conversion ppm to volume fraction ',beta 

c 

 nsteps= tmax/dt 

 nsteps= nsteps/nout+1 

 print*,'Number of output steps ',nsteps 

c 

 allocate(dconc(nmax)) 

 allocate(dconst(nmax)) 

 allocate(efield(nmax)) 

 allocate(potential(nmax)) 

 allocate(charge(nmax)) 

 allocate(earray(nmax,nsteps)) 

 allocate(carray(nmax,nsteps)) 

c 

 t=0.0 

 do i=1,nmax 

 dconc(i)=0.0 

 dconst(i)=p1 

 efield(i)=0.0 

 potential(i)=0.0 

 charge(i)=0.0 

 end do 

 do i=1,nmax 

 do j=1,nsteps 

 earray(i,j)=0.0 

 carray(i,j)=0.0 

 end do 

 end do 

 

 Ldot=0.0 

 LC=L0 

 kount=nout 

 count=nout 

 ke=0 

 kc=0 

c 

 100 if(t.le.tmax) then 

c 

 dx=LC/(nmax-1.) 

c 
23

 if(problem_type.eq.'maxwell')then 

 call elliptic(dconc,dx,temp,p1,p2,z,e,charge,dt,pfree, 

 1 alpha,t,dconst,nmax,efield,potential, 

 2 nout,count,tmax,earray,nsteps,ke) 

 end if 

c 

 if(problem_type.eq.'fick')then 

 call parabolic(dconc,dconst,d1,d2,temp,nmax,dx,Ldot,s,cp, 

 1 beta,t,dt,z,e,kb,pfree,alpha,efield,charge,potential,tmax, 

 2 nout,kount,cave0,carray,nsteps,kc) 

 end if 

c 

 if(problem_type.eq.'mixed')then 

 call elliptic(dconc,dx,temp,p1,p2,z,e,charge,dt,pfree, 

 1 alpha,t,dconst,nmax,efield,potential, 

 2 nout,count,tmax,earray,nsteps,ke) 

c 

 call parabolic(dconc,dconst,d1,d2,temp,nmax,dx,Ldot,s,cp, 

 1 beta,t,dt,z,e,kb,pfree,alpha,efield,charge,potential,tmax, 

 2 nout,kount,cave0,carray,nsteps,kc) 

 end if 

c 

 t = t + dt 

 LC = LC+(Ldot*dt) 

 go to 100 

 end if 

c 

 open(unit=5,file='earray.txt',status='unknown') 

 open(unit=6,file='carray.txt',status='unknown') 

 x=0.0 

 do i=1,nmax 

 write(5,500)x,(earray(i,j),j=1,nsteps) 

 write(6,500)x,(carray(i,j),j=1,nsteps) 

 x=x+dx*1.0d9 

 end do 

 close(5) 

 close(6) 

c 

 deallocate(dconc) 

 deallocate(dconst) 

 deallocate(efield) 

 deallocate(potential) 

 deallocate(charge) 

 deallocate(earray) 

 deallocate(carray) 

c 

 end program semiconductor_1D 

c 

c################################################################### 

c 

 subroutine elliptic(dconc,dx,temp,p1,p2,z,e,charge,dt, 

 1 pfree,alpha,t,dconst,nmax,efield,potential, 

 2 nout,count,tmax,earray,nsteps,ke) 

c 

 implicit none 

c 

 character(len=50)bcleft,bcright,dummychar 
24

 integer i,j,k,l,m,n 

 integer nmax,adim,nout,count,nsteps,ke 

 integer,pointer,dimension(:)::index 

 real*8 dx,temp,p1,p2,z,e,pfree,alpha,t,dt,tmax,tout 

 real*8 valleft,valright 

 real*8 dconc(nmax) 

 real*8 dconst(nmax) 

 real*8 charge0(nmax) 

 real*8 charge(nmax) 

 real*8 potential(nmax) 

 real*8 efield(nmax) 

 real*8 x(nmax) 

 real*8 earray(nmax,nsteps) 

 real*8 factor,energy,d 

 real*8,pointer,dimension(:,:)::amatrix 

 real*8,pointer,dimension(:)::bvector 

 real*8,pointer,dimension(:)::xvector 

 200 format(i6,5e14.5) 

c 

 open(unit=2,file='electrostatic.inp',status='unknown') 

 read(2,*)bcleft 

 read(2,*)valleft 

 read(2,*)bcright 

 read(2,*)valright 

 read(2,*)dummychar 

 do i=1,nmax 

 read(2,*)charge0(i) 

 end do 

 close(2) 

c 

 adim=nmax 

c 

 if(bcleft.eq.'potential')then 

 adim=adim-1 

 end if 

c 

 if(bcright.eq.'potential')then 

 adim=adim-1 

 end if 

c 

 allocate(amatrix(adim,adim)) 

 allocate(bvector(adim)) 

 allocate(xvector(adim)) 

 allocate(index(adim)) 

c 

 do i=1,nmax 

 charge(i)=charge0(i)+e*z*dconc(i)*alpha 

 x(i)=dfloat(i-1)*dx 

 end do 

c 

 do i=1,adim 

 do j=1,adim 

 amatrix(i,j)=0.0 

 end do 

 amatrix(i,i)=-2.0 

 k=i+1 

 m=i-1 
25

 if(i.lt.adim)amatrix(i,k)=1.0 

 if(i.gt.1)amatrix(i,m)=1.0 

 end do 

c 

 factor=-dx*dx/pfree 

 do i=1,adim 

 j=i 

 if(bcleft.eq.'potential')j=j+1 

 bvector(i)=charge(j)*factor/dconst(j) 

 end do 

c 

 if(bcleft.eq.'potential') bvector(1)=bvector(1)-valleft 

c 

 if(bcright.eq.'potential')then 

 bvector(adim)=bvector(adim)-valright 

 end if 

c 

 if(bcleft.eq.'flux') then 

 bvector(1)=bvector(1)-2.*valleft*dx 

 amatrix(1,2)=2.0 

 end if 

c 

 if(bcright.eq.'flux') then 

 bvector(adim)=bvector(adim)+2.*valright*dx 

 i=adim-1 

 amatrix(adim,i)=2.0 

 end if 

c 

 call LU_Decomp(adim,amatrix,index) 

 call LU_BackSub(adim,amatrix,index,bvector) 

c 

 do i=1,adim 

 j=i 

 if(bcleft.eq.'potential')then 

 potential(1)=valleft 

 j=i+1 

 end if 

 potential(j)=bvector(i) 

 if(bcright.eq.'potential') potential(nmax)=valright 

 end do 

c 

 k=nmax-1 

 do i=2,k 

 j=i-1 

 m=i+1 

 efield(i)=(potential(j)-potential(m))/(2.*dx) 

 end do 

 efield(1)=(potential(1)-potential(2))/dx 

 efield(nmax)=(potential(k)-potential(nmax))/dx 

 if(bcleft.eq.'flux')efield(1)=valleft 

 if(bcright.eq.'flux')efield(nmax)=valright 

c 

 energy=0.0 

 do i=1,k 

 energy=energy+dconst(i)*efield(i)*efield(i)*dx 

 end do 

c 
26

 d=energy*dx*dfloat(nmax-1)/(potential(1)-potential(nmax))**2. 

c 

 do i=1,nmax 

 dconst(i)=p1/((1.+p2*efield(i)*efield(i))**(0.3333333333)) 

 end do 

c 

 tout=tmax-dt 

c 

 if(count.eq.nout.or.t.gt.tout)then 

 ke=ke+1 

 count=0 

 if(t.lt.dt)then 

 open(unit=3,file='outelec.txt', 

 1 status='unknown') 

 else 

 open(unit=3,file='outelec.txt',position='append', 

 1 status='unknown') 

 end if 

 do i=1,nmax 

 write(3,200)i,t,x(i),potential(i),efield(i),charge(i) 

 earray(i,ke)=efield(i)/1.0d5 

 end do 

 close(3) 

 print*,'Effective dielectric constant ',d 

 end if 

c 

 count=count+1 

c 

 deallocate(amatrix) 

 deallocate(bvector) 

 deallocate(xvector) 

 deallocate(index) 

c 

 return 

 end 

c 

c-################################################################### 

c 

 subroutine parabolic(dconc,dconst,d1,d2,temp,nmax,dx,Ldot,s,cp, 

 1 beta,t,dt,z,e,kb,pfree,alpha,efield,charge,potential,tmax, 

 2 nout,kount,cave0,carray,nsteps,kc) 

c 

 implicit none 

c 

 character(len=50)bcleft,bcright,constraint 

 integer i,j,k,l,m,n 

 integer nmax,nout,kount,nsteps,kc 

 real*8 dconc(nmax) 

 real*8 dconst(nmax) 

 real*8 efield(nmax),charge(nmax),potential(nmax) 

 real*8 d1,d2,temp,dx,t,dt,z,e,kb,pfree,beta 

 real*8 alpha,Ldot,s,cp,tmax,tout 

 real*8 valleft,valright 

 real*8 conc0(nmax) 

 real*8 conc(nmax) 

 real*8 cdot(nmax),gradc(nmax) 

 real*8 diff1(nmax),diff2(nmax),q(nmax),x(nmax) 
27

 real*8 denom,phi,dphi 

 real*8 c,cdotmax,cdottemp,cave,cave0,L0,LC 

 real*8 test,scale 

 real*8 carray(nmax,nsteps) 

 300 format(i6,5e14.5) 

c 

 open(unit=4,file='diffusion.inp',status='unknown') 

 read(4,*)bcleft 

 read(4,*)valleft 

 read(4,*)bcright 

 read(4,*)valright 

 read(4,*)constraint 

 do i=1,nmax 

 read(4,*)conc0(i) 

 end do 

 close(4) 

c 

 if(t.eq.0.0)then 

 cave0=0.0 

 L0=0.0 

 do i=1,nmax 

 cave0=cave0+dx*conc0(i) 

 L0=L0+dx 

 end do 

 cave0=cave0/L0 

 end if 

c 

 test = 1./6.*dx*dx/(d1*kb*temp) 

c if(dt.gt.test)print*,'Warning: dt may be unstable!',test 

c 

 do i=1,nmax 

 x(i)=dfloat(i-1)*dx 

 conc(i)=conc0(i)+dconc(i) 

 c=conc(i) 

 call chem_potential(c,temp,kb,alpha,s,cp,phi,dphi) 

 diff1(i)=-d1*dphi 

 diff2(i)=d2*e*z*efield(i) 

 end do 

c 

 k=nmax-1 

 do i=2,k 

 m=i+1 

 n=i-1 

 gradc(i)=(diff1(m)-diff1(n))/(2.*dx) 

 q(i)= gradc(i)+diff2(i) 

 end do 

c 

 if(bcleft.eq.'potential')then 

 gradc(1)=(-3.*diff1(1)+4.*diff1(2)-diff1(3))/(2.*dx) 

 q(1)= gradc(1)+diff2(1) 

 end if 

c 

 if(bcright.eq.'potential')then 

 i=nmax-1 

 j=nmax-2 

 gradc(nmax)=(3.*diff1(nmax)-4.*diff1(i)+diff1(j))/(2.*dx) 

 q(nmax)= gradc(nmax)+diff2(nmax) 
28

 end if 

c 

 if(bcleft.eq.'flux')q(1)=valleft 

 if(bcright.eq.'flux')q(nmax)=valright 

c 

 if(bcleft.eq.'special')then 

 c=conc(1) 

 call chem_potential(c,temp,kb,alpha,s,cp,phi,dphi) 

 q(1)=valleft/(1.-conc(1)*beta)* 

 1 (phi+potential(1)*charge(1)+0.5*efield(1)*efield(1) 

 2 *dconst(1)*pfree+(1./beta-conc(1))*dphi/beta+e*z*potential(1)* 

 3 (1.-conc(nmax)*beta)) 

 end if 

c 

 if(bcright.eq.'special') then 

 c=conc(nmax) 

 call chem_potential(c,temp,kb,alpha,s,cp,phi,dphi) 

 q(nmax)=valright/(1.-conc(nmax)*beta)*(phi+potential(nmax)* 

 1 charge(nmax)+0.5*efield(nmax)*efield(nmax)*dconst(nmax)*pfree 

 2 +(1./beta-conc(nmax))*dphi/beta+e*z*potential(nmax)* 

 3 (1.-conc(nmax)*beta)) 

 end if 

c 

 k=nmax-1 

 do i=2,k 

 m=i+1 

 n=i-1 

 cdot(i)=-(q(m)-q(n))/(2.*dx) 

 end do 

c 

 cdot(1)=-(-3.*q(1)+4.*q(2)-q(3))/(2.*dx) 

 i=nmax-1 

 j=nmax-2 

 cdot(nmax)=-(3.*q(nmax)-4.*q(i)+q(j))/(2.*dx) 

c 

 if(bcleft.eq.'potential')cdot(1)=0.0 

 if(bcright.eq.'potential')cdot(nmax)=0.0 

c 

 do i=1,nmax 

 cdottemp=cdot(i) 

 conc(i)=conc(i)+cdot(i)*dt 

 if (conc(i).lt.1.d-15)cdot(i)=0.0 

 conc(i)=conc(i)-cdottemp*dt 

 end do 

c 

 if (constraint.eq.'yes')then 

 cave=0.0 

 LC=0.0 

 do i=1,nmax 

 dconc(i)=dconc(i)+cdot(i)*dt 

 cave=cave+(conc0(i)+dconc(i))*dx 

 LC=LC+dx 

 end do 

 cave=cave/LC 

 scale=cave0/cave 

 do i=1,nmax 

 dconc(i)=dconc(i)*scale-conc0(i)*(1.-scale) 
29

 end do 

 else 

 do i=1,nmax 

 dconc(i)=dconc(i)+cdot(i)*dt 

 end do 

 end if 

c 

 cdotmax = 0.0 

 cave = 0.0 

 n=0 

 do i=1,nmax 

 if(abs(cdot(i)).gt.cdotmax)then 

 cdotmax=abs(cdot(i)) 

 n=i 

 end if 

 cave = cave + conc(i) 

 end do 

 cave = cave/dfloat(nmax) 

c 

 if(bcright.eq.'special'.or.bcleft.eq.'special')then 

 Ldot=q(1)/(1.-conc(1)*beta)+q(nmax)/(1.-conc(nmax)*beta) 

 else 

 Ldot=0.0 

 endif 

c 

 tout = tmax-dt 

c 

 if(kount.eq.nout.or.t.gt.tout)then 

 kount=0 

 kc=kc+1 

 if(t.lt.dt) then 

 open(unit=4,file='outdiff.txt', 

 1 status='unknown') 

 else 

 open(unit=4,file='outdiff.txt',position='append', 

 1 status='unknown') 

 end if 

 print*,'Time',t,'C average',cave 

 print*,'Max conc. rate', cdotmax, n 

 if(constraint.eq.'yes')print*,'Constrained rate:',test 

 if(bcright.eq.'special'.or.bcleft.eq.'special')then 

 print*,'Ldot',Ldot 

 end if 

 do i=1,nmax 

 write(4,300)i,t,x(i),conc(i),cdot(i),q(i) 

 carray(i,kc)=conc(i) 

 end do 

 close(4) 

 end if 

c 

 kount=kount+1 

c 

 return 

 end 

c 

c-#################################################################### 

c 
30

 subroutine chem_potential(c,temp,kb,alpha,s,cp,phi,dphi) 

c 

 implicit none 

c 

 real*8 c,temp,kb,alpha,cp,s,phi,dphi 

 real*8 gtotal 

 real*8 xv 

 real*8 avagadro,r 

 real*8 nt 

c 

 avagadro = 6.02214199447d23 

 r = 8.314472 

c 

 nt = alpha*1.0d6 

 xv = c*alpha/nt 

c 

 gtotal = temp*(cp-s+r*log(xv)) 

c 

 phi = nt/avagadro*gtotal 

c dphi = kb*temp*c 

 dphi = kb*temp*log(c) 

c 

 return 

 end 

c 

c-#################################################################### 

c 

 subroutine LU_Decomp(n,a,index) 

c 

 implicit double precision (a-h,o-z) 

 dimension a(n,n), index(n), v(n) 

c 

 tiny = 1.0e-20 

c 

 do i = 1,n 

 a_max = 0.0 

 do j = 1,n 

 a_max = max(a_max,abs(a(i,j))) 

 end do !j 

 v(i) = 1.0 / a_max 

 end do !i 

c 

 do j = 1,n 

c 

 do i = 1,j-1 

 sum = a(i,j) 

 do k = 1,i-1 

 sum = sum - a(i,k) * a(k,j) 

 end do 

 a(i,j) = sum 

 end do 

c 

 a_max = 0.0 

 do i = j,n 

 sum = a(i,j) 

 do k = 1,j-1 

 sum = sum - a(i,k) * a(k,j) 
31

 end do 

 a(i,j) = sum 

 dummy = v(i) * abs(sum) 

 if ( dummy .gt. a_max ) then 

 imax = i 

 a_max = dummy 

 end if 

 end do 

c 

 if ( j .ne. imax ) then 

 do k = 1,n 

 dummy = a(imax,k) 

 a(imax,k) = a(j,k) 

 a(j,k) = dummy 

 end do 

 v(imax) = v(j) 

 end if 

 index(j) = imax 

c 

 if ( a(j,j) .eq. 0.0 ) a(j,j) = tiny 

 if ( j .ne. n ) then 

 dummy = 1.0 / a(j,j) 

 do i = j+1,n 

 a(i,j) = a(i,j) * dummy 

 end do 

 end if 

 end do !j 

c 

 return 

 end 

c 

c-################################################################## 

c 

 subroutine LU_BackSub(n,a,index,b) 

c 

 implicit double precision (a-h,o-z) 

 dimension a(n,n), index(n), b(n) 

c 

 ii = 0 

c 

 do i = 1,n 

 m = index(i) 

 sum = b(m) 

 b(m) = b(i) 

 if ( ii .ne. 0 ) then 

 do j = ii,i-1 

 sum = sum - a(i,j) * b(j) 

 end do 

 else if ( sum .ne. 0.0 ) then 

 ii = i 

 end if 

 b(i) = sum 

 end do 

c 

 do i = n,1,-1 

 sum = b(i) 

 if ( i .lt. n ) then 
32

 do j = i+1,n 

 sum = sum - a(i,j) * b(j) 

 end do 

 end if 

 b(i) = sum / a(i,i) 

 end do 

c 

 return 

 end 

c 

c==================================================================== 
NO. OF 

COPIES ORGANIZATION 

33

 1 DEFENSE TECHNICAL 

 (PDF INFORMATION CTR 

 ONLY) DTIC OCA 

 8725 JOHN J KINGMAN RD 

 STE 0944 

 FORT BELVOIR VA 22060-6218 

 1 US ARMY RSRCH DEV & 

 ENGRG CMD 

 SYSTEMS OF SYSTEMS 

 INTEGRATION 

 AMSRD SS T 

 6000 6TH ST STE 100 

 FORT BELVOIR VA 22060-5608 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 IMNE ALC IMS 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRD ARL CI OK TL 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

ABERDEEN PROVING GROUND

 1 DIR USARL 

 AMSRD ARL CI OK TP (BLDG 4600) 
NO. OF 

COPIES ORGANIZATION 

34

 38 DIR USARL 

 AMSRD ARL CI HC 

 P CHUNG 

 R NAMBURU 

 D GROVE 

 AMSRD ARL WM 

 J MCCAULEY 

 T WRIGHT 

 AMSRD ARL WM TA 

 S SCHOENFELD 

 AMSRD ARL WM TC 

 M FERMAN COKER 

 R COATES 

 AMSRD ARL WM TD 

 S BILYK 

 T BJERKE 

 D CASEM 

 J CLAYTON (5 CPS) 

 T CLINE 

 D DANDEKAR 

 W EDMANSON 

 M GREENFIELD 

 C GUNNARSSON 

 Y HUANG 

 K IYER 

 R KRAFT 

 B LOVE 

 S MCNEILL 

 H MEYER 

 R MUDD 

 M RAFTENBERG 

 E RAPACKI 

 M SCHEIDLER 

 S SEGLETES 

 T WEERASOORIYA 

 AMSRD ARL WM MA 

 M COLE 

 J ANDZELM 

 M VANLANDINGHAM 

 W NOTHWANG 

 AMSRD ARL WM MD 

 G GAZONAS