Documentation

This file provides general information about the workings of the code, for those who want to make modifications. No guarantee that any of this is correct; it was compiled while reading through the program to understand what it is doing, not to check the program. Also, some of it was based on an abbreviated 1985 version of the program and may not apply to the full 1986 version.

General notes

The code uses Angstrom, picoseconds, atomic mass, and electron volt units. So it is frequently converting between eV and AMU A²/ps². See ECONV in the variable name table.

Real*8 numerical precision is needed, or bcuint has notable round-off error. The same is almost certainly true for the other spline interpolation subroutines.

For the ideas behind the code, see Phys. Rev. B, Vol 42, pp 9458-9471 for formulae and their symbols. Also note, however

The values found for AD, AXL, ... for V_A in param.f are grossly incompatible with those in Brenner's paper.
Moreover, it appears that the potential V_A for CC in the program is made up of *three* of the exponential expressions as listed in the paper, (AD, AXL/ BD, BXL / CD, CXL), with expression A dominating for small r, being overtaken by C at 1.8, and B similar to C.
The repulsive potential seems to have an additional factor (1 + const/r), presumably to strongly prohibit atoms penetrating to the center of another atom. The const is small, about 0.3.
N_i^(C), N_i^(H), N_i^(t) are really N_ij values since j is excluded in the sums.
G_i for carbon has been made dependent on N_i^(t): the program switches smoothly to a different representation for G_i between 4.2 and 4.7 neighbors to the i-atom (including j), i.e. it switches in the range 3.2 < N_i^(t) < 3.7
δ_ij has been hardcoded in pibond to be 0.5 for CC, CH, HC, HH
The program seems to square the sums in N_ij^conj
The program includes a "torsional interaction" for CC bonds. The added torsional potential is:
V^tor = V_A(|r_ij|) A_ij(N_i^(t),N_j^(t),N_ij^conj) Σ_kl sin²(φ_ijkl) f_ik f_jl
where φ_ijkl is the angle between the ijk and the ijl planes, and the sum extends over k≠ij, l≠ij, |sin(θ_ijk)| > 0.1, and |sin(θ_jil)| > 0.1. The f-values are hardcoded to special values for H.
The conditions on the sines in V^tor imply that energy is NOT conserved when three atoms become almost along the same line. So, if there is a problem with energy conservation, try turning of the torsional potential by setting ndihed to 10 in Subroutines/Bond_Order/param.f.
Alternatively, the subroutine has been modified to allow the restrictions to be smoothed. Setting SINHI to 0.2 and SINLO to 0.0 smooths the cut-off at 0.1 over that range and energy should again be conserved:
    V^tor = V_A(|r_ij|) A_ij(N_i^(t),N_j^(t),N_ij^conj) Σ_kl sin²(φ_ijkl) f_ik f_jl Z_ijk Z_jil
where in the smoothing range
    Z_ijk = 1 + cos(π (sin(θ_ijk)-SINHI)/(SINHI-SINLO))
    Z_jil = 1 + cos(π (sin(θ_jil)-SINHI)/(SINHI-SINLO))
The program includes a Lennard-Jones potential: V_LJ= 4 ε_ij [(σ_ij/r_ij)¹²- (σ_ij/r_ij)⁶]
where ε is read in in K (corresponding to energy k K).
Interactions between two different type atoms use a geometrically average of the ε values of each atom and a straigth average σ.

The Lennard-Jones potential is set to zero for r_ij > 2.5 σ_ij and in the range r_ij < R_ij⁽²⁾ in which the bond order potential is nonzero. In the range R_ij⁽²⁾ < r_ij < 0.95 σ_ij;, the potential is furthermore replaced by a cubic which vanishes quadratically at R_ij⁽²⁾ and which meets the Lennard-Jones expression with continuous first derivative at 0.95 σ_ij.

The discontinuity at 2.5 σ_ij means that energy is not preserved. I put in a variable ISMTH which can be changed to a nonzero value to substract the potential value at the cut-off, thus making the potential continuous. Note however that this changes the forces in the range described by the cubic.

Common blocks

The program uses common blocks for most variables which are loaded into each subroutine through common_files.inc:

common_ch.inc: hydrocarbon potential common blocks
common_lj_new.inc: Lennard-Jones common blocks
common_tb.inc: tight binding
common_md.inc: general
common_leon.inc: added by me

I/O units

The following Fortran I/O units are in use by the code:

1: xmol.d (to be post converted to xmol.xyz)
9: output.d
11: coord.d
13: input.d
14-18: REBO potential splines
23: for bombard.f
41: (used in the periodic box straining modification)
50: overwrite.d
51: load.d
53: fort.53 debug file
55: max_ke.d
85: pair_energy.d eigenvectors.d
86: eigen_energies.d
87: dos.d
88: ldos.d

I moved pair_energy.d to I/O unit 84.

Files are opened in open.inc, closed in close.inc.

Subroutines

Here is a list of subroutines and what they do:

Bcuint: Bicubic interpolation of CLM; also returns derivatives.
Bere: Berendsen thermostat: KFLAG=1 damping; KFLAG=other(=5): rescale velocity.
Bombard: Adds a single atom from file bombard.d.
Caguts: Finds forces on atoms from pair potentials.
Ccorr: Call the corrector.
Cpred: Calls Pred (if there are nonrigid molecules) and increments time.
Damage: Checks whether the number of carbons less than 1.8 apart changes, indicating damage..
Gleq: Langevin (friction and random force).
Hoov: Evans-Hoover thermostat.
Ljcont: Seems to give net Lennard-Jones potential exerted by e semi-infinite continuum with a normal in the 'ndir' direction, at a location 'surf'.
Model: Calls Caguts (if not tight binding).
Mtable: Creates tables of the potentials V_A and V_R from the Brenner paper..
Param: C and H potentials.
Pibond: Finds bonding forces involving C and H atoms.
Pred: Nordsieck scheme prediction of z_n+1 from z_n only.
Radic: Tricubic interpolation of CLMN.
Rannum: Random number generator. Its argument is unused..
Read_data: Reads input.d and coord.d.
Set_md: Initialize time and such, also sets ENPR.
Setgle: Set langevin parameters..
Setin: Initialize kt, xmass, noa, sig, eps, pi, bolz, avo, epsi, econv.
Setpc: Set Nordsieck parameters (except F22=1).
Setpp: Calls param and mtable.
Setran: Initializes RNG, called once at the start..
Setrn: Random number initialization. Only called once, inside setran.
Thermos: Chooses thermostat routine.
Tor: Torsional interaction TLMN tricubic interpolation.
Write_data1: Writes to output.d.
Write_data2: Writes to unit 9, output.d, the step number, average energy per atom, some kinetic energy; writes to unit 6 (the screen) total potential energy. Modified by me to write more output to the screen less frequently.
Write_data3: Overwrites 11, coord.d, .
Xmol: Writes to 1, xmol.d data to be converted later to .xyz files..
Zero: Zero velocities.
vscale: Seems to be some attempt to select the scalings that give minimum potential energy.

Variable Names (General)

The meaning of variable names, as well as I can make out follows. See Phys. Rev. B, Vol 42, pp 9458-9471 for formulae and their symbols. Also, I define the additional variables:

e_ijk = exp(α_ijk [(r_ij - R_ij^e) - (r_ik-R_ik^e)])
f_ij = f_ij(|r_ij|)
V_ijkl^tor = V_A(|r_ij|) A_ij(N_i^(t),N_j^(t),N_ij^conj) sin²(φ_ijkl) f_ik f_jl (a term in the torsional potential)
b_ij is defined by the relation B_ij = b_ij^-δ

AD(KT1,KT2): D_ij^e S_ij exp(AXL R_ij^e) /(1 - S_ij) value 1
ATABLE(KT1,KT2,1 + r_ij/DDtab): attractive potential V_A/2 function table
ATT: = 3.2; the value of N_i^(t) at which the evaluation of G_i for 0<θ<109 starts to smoothly switch to a different polynomial
AXL(KT1,KT2): sqrt(2/S_ij) β_ij value 1
BD(KT1,KT2): D_ij^e S_ij exp(AXL R_ij^e) /(1 - S_ij) value 2
BXL(KT1,KT2): sqrt(2/S_ij) β_ij value 2
CD(KT1,KT2): D_ij^e S_ij exp(AXL R_ij^e) /(1 - S_ij) value 3
CLM(KT2,1+N_i^(H),1+N_i^(C),term): bicubic splines H_CC and H_CH in inter2d_iv.d
CLMN(1,1+N_i^(t),1+N_j^(t),N_ij^conj,term): tricubic splines 2 F_CC from inter3d_iv.d
CLMN(2,1+N_i^(t),1+N_j^(t),N_ij^conj,term): tricubic splines 2 F_CH from inter3d_ch.d
CLMN(3,1+N_i^(t),1+N_j^(t),N_ij^conj,term): tricubic splines 2 F_HH from inter3d_h.d
COM: (in Corr) velocity of center of gravity
COR(K,.): R0(I,.) - R0(j,.) by interacting atom pair index K
CUBE(3): cube dimensions, set to 1e20 if not periodic
CUBE2(3): cube dimensions/2, set to .5e20 if not periodic (never used)
CXL(KT1,KT2): sqrt(2/S_ij) β_ij value 3
DATABLE(KT1,KT2,1 + r_ij/DDtab): attractive potential's derivative * 2 table
DD(KT1,KT2): D_ij^e exp(DXL R_ij^e) /(1 - S_ij)
DDTAB(KT1,KT2): r_ij-spacing of the potential tables
DELLJ: table spacing for Lennard-Jones potential, 0.001 A
DELTA: time step
DELTSQ: DELTA²/2
DEXX1(K): derivative of the attractive potential V_A / 2 of a pair K
DRTABLE(KT1,KT2,1 + r_ij/DDtab): repulsive potential's derivative table
DWW(K): derivative of f_ij(r_ij) for a pair of neighbors
DXL(KT1,KT2): sqrt(2 S_ij) β_ij
EATOM(I): potential energy of an atom
ECONV: Conversion factor: 103.6434 eV per (amu A²/fs²)
EN(NSMAX): never used, intended to store an energy?
EPS(KT1,KT2): Initially depth ε of Lennard-Jones potential in K (i.e. energy kK), converted internally to 4 eV units
EPSS: strength constant in the LJ potential induced by a semi-infinite solid
ETOT(NSMAX): never used, probably intended to store the time history of the total energy
EXX1(K): attractive potential V_A / 2 of a pair of neighbors K
I2D: potential version number in inter2d_iv.d
I3D: potential version number in inter3d_iv.d, inter3d_h.d, inter3d_ch.d
IDUM: set to 0 by make_tube, set to 3 by xmol, maybe to indicate that only one parameter, position, exists in the .d file.
IGC: Interpolation index, 4 for 1>=c>-1/3, 3 for -1/3>=c>-1/2, 2 for -1/2>=c>-2/3, 1 for -2/3>=c, with c=cos(θ), used for the G_i function of carbon
IGH: Interpolation index, 3 for 1>=c>-1/2, 2 for -1/2>=c>-5/6, 1 for -5/6>=c>-2/3, with c=cos(θ), used for the G_i function of hydrogen
IN2(term,variable): powers of the terms of the 2D spline
IN2(term,variable): powers of the terms of the 3D spline
IPIB: number of times pibond has been called
IPOT: 1=REBO, 2=Tight binding
IT: in Caguts, the interpolation interval
ITD: potential version number in inter3dtors.d
ITR: 2 is nonmoving, 1=moving and thermostated, other=moving only
IVCT2B: first atom of the pair in the list of interacting atom pairs
JVCT2B: second atom of the pair in the list of interacting atom pairs (≠I)
KEND: number of interacting atom pairs
KFLAG: thermostat flag: -1 Langevin; 1 Berendsen; 2 zero velocity; 3 Evans-Hoover; 5 Berendsen after corrector; 6 minimize energy; 8 rescale to minimize energy
KT2: inverse array to KT
KT: atom code indexed by atom number (C=1, H=2, ...)
KTYPE: atom code kt, but by atom index
KVC: number of steps to take
LCHECK: 1=both C|H; 2=both Si|Ge; 0=one C|H, one Si|Ge, or distance > RMAX
LCHK: 1=update neighbors list; 2=still OK
LIST: Same as JVCT2B
LSTEP: step number
MAXKB: steps between data writing using write_data2
MLIST: atoms that get advanced in time (all atoms with itr≠2)
NABORS: start in (IVCT2B,JVCT2B,LIST) of the list of neigbors of an atom
NATX: number of atoms used in tight-binding matrices. MUST be set equal to np before compiling.
NDIR: direction of the surface normal of a neighboring semi-infinite continuum
NLA: never used or assigned a value. Probably number of Langevin atoms, to be ignored in adaptive time step modifications
NLIST: atoms that get damping and random forces added in Gleq, damping in Bere
NLMAX: maximum number elements in neighbor list
NMA: number of atoms in MLIST
NOA: number of atoms of each type
NP: number of atoms
NPMAX: maximum number of atoms
NRA: never used or assigned a value. Probably number of rigid atoms, to be ignored in adaptive time step modifications
NSMAX: (or NSMAS) maximum number of steps. Derived arrays never used.
NTA: number of atoms in NLIST
NTAB: array size for potential table look-up
NTYPES: maximum number of different types of atoms allowed for Len-Jon potential
NXMOL: number of steps between xmol writes
PSEED: random number seed
R0: position / Nordsieck component 0
R0L: position when the neighborhood list was updated
R1: velocity when read from 11, internally converted into Nordsieck component 1
R2: Nordsieck component 2
R3: Nordsieck component 3
R4: Read from coord.d and written to it but never used anywhere, all zeros
RB1: distance R_ij⁽¹⁾ at the start of the bond order potential decay
RB2: distance R_ij⁽²⁾ at which the bond order potential is gone
RCOR: distance between a pair of interacting atoms
REG: exp(α_ijk [R_ik^e - R_ij^e])
RLIST: square interaction distance between two atom types, used to create neighbor list = (R_ij⁽²⁾)²
RLL: .5*RLL is the distance over which molecules can move without changing the neighbor list; RLL must be less than min(RB2)
RMAX: largest square distance at which interaction is computed = (R_ij⁽²⁾)²
RMAXLJ: largest square distance at which LJ interaction is computed; (2.5 σ_ij)² or 0 if ε=0
RNP: force on the atoms
RPP: gradient of potential energy for a pair
RSLJ: square distance within which atoms might interact through LJ; (sqrt(RMAXLJ)+RLL)² or 0 if ε=0
RSPL: start of the cubic LJ potential range; 0, or for C-H or Si-G pairs with nonzero ε set to 0.95 σ_ij
RSPLS: RSPL² or undefined if ε=0
RTABLE(KT1,KT2,1 + r_ij/DDtab): repulsive potential V_R function table
SIG(KT1,KT2): equilibrium distance σ of Lennard-Jones potential, converted internally to its square
SIGS: scale constant in the LJ potential induced by a semi-infinite solid
SPGC: G_C(θ) = Σ_p=0⁵ SPGC(p+1,IGC) cos^p(θ) where IGC depends on the interval in which cos(θ) is, with IGC=5 an alternate version for the last interval switched to for large enough N_i^(t)
SPGH: G_H(θ) = Σ_p=0⁵ SPGC(p+1,IGH) cos^p(θ) where IGH depends on the interval in which cos(θ) is
RSPL(KT1,KT2): set to 0.95 σ
SURF: NDIR-location of the surface of a neighboring continuum
TABDFC(KT1,KT2,1 + r_ij/DDtab): cut-off function's derivative table
TABFC(KT1,KT2,1 + r_ij/DDtab): cut-off function table
TEM: temperature
TLMN(1+N_i^t,1+N_j^t,N_ij^conj,term): tricubic splines for 2 A_ij from inter3dtors.d
TOTE: total potential energy
TTCONV: 2/(3 NP)
TTIME: total time
VOL: box volume (never used)
WW(K): f_ij(r_ij) for a pair of neighbors
XDB: α_ijk
XH(KT2,1+N^(H),1+N^(C)): values of H_CC and H_CH at the spline collocation points
XH1: derivative of XH wrt N^C
XH2: derivative of XH wrt N^H
XHC(I,1): 1 + sum of atom i's C neighbors' f_ij values
XHC(I,2): 1 + sum of atom i's H neighbors' f_ij values
XMASS: atomic mass, ordered by KT
XM(KT1,KT2): sqrt(XMM) or undefined if ε=0
XMM(KT1,KT2): minimum square distance for LJ interaction; 0, or (R_ij⁽²⁾)² for C-H or Si-G pairs with nonzero ε
XMMS(KT1,KT2): square radius of the inner range in which we will definitely not compute a LJ potential; corrected to max(0,sqrt(XMM)-RLL)² or undefined if ε=0
XMT: total mass
XQM: = 3.7; the value of N_i^(t) at which the evaluation of G_i for 0<θ<109 ends switching to a different polynomial

Variable Names (pibond.f)

Since subroutine pibond is probably the most impenetrable one, I include a separate list of local variables in pibond:

AA: d(V_ijkl^tor)/d(T2) keeping V_A, A_ij, T1, f_ik, f_jl, Z_ijk, Z_jil fixed
AAA1: V_ijkl^tor/f_ik f_jl Z_ijk Z_jil
AT2: - T1 d(V_ijkl^tor)/d(T1) keeping V_A, A_ij, T2, f_ik, f_jl, Z_ijk, Z_jil fixed
ATOR: A_ij(N_i^(t),N_j^(t),N_ij^conj) * 2 in the torsional potential
BIJ: B_ij
BJI: B_ji
BTOR: Σ_kl sin²(φ_ijkl) f_ik f_jl in the torsional potential
BTOT: bar B_ij * 2
CFUNI(NK): F(x_ik)
CFUNJ(NL): F(x_jl)
CJ: r_i - r_j = r_ij
CK: r_i - r_k = r_ik
CL: r_j - r_l = r_jl
CONJUG: N_ij^conj
CONK: Σ_k f_ik F(x_ik)
CONL: Σ_l f_jl F(x_jl)
COSK(NK): cos(θ_ijk)
COSL(NK): cos(θ_jil)
CRKX: x-component of r_ik x r_ij
CRKY: y-component of r_ik x r_ij
CRKZ: z-component of r_ik x r_ij
CRLX: x-component of r_ij x r_jl
CRLY: y-component of r_ij x r_jl
CRLZ: z-component of r_ij x r_jl
DALDIK: d(G_i)/d(N_i^(t))
DALDJL: d(G_j)/d(N_j^(t))
DATORC: d(A_ij)/d(N_ij^conj)
DATORI: d(A_ij)/d(N_i^(t))
DATORJ: d(A_ij)/d(N_j^(t))
DBDZI: d(B_ij)/d(b_ij)
DBDZJ: d(B_ji)/d(b_ji)
DBTOR*: never used
DCFUNI: f_ik d(F(x_ik))/d(x_ik)
DCFUNJ: f_jl d(F(x_jl))/d(x_jl)
DCTIJ(NK): d(cos(θ_ijk))/d(|r_ij|²/2)
DCTIK(NK): d(cos(θ_ijk))/d(|r_ik|²/2)
DCTIL(NL): d(cos(θ_jil))/d(|r_il|²/2)
DCTJI(NL): d(cos(θ_jil))/d(|r_ij|²/2)
DCTJK(NK): d(cos(θ_ijk))/d(|r_jk|²/2)
DCTJL(NL): d(cos(θ_jil))/d(|r_jl|²/2)
DEXNI(1): d(H_ij)/d(N_i^(C))
DEXNI(2): d(H_ij)/d(N_i^(H))
DEXNJ(1): d(H_ji)/d(N_j^(C))
DEXNJ(2): d(H_ji)/d(N_j^(H))
DFCK: d(FCK)/d(|r_ik|)
DFCL: d(FCL)/d(|r_jl|)
DGDTHET: d(G_i)/d(cos θ_ijk) respectively d(G_j)/d(cos θ_jil)
DRADI: d(bar B_ij)/d(N_i^(t)) * 2
DRADJ: d(bar B_ij)/d(N_j^(t)) * 2
DRDC: d(bar B_ij)/d(N_ij^conj) * 2
DT1DIJ: d(T1)/d(|r_ij|²/2)/T1
DT1DIK: d(T1)/d(|r_ik|²/2)/T1
DT1DIL: d(T1)/d(|r_il|²/2)/T1
DT1DJK: d(T1)/d(|r_jk|²/2)/T1
DT1DJL: d(T1)/d(|r_jl|²/2)/T1
DT2DIJ(.): d(T2)/d(r_ij(.))
DT2DIK(.): d(T2)/d(r_ik(.))
DT2DJL(.): d(T2)/d(r_jl(.))
EXNIJ: H_ij
EXNJI: H_ji
EXX: e_ijk = exp(α_ijk[(r_ij - R_ij^e) - (r_ik-R_ik^e)]) respectively e_jil
FCK: f_ik (modified for H to R_ij¹ = 1.3, R_ij² = 1.6)
FCL: f_jl (modified for H to R_ij¹ = 1.3, R_ij² = 1.6)
GANGLE1: alternate value of GANGLE used for high N_i^(t) or N_j^(t)
GANGLE: G_i(θ_ijk) respectively G_j(θ_jil)
IG: G_i-interpolation interval index
J: pair list index of pairs with i the first atom
JN: j value, with iKI: KTYPE(i)
KIKJ: KI+KJ
KJ: KTYPE(j)
KK: KTYPE(k)
KN: second neighbor k≠j to i
L: pair list index of pairs with j the first atom
NK: counts number of secondary neighbors k of i
NL: counts number of secondary neighbors l of j
QI: N_i^(t) = N_i^(C) + N_i^(H) (does not include j)
QJ: N_j^(t) = N_j^(C) + N_j^(H) (does not include i)
RCK: |r_ik|
RCL: |r_jl|
RSQ2: |r_jk|² respectively |r_il|²
RSQ3: |r_ik|² respectively |r_jl|²
RSQIJ: |r_ij|²
S3: |r_ik| respectively |r_jl|
SDALIK: d(Σ_k G_i f_ij e_ijk)/d(N_i^(t))
SDALJL: d(Σ_l G_j f_ji e_jil)/d(N_j^(t))
SIJ: |r_ij|
SINK(NK): sin(θ_ijk)
SINL(NK): sin(θ_jil)
SSUMK: Σ_k G_i f_ij e_ijk
SSUML: Σ_l G_j f_ji e_jil
T1: |r_ik x r_ij| |r_jl x r_ji|
T2: (r_ik x r_ij) . (r_jl x r_ji); negative if k and l point to the same side of i-j
VATT: V_A / 2
VDBDI: d(V)/d(b_ij)
VDBDJ: d(V)/d(b_ji)
VDRDC: d(V)/d(N_ij^conj)
VDRDI: d(V)/d(N_i^(t))
VDRDJ: d(V)/d(N_j^(t))
XK: r_j - r_k = r_jk
XL: r_i - r_l = r_il
XNI(1): = 1 + N_i^(C) = 1 + atom i's C neighbors f_ik values excluding j
XNI(2): = 1 + N_i^(H) = 1 + atom i's C neighbors f_ik values excluding j
XNJ(1): = 1 + N_j^(C) = 1 + atom j's C neighbors f_jl values excluding i
XNJ(2): = 1 + N_j^(H) = 1 + atom j's C neighbors f_jl values excluding i
XNT1: N_i^(t)+1
XNT2: N_j^(t)+1
XSIJ: d(Σ_k G_i f_ij e_ijk)/d(|r_ij|²/2)
XSIK(NK): d(G_i f_ij e_ijk)/d(|r_ik|²/2)
XSIL(NL): d(G_j f_ji e_jil)/d(|r_il|²/2)
XSJI: d(Σ_l G_j f_ji e_jil)/d(|r_ij|²/2)
XSJK(NK): d(G_i f_ij e_ijk)/d(|r_jk|²/2)
XSJL(NL): d(G_j f_ji e_jil)/d(|r_jl|²/2)
XX: x_ik respectively
ZF1: Z_ijk
ZF1DCT: d(Z_ijk)/d(cos(θ_ijk))
ZF2: Z_jil
ZF2DCT: d(Z_jil)/d(cos(θ_jil))

Return