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 A2/ps2. 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 VA in param.f are
grossly incompatible with those in Brenner's paper.
- Moreover, it appears that the potential VA 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.
- Ni(C), Ni(H),
Ni(t) are really Nij values since j
is excluded in the sums.
- Gi for carbon has been made dependent on
Ni(t): the program switches smoothly to a
different representation for Gi between 4.2 and 4.7
neighbors to the i-atom (including j), i.e. it switches in the range
3.2 < Ni(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 Nijconj
- The program includes a "torsional interaction" for CC bonds. The
added torsional potential is:
Vtor =
VA(|rij|)
Aij(Ni(t),Nj(t),Nijconj)
Σkl sin2(φijkl) fik
fjl
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 Vtor 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/BondOrder/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:
Vtor = VA(|rij|)
Aij(Ni(t),Nj(t),Nijconj)
Σkl sin2(φijkl) fik
fjl Zijk Zjil
where in the smoothing range
Zijk = 1 + cos(π (sin(θijk)-SINHI)/(SINHI-SINLO))
Zjil = 1 + cos(π (sin(θjil)-SINHI)/(SINHI-SINLO))
- The program includes a Lennard-Jones potential:
VLJ=
4 εij [(σij/rij)12-
(σij/rij)6]
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
rij > 2.5 σij and in the range
rij < Rij(2) in which the bond order
potential is nonzero.
In the range Rij(2) < rij
< 0.95 σij;, the potential
is furthermore replaced by a cubic which vanishes quadratically at
Rij(2) 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 VA and VR from the Brenner paper..
- Param
- C and H potentials.
- Pibond
- Finds bonding forces involving C and H atoms.
- Pred
- Nordsieck scheme prediction of zn+1 from zn 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:
- eijk = exp(αijk [(rij -
Rije) -
(rik-Rike)])
- fij = fij(|rij|)
- Vijkltor = VA(|rij|)
Aij(Ni(t),Nj(t),Nijconj)
sin2(φijkl) fik fjl (a
term in the torsional potential)
- bij is defined by the relation Bij
= bij-δ
- AD(KT1,KT2): Dije Sij exp(AXL Rije) /(1 - Sij) value 1
- ATABLE(KT1,KT2,1 + rij/DDtab): attractive potential VA/2 function table
- ATT: = 3.2; the value of Ni(t) at which the evaluation of Gi for 0<θ<109 starts to smoothly switch to a different polynomial
- AXL(KT1,KT2): sqrt(2/Sij) βij value 1
- BD(KT1,KT2): Dije Sij exp(AXL Rije) /(1 - Sij) value 2
- BXL(KT1,KT2): sqrt(2/Sij) βij value 2
- CD(KT1,KT2): Dije Sij exp(AXL Rije) /(1 - Sij) value 3
- CLM(KT2,1+Ni(H),1+Ni(C),term): bicubic splines HCC and HCH in inter2div.d
- CLMN(1,1+Ni(t),1+Nj(t),Nijconj,term): tricubic splines 2 FCC from inter3div.d
- CLMN(2,1+Ni(t),1+Nj(t),Nijconj,term): tricubic splines 2 FCH from inter3dch.d
- CLMN(3,1+Ni(t),1+Nj(t),Nijconj,term): tricubic splines 2 FHH from inter3dh.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/Sij) βij value 3
- DATABLE(KT1,KT2,1 + rij/DDtab): attractive potential's derivative * 2 table
- DD(KT1,KT2): Dije exp(DXL Rije) /(1 - Sij)
- DDTAB(KT1,KT2): rij-spacing of the potential tables
- DELLJ: table spacing for Lennard-Jones potential, 0.001 A
- DELTA: time step
- DELTSQ: DELTA2/2
- DEXX1(K): derivative of the attractive potential VA / 2 of a pair K
- DRTABLE(KT1,KT2,1 + rij/DDtab): repulsive potential's derivative table
- DWW(K): derivative of fij(rij) for a pair of neighbors
- DXL(KT1,KT2): sqrt(2 Sij) βij
- EATOM(I): potential energy of an atom
- ECONV: Conversion factor: 103.6434 eV per (amu A2/fs2)
- 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 VA / 2 of a pair of neighbors K
- I2D: potential version number in inter2div.d
- I3D: potential version number in inter3div.d, inter3dh.d, inter3dch.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 Gi 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 Gi 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 writedata2
- 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 Rij(1) at the start of the bond order potential decay
- RB2: distance Rij(2) at which the bond order potential is gone
- RCOR: distance between a pair of interacting atoms
- REG: exp(αijk [Rike - Rije])
- RLIST: square interaction distance between two atom types, used to create neighbor list = (Rij(2))2
- 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 = (Rij(2))2
- RMAXLJ: largest square distance at which LJ interaction is computed; (2.5 σij)2 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)2 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: RSPL2 or undefined if ε=0
- RTABLE(KT1,KT2,1 + rij/DDtab): repulsive potential VR 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: GC(θ) = Σp=05 SPGC(p+1,IGC) cosp(θ) 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 Ni(t)
- SPGH: GH(θ) = Σp=05 SPGC(p+1,IGH) cosp(θ) 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 + rij/DDtab): cut-off function's derivative table
- TABFC(KT1,KT2,1 + rij/DDtab): cut-off function table
- TEM: temperature
- TLMN(1+Nit,1+Njt,Nijconj,term): tricubic splines for 2 Aij from inter3dtors.d
- TOTE: total potential energy
- TTCONV: 2/(3 NP)
- TTIME: total time
- VOL: box volume (never used)
- WW(K): fij(rij) for a pair of neighbors
- XDB: αijk
- XH(KT2,1+N(H),1+N(C)): values of HCC and HCH at the spline collocation points
- XH1: derivative of XH wrt NC
- XH2: derivative of XH wrt NH
- XHC(I,1): 1 + sum of atom i's C neighbors' fij values
- XHC(I,2): 1 + sum of atom i's H neighbors' fij 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 (Rij(2))2 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)2 or undefined if ε=0
- XMT: total mass
- XQM: = 3.7; the value of Ni(t) at which the evaluation of Gi 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(Vijkltor)/d(T2) keeping VA, Aij, T1, fik, fjl, Zijk, Zjil fixed
- AAA1: Vijkltor/fik fjl Zijk Zjil
- AT2: - T1 d(Vijkltor)/d(T1) keeping VA, Aij, T2, fik, fjl, Zijk, Zjil fixed
- ATOR: Aij(Ni(t),Nj(t),Nijconj) * 2 in the torsional potential
- BIJ: Bij
- BJI: Bji
- BTOR: Σkl sin2(φijkl) fik fjl in the torsional potential
- BTOT: bar Bij * 2
- CFUNI(NK): F(xik)
- CFUNJ(NL): F(xjl)
- CJ: ri - rj = rij
- CK: ri - rk = rik
- CL: rj - rl = rjl
- CONJUG: Nijconj
- CONK: Σk fik F(xik)
- CONL: Σl fjl F(xjl)
- COSK(NK): cos(θijk)
- COSL(NK): cos(θjil)
- CRKX: x-component of rik x rij
- CRKY: y-component of rik x rij
- CRKZ: z-component of rik x rij
- CRLX: x-component of rij x rjl
- CRLY: y-component of rij x rjl
- CRLZ: z-component of rij x rjl
- DALDIK: d(Gi)/d(Ni(t))
- DALDJL: d(Gj)/d(Nj(t))
- DATORC: d(Aij)/d(Nijconj)
- DATORI: d(Aij)/d(Ni(t))
- DATORJ: d(Aij)/d(Nj(t))
- DBDZI: d(Bij)/d(bij)
- DBDZJ: d(Bji)/d(bji)
- DBTOR*: never used
- DCFUNI: fik d(F(xik))/d(xik)
- DCFUNJ: fjl d(F(xjl))/d(xjl)
- DCTIJ(NK): d(cos(θijk))/d(|rij|2/2)
- DCTIK(NK): d(cos(θijk))/d(|rik|2/2)
- DCTIL(NL): d(cos(θjil))/d(|ril|2/2)
- DCTJI(NL): d(cos(θjil))/d(|rij|2/2)
- DCTJK(NK): d(cos(θijk))/d(|rjk|2/2)
- DCTJL(NL): d(cos(θjil))/d(|rjl|2/2)
- DEXNI(1): d(Hij)/d(Ni(C))
- DEXNI(2): d(Hij)/d(Ni(H))
- DEXNJ(1): d(Hji)/d(Nj(C))
- DEXNJ(2): d(Hji)/d(Nj(H))
- DFCK: d(FCK)/d(|rik|)
- DFCL: d(FCL)/d(|rjl|)
- DGDTHET: d(Gi)/d(cos θijk) respectively d(Gj)/d(cos θjil)
- DRADI: d(bar Bij)/d(Ni(t)) * 2
- DRADJ: d(bar Bij)/d(Nj(t)) * 2
- DRDC: d(bar Bij)/d(Nijconj) * 2
- DT1DIJ: d(T1)/d(|rij|2/2)/T1
- DT1DIK: d(T1)/d(|rik|2/2)/T1
- DT1DIL: d(T1)/d(|ril|2/2)/T1
- DT1DJK: d(T1)/d(|rjk|2/2)/T1
- DT1DJL: d(T1)/d(|rjl|2/2)/T1
- DT2DIJ(.): d(T2)/d(rij(.))
- DT2DIK(.): d(T2)/d(rik(.))
- DT2DJL(.): d(T2)/d(rjl(.))
- EXNIJ: Hij
- EXNJI: Hji
- EXX: eijk = exp(αijk[(rij - Rije) - (rik-Rike)]) respectively ejil
- FCK: fik (modified for H to Rij1 = 1.3, Rij2 = 1.6)
- FCL: fjl (modified for H to Rij1 = 1.3, Rij2 = 1.6)
- GANGLE1: alternate value of GANGLE used for high Ni(t) or Nj(t)
- GANGLE: Gi(θijk) respectively Gj(θjil)
- IG: Gi-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: Ni(t) = Ni(C) + Ni(H) (does not include j)
- QJ: Nj(t) = Nj(C) + Nj(H) (does not include i)
- RCK: |rik|
- RCL: |rjl|
- RSQ2: |rjk|2 respectively |ril|2
- RSQ3: |rik|2 respectively |rjl|2
- RSQIJ: |rij|2
- S3: |rik| respectively |rjl|
- SDALIK: d(Σk Gi fij eijk)/d(Ni(t))
- SDALJL: d(Σl Gj fji ejil)/d(Nj(t))
- SIJ: |rij|
- SINK(NK): sin(θijk)
- SINL(NK): sin(θjil)
- SSUMK: Σk Gi fij eijk
- SSUML: Σl Gj fji ejil
- T1: |rik x rij| |rjl x rji|
- T2: (rik x rij) . (rjl x rji); negative if k and l point to the same side of i-j
- VATT: VA / 2
- VDBDI: d(V)/d(bij)
- VDBDJ: d(V)/d(bji)
- VDRDC: d(V)/d(Nijconj)
- VDRDI: d(V)/d(Ni(t))
- VDRDJ: d(V)/d(Nj(t))
- XK: rj - rk = rjk
- XL: ri - rl = ril
- XNI(1): = 1 + Ni(C) = 1 + atom i's C neighbors fik values excluding j
- XNI(2): = 1 + Ni(H) = 1 + atom i's C neighbors fik values excluding j
- XNJ(1): = 1 + Nj(C) = 1 + atom j's C neighbors fjl values excluding i
- XNJ(2): = 1 + Nj(H) = 1 + atom j's C neighbors fjl values excluding i
- XNT1: Ni(t)+1
- XNT2: Nj(t)+1
- XSIJ: d(Σk Gi fij eijk)/d(|rij|2/2)
- XSIK(NK): d(Gi fij eijk)/d(|rik|2/2)
- XSIL(NL): d(Gj fji ejil)/d(|ril|2/2)
- XSJI: d(Σl Gj fji ejil)/d(|rij|2/2)
- XSJK(NK): d(Gi fij eijk)/d(|rjk|2/2)
- XSJL(NL): d(Gj fji ejil)/d(|rjl|2/2)
- XX: xik respectively
- ZF1: Zijk
- ZF1DCT: d(Zijk)/d(cos(θijk))
- ZF2: Zjil
- ZF2DCT: d(Zjil)/d(cos(θjil))
Return