GULP help file

cell => don't output the unit cell parameters for vectors coordinates => don't output the coordinates of the atoms structure => don't output the structure (i.e. cell & coordinates) potentials => don't output the interatomic potentials derivatives => don't output the final derivatives

nebspring 10.0

For a spring constant that varies between 10.0 and 1.0 the input would be;

nebspring vary 10.0 1.0

1 => central finite difference between replicas 2 => bisect vectors 3 => new algorithm favouring higher energy replica neighbour

plane_lj 10 4 C core -0.75 1000.0 14.0 0.0 12.0

E = sum(i) 1/2 K*[(x_i - x)**2 + (y_i - y)**2 + (z_i - z)**2]

where the sum is over all atoms, i.

F = forceA*cos[2*pi*(t*forceB + forceC)]

where t is the time within the simulation. Example 1: 1 z 3.5 0.5 0.5 3 z 3.5 0.5 0.0 Example 2: 1 z -3.5 0.5 0.0 3 z 3.5 0.5 0.0 Note: both of the other examples would cause the atoms to oscillate exactly out of phase with each other and are equivalent. The maximum force experienced would be +/- 3.5 eV/Angstrom

all - copy all of the below items conditions - copy temperature / pressure data md - copy all molecular dynamics related parameters solvent - copy solvent related data structure - copy structure related data

The parameter, configuration_no, gives the configuration number to be copied when creating the present new structure. This option is useful when wanting to run a sequence of calculations on the same structure with different conditions. Note: The coordinates cannot be presently overwritten.

Example :

omega 500.0 100.0 5

the above would calculate properties at a frequencies of 500, 600, 700, 800, 900 and 1000 cm-1. A no_of_steps = 0 implies a calculation at a single frequency.

NOTE : This option only applies to a gamma point calculation

NOTE : This option only applies to a gamma point calculation

By default the directions are assumed to be in Cartesian space unless the "frac" sub-option is given in which case they will be taken relative to the crystal axes.

The format of the file is as follows. Each configuration is written as follows:

write(31)numat write(31)energy write(31)(atomicno(i),i=1,numat) write(31)(atomictype(i),i=1,numat) write(31)(xcoordinate(i),i=1,numat) write(31)(ycoordinate(i),i=1,numat) write(31)(zcoordinate(i),i=1,numat)

where the variables are as follows :

integer(i4) numat = total number of atoms integer(i4) atomicno = atomic number of atom (1-maxele) integer(i4) atomictype = atomic type of atom (0-9999) real(dp) energy = total energy of system (eV) real(dp) xcoordinate = Cartesian coord in X direction (Angs) real(dp) ycoordinate = Cartesian coord in Y direction (Angs) real(dp) zcoordinate = Cartesian coord in Z direction (Angs)

and the datatypes are given by :

integer, parameter :: i4 = selected_int_kind(9) integer, parameter :: dp = kind(1.0d0)

gcmcexistingmolecules 2 3

Where there are many molecules a range can also be specified as follows:

gcmcexistingmolecules 4 to 8

e.g. lowest 4 9

e.g. keyword opti conp prop

NB Electronegativity equalisation schemes should NOT be used in combination with Coulomb subtraction of any form otherwise the calculation of the charges and the energy will not be self-consistent. This leads to all derivatives being incorrect as dE/dQ is no longer zero. If Coulomb subtracted potentials are to be used then charges must be calculated for the initial geometry and then frozen.

Note: it is important to investigate the effect of qeqradius on the degree of convergence and CPU time.

Note: Phonon calculations are limited with the QEq keyword to the gamma point without the LO/TO splitting.

Note: Only region 1 atoms are included in QEq and fixed charges will be taken for region 2.

E = -A.[1+B*((r/r0)-1)].exp(-B*((r/r0)-1))

The densities at each site are determined by the "eam_density" option and the functional dependance on the total density by "eam_functional".

Note that although for simplicity the manybody potential appears as part of the two-body potentials, it is infact many body at short-range, but tends to an effective pair potential at long range. Also note that this potential type is NOT compatible with <intra/inter/molmec> directives.

if r < d:

E = a14*(r-d)**4 + a13*(r-d)**3 + a12*(r-d)**2 + a11*(r-d) + a10

if d < r =< rmax:

E = a26*(r-d)**6 + a25*(r-d)**5 + a24*(r-d)**4 + a23*(r-d)**3 + a22*(r-d)**2 + a21*(r-d) + a20

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction.

Square_root:

E = - sum(i) A(i)*(rho(i))**1/2

this is the most common functional, as used in the Sutton-Chen potential

Power:

E = - sum(i) A(i)*(rho(i))**1/n

this is just a generalisation of the above case

Banerjea_smith:

E = - sum(i) F0 [1-ln(r)/n]*r**1/n + F1*r

where r = rho(i)/rho0(i)

this is the functional of Banerjea and Smith (Phys. Rev. B, 37, 6632 (1988)) - note that in this case that are atom dependant parameters also to be specified (F0, F1, rho0) where rho0 is the electron density at equilibrium.

Johnson:

E = - sum(i) F0 [1-ln(x)]*x + F1*y

where x = (rho(i)/rho0(i))**(alpha/beta) and y = (rho(i)/rho0(i))**(gamma/beta)

this functional is similar to that of Banerjea & Smith and is due to Johnson (PRB, 39, 12554 (1989)).

Glue:

if rho < rho1 E = c1_4*(rho-rho1)**4 + c1_3*(rho-rho1)**3 + c1_2*(rho-rho1)**2 + c1_1*(rho-rho1) + c1_0

if rho1 =< rho < rho2 E = c2_4*(rho-rho2)**4 + c2_3*(rho-rho2)**3 + c2_2*(rho-rho2)**2 + c2_1*(rho-rho2) + c2_0

if rho2 =< rho E = c3_3*(rho-rho2)**3 + c3_2*(rho-rho2)**2 + c3_1*(rho-rho2) + c3_0

this is the functional from Ercolessi et al, Phil. Mag. A, 58, 213 (1988). NB: At present there are no fitting parameters since there are constraints on the coefficients to ensure a smooth and continuous function.

Foiles:

E = sum(i) F0*rho(i)**2 + F1*rho(i) + F2*(rho(i)**(5/3)/(F3 + rho(i)))

Mei-Davenport:

E = sum(i) - Ec*[1-(alpha/beta)*ln(rho(i))]*rho(i)**(alpha/beta) + sum(m=1->3) 0.5*phi0*s_m*exp(-(sqrt(m)-1)*gamma)* [1+(sqrt(m)-1)*delta-sqrt(m)*delta*ln(rho(i))/beta]* rho(i)**(sqrt(m)*gamma/beta)

Power Law:

rho(i) = C*rij**(-n)

e.g. eam_density power 6 Ni core 729.7

Exponential:

rho(i) = A*(rij**n)*exp(-B(rij-r0))

e.g. eam_density exponential 0 Ni core 500.0 4.0 3.52

Gaussian:

rho(i) = A*(rij**n)*exp(-B(rij-r0)**2)

e.g. eam_density gaussian 2 Ni core 400.0 3.0 3.52

Quadratic:

rho(i) = A*(rij-r0)**2 if r < r0, else = 0

Cubic:

rho(i) = A*(rij-r0)**3 if r < r0, else = 0

Quartic:

rho(i) = A*(rij-r0)**4 if r < r0, else = 0

Voter:

rho(i) = A*r**6*(exp(-beta*r) + 2**9*exp(-2*beta*r))

eVoter:

rho(i) = A*r**6*(exp(-beta*r) + 2**9*exp(-2*beta*r))*exp(-1/(rmax - r))

Mei-Davenport:

rho(i) = sum(l=0->5) (c_l/12)*(r/r0)**l

In the above, rmax is the cutoff for the potential from the manybody option that controls the range of the density.

Glue:

if r < r1 rho(i) = b1_3*(r-r1)**3 + b1_2*(r-r1)**2 + b1_1*(r-r1) + b1_0

if r1 =< r < r2 rho(i) = b2_3*(r-r2)**3 + b2_2*(r-r2)**2 + b2_1*(r-r2) + b2_0

if r2 =< r < rm rho(i) = b3_3*(r-rm)**3 + b3_2*(r-rm)**2 + b3_1*(r-rm) + b3_0

Note that the cut-offs are set by the manybody potential except for the glue model where the maximum cutoff, rm, is a parameter of the model.

rho'(i) = scale(i)*rho(i)

F'(sum(rho(i))) = F(sum(rho(i))/scale(i))

where F is the unscale function of the density.

The form of the potential is

E = 2.0*g*r**6*(exp(-beta*r) + 2**9*exp(-2*beta*r))

The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance.

E = 2*pi*D**2 / 3.V

where D is the dipole per unit cell and V is the volume. By default the Ewald sum assumes that there is no dipole moment across the crystal, while this term is applicable to cases where there is a permenant dipole. Note that the dipole is ambiguous in many cases because it depends on the termination of the crystal at the surface and hence this correction should be used with care! Note that defect calculations cannot be performed when this term is present.

shellmass Al 0.25 O 0.12

Example:

name alumina

translate x y z nstep

where x, y, z are the components of the vector between the the initial and final positions of the atoms. If the system is three-dimensional then x, y and z are assumed to be in fractional units. If the system is a cluster then they are assumed to be in angstroms. nstep is the number of points to be sampled along the translation vector (this leads to nstep+1 calculations including the first and last points). The subset of atoms to which the translation is to be applied is defined by adding a "T" flag to the end of the coordinate lines of these atoms.

NOTE: cmm cannot be used in conjunction with EEM or QEq at the moment.

Valid minimisers are: bfgs (default) rfo (rational function optimisation) unit (bfgs, but starting from a unit hessian) nume (bfgs, but starting from numerical diagonal hessian) conj (conjugate gradients)

Examples:

switch rfo cycle 10

=> Change to rfo method after 10 cycles of optimisation. Note that the order of words after switch is irrelevant and additional words can be inserted to make the line more readable. For example this line could be written as:

switch to rfo method after 10 cycles of optimisation

switch conj gnorm 0.123

=> Change to conjugate gradient method when Gnorm is less than 0.123. Note that if cycle and gnorm are omitted then GULP assumes gnorm as a default.

If "reverse" is specified then order is: V x y z <weight>

buck Si core O shel 1280.0 0.3 0.0 0.0 12.0 1 0 0

or

buck Si core O shel & 1280.0 0.3 0.0 0.0 12.0 & 1 0 0

atom_symbol/atomic_number <core/shel> K r0 <2 x flags>

if exponential form:

atom_symbol/atomic_number <core/shel> K rho r0 <3 x flags>

if single_exponential form:

atom_symbol/atomic_number <core/shel> K rho r0 <3 x flags>

E(bs) = 1/2 * K * (r - r0)**2

or the constants of the exponential restoring term:

E(bs) = K * [exp(rho*(r-r0)) + exp(-rho*(r-r0))]

or the constants of the single_exponential restoring term:

E(bs) = K * exp(rho*(r-r0))

momentum_correct 100 10

mechanical: This implies that the interactions within the MM regions and between the QM and MM regions are calculated, but not within the QM region.

electronic: This implies the same as mechanical embedding, except that the Coulomb interactions between the QM and MM regions are also excluded since these are assumed to be accounted for by the QM portion of the calculation.

NB: At present there are restrictions on which facilities will work in combination with this option. Only 0-D calculations with up to first derivatives are generally enabled for all models.

By default the program uses method 4 for charged defects as this offers the best compromise between accuracy and computational effort. It is recommended that a comparison is made with mode2a=1 as a check. For optimisation the most efficient approach could be to optimise first with a higher mode and then restart in mode 1.

1 => xx 2 => yy 3 => zz 4 => yz 5 => xz 6 => xy

Note that stresses are dependent on the cell orientation and so take care when specifying.

NB Electronegativity equalisation schemes should NOT be used in combination with Coulomb subtraction of any form otherwise the calculation of the charges and the energy will not be self-consistent. This leads to all derivatives being incorrect as dE/dQ is no longer zero. If Coulomb subtracted potentials are to be used then charges must be calculated for the initial geometry and then frozen.

Note: Phonon calculations are limited with the EEM keyword to the gamma point without the LO/TO splitting.

Note: Only region 1 atoms are included in EEM and fixed charges will be taken for region 2.

qincrement O core C core -0.24

NB: The order of the atoms is important in that the sign of delta changes if the order is reversed. Consequently if atom1 and atom2 are of the same type then delta must be zero. Hence, the value of delta is applied with the input sign to atom1 and the negative of this to atom2.

e.g. broaden 0.3

The expression used for the broadening is:

w(f) = b/[pi*(1 + (b*(f - f0))**2)]

where w is the fractional weight at frequency f of a peak at frequency f0. b is the broadening factor and is the inverse of the half-width at half-maximum parameter.

Note, the smaller the factor, the greater the broadening. The range of the plot is automatically extended when broadening is turned on to allow the highest frequency peak to decay to 1/100th of its original magnitude.

cart region 3 rigid z

If a "T" is specified then the atom is marked for the translate option If a "%" is specified then the atom is part of a growth slice if it is in region 1 of a surface calculation.

sfrac region 3 rigid z

If a "T" is specified then the atom is marked for the translate option If a "%" is specified then the atom is part of a growth slice if it is in region 1 of a surface calculation.

The form of the Voter taper is:

E_smooth(r) = E(r) - E(rcut) + (rcut/m)*(1 - (r/rcut)**m)*(dE(rcut)/dR)

The form of the exponential taper is:

E_smooth(r) = E(r)*f_cut(r)

where f_cut(r) = exp(-1/(rcut-r)) for r < rcut and = 0 for r >= rcut

Note that in the case of the Voter and exponential tapers there is no taper range since it applies over the whole range.

The form of the MDF (Mei-Davenport-Fernando) taper is:

E_smooth(r) = E(r)*f_cut(r)

where f_cut(r) = 1.0 for r < r_m and = 0 for r >= rcut. In between it takes the value f_cut(r) = ((1 - x)**3)*(1+3x+6x**2), where x = (r-r_m)/(rcut-r_m). Here r_m is the start of the taper range.

Atom 1: Q = sum [ - deltaQ*H(r) ] Atom 2: Q = sum [ + deltaQ*H(r) ]

Here the sum is over all atoms of the approriate type within the cutoff radius, Rmax. H(r) is taper function that acts between Rmin and Rmax. In the original paper of Jiang and Brown a cosine taper is used and this is the default. However, the sub-option allows the sine taper of Watanabe et al (Jpn. J. Appl. Phys. 38, L367 (1999)) to be used which has better behaviour for the second derivatives at the cutoff.

E = Kq*exp(-rho/(q - q0)) if q > q0 for q0 > 0

or

E = Kq*exp(rho/(q - q0)) if q > q0 for q0 < 0

else E = 0

or

botwobody combine atom1 atom2 chiR chiA <2*flags>

E = f(r)[A.exp(-za.r).(BOr) - B.exp(-zb.r).(BOa)]

where f(r) is a taper function that smooths the decay from rtaper to rmax, BOr and BOa are the bond orders for the repulsive and attractive terms, respectively. These terms are set by separate options. If the combine sub-option is specified then the parameters are generated using Tersoff's combination rules from the values for the corresponding element's self-self interaction. The chi values then scale the repulsive and attractive terms.

or

borepulsive theta atom1 alpha m n lambda c d h <6*flags> NB The value of m cannot be fitted and therefore there is no flag

BOr = (1 + (alpha*zeta)**n)**(-1/2n)

For case without "theta" sub-option;

zeta = Sum(rik) [f(rik).exp(lambda**m.(rij-rik)**m)]

else;

zeta = Sum(rik) [f(rik).g(theta).exp(lambda**m.(rij-rik)**m)]

f(rik) = cosine taper function

g(theta) = 1 + (c/d)**2 - c**2/[d**2 + (h - cos(theta))**2]

Note that the cut-offs are derived from the values for the corresponding botwobody potential.

or

boattractive theta atom1 alpha m n lambda c d h <6*flags> NB The value of m cannot be fitted and therefore there is no flag

BOa = (1 + (alpha*zeta)**n)**(-1/2n)

For case without "theta" sub-option;

zeta = Sum(rik) [f(rik).exp(lambda**m.(rij-rik)**m)]

else;

zeta = Sum(rik) [f(rik).g(theta).exp(lambda**m.(rij-rik)**m)]

f(rik) = cosine taper function

g(theta) = 1 + (c/d)**2 - c**2/[d**2 + (h - cos(theta))**2]

Note that the cut-offs are derived from the values for the corresponding botwobody potential.

E = A.exp(-r/rho) - C/r**6

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential, cutoffs are omitted from input. The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance.

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential, cutoffs are omitted from input.

U = [Q1.Q2 - q1.q2]g(r)/r

where q1 & q2 are the regular charges specified elsewhere, and g(r) is:

Form1: g(r) = (1+z.r)exp(-2z.r)

Form2: g(r) = (1 + 11(z.r)/8 + 3(z.r)**2/4 + (z.r)**3/6)exp(-2z.r)

where z is used as a short-hand for zeta

E = A/r**m - B/r**n

The exponents are by default 12 and 6, but these can be changed by specifying values after the option word lennard. By specifying epsilon after lennard this means that the input is in terms of epsilon and sigma instead of A and B.

E = epsilon*(c1*(sigma/r)**m - c2*(sigma/r)**n)

Specifying "zero" allows the user to chose between sigma defined as the potential energy minimum distance (default) or the distance at which the potential energy goes to zero.

r = sigma => E = epsilon

c1 = (n/(m-n)) c2 = (m/(m-n))

r = sigma => E = zero

c1 = (n/(m-n))*(m/n)**(m/(m-n)) c2 = (m/(m-n))*(m/n)**(n/(m-n))

If combine is specified as well as epsilon, then combination rules are used to obtain the epsilon and sigma parameters based on the values for individual species. The combination rules are as follows:

epsilon = 2*sqrt(e1*e2)(s1**3.s2**3)/(s1**6+s2**6) sigma = ((s1**6+s2**6)/2)**1/6

where e1, e2, s1 and s2 are the atom related parameters.

Alternatively, if geometric combination rules are specified then this becomes:

epsilon = sqrt(e1*e2) sigma = (s1+s2)/2

If combine is specified for a standard lennard-jones potential then the species values given in the "atomab" command are used to obtain potential parameters using the following rules:

Aij = sqrt(Ai.Aj) Bij = sqrt(Bi.Bj)

If esff is specified then this implies that the ESFF combination rules will be used to calculate the potential parameters. This automatically implies "combine" and that the functional form is 9/6. The parameters are then obtained from epsilon and sigma as follows:

Aij = Ai.Bj + Aj.Bi Bij = 3 * Bi.Bj

where Ai = sqrt(epsilon)*sigma**6 Bi = sqrt(epsilon)*sigma**3

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction.

If "all" is specified as a sub-option, then Lennard-Jones potentials will be automatically created and added for all combinations of A/B or epsilon/sigma, as appropriate, based on the values that have been specified up to that point. In this case, the line following should only have the rmin and rmax to be used for all the potentials.

e.g. lennard 12 6 combine all 0.0 12.0

When specified as bonded potential, cutoffs are omitted from input. The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance.

E = - A.exp(-b(r-r0)**2)

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input.

E = -D*exp(-a*(r-r0)**2/(2r))

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input. The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance.

atom1 and atom2 may be specified either by atomic number or symbol which in the latter case can be followed by a species type. If no species type is given then type is assumed to be core.

E = De.((1-exp(-a(r-r0)))**2 - 1.0) - coul.qi.qj/r

If "etaper" is specified then the potential is scaled by:

exp(-1/(rmax-r))

to ensure that the potential (excluding the Coulomb subtraction) goes to zero at the cutoff. In this case, the energy/gradient options cannot be used. Note that the potential can also be written as:

E = De.(exp(-2a(r-r0)) - 2exp(-a(r-r0))) - coul.qi.qj/r

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input. The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance.

E = a/(1+exp(b(r-r0)))

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input. The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance.

E_ij = A.erf(r/alpha).erfc(r/beta)/r

E_ij = A.erf(r/alpha)/r

E_ij = A.erfc(r/beta)/r

E_ij = - phi0*(1 + delta*(r/r0 - 1))*exp(-gamma*(r/r0 - 1))

E = A/(r + r0)**m - B/(r + r0)**n

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input.

E = [qi.qj/r].f(r) + C.(1-f(r)), where f(r) = polynomial taper

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input.

E = [qi.qj/r].erfc(r/rho)

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input. The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance. NOTE : The keyword "noelectrostatics" should be included when using this potential to turn off the normal Ewald sum, otherwise the Coulomb interaction will be double counted.

E = [qi.qj/r**2]

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then potential will only act between bonded atoms. x12 => exclude 1-2 interaction and x13 => exclude 1-2 and 1-3 interaction. When specified as bonded potential cutoffs are omitted from input. The words energy or gradient may be given on the first line resulting in energy or gradient offsets being applied such that the specified quantity goes to zero at the cutoff distance. NOTE : The keyword "noelectrostatics" should be included when using this potential to turn off the normal Ewald sum, otherwise the Coulomb interaction will be double counted. NOTE : This expression is used for the Coulomb term in Accelrys's implementation of Dreiding

Eij = [qi.qj/r].erfc(eta*r) - lim{r->rmax}([qi.qj/rmax].erfc(eta*rmax))

Ei = - (erfc(eta*rmax)/(2*rmax) + eta/sqrt(pi))*qi**2

Note: it is up to the user to set the eta and rmax values such that the Ewald limit is recovered. At present the Wolf sum cannot be used with a defect calculation.

space P 1 21/A 1

Exponentially decaying form is also available: E(three) = 1/2 * k * (theta-theta0)**2.exp(-r12/rho1).exp(-r13/rho2) Format: three exponential atom1 atom2 atom3 k theta0 rho1 rho2 rmax12 rmax13 rmax23 <4xflags>

Exponentially decaying form is also available in Vessal form: E(three) = 1/4*A*(B**2).exp(-r12/rho1).exp(-r13/rho2) A = k/(2*(theta0-pi)**2) B = (theta0-pi)**2 - (theta-pi)**2 Format: three vessal atom1 atom2 atom3 k theta0 rho1 rho2 rmax12 rmax13 rmax23 <4xflags>

k3 and k4 terms can also be included: E(three) = 1/2*k*(theta-theta0)**2 + 1/6*k3*(theta-theta0)**3 + 1/24*k4*(theta-theta0)**4 Formats: three k3 atom1 atom2 atom3 k k3 theta0 rmax12 rmax13 rmax23 <3xflags> three k4 atom1 atom2 atom3 k k4 theta0 rmax12 rmax13 rmax23 <3xflags> three k3 k4 atom1 atom2 atom3 k k3 k4 theta0 rmax12 rmax13 rmax23 <4xflags>

Cosine of theta can also be used instead of theta, if the cosine option is used :

E(three) = 1/2*k*(cos(theta)-cos(theta0))**2

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species where 1-2 and 1-3 are bonded.

The type of bond sub-option can be used to control which bonded atoms the potential applies to. For a three-body potential, the bond type for the two bonds connected to the pivot atom can be supplied and then checked before applying the potential. Valid type of bond options are: single, double, triple, quadruple, resonant, amide. In addition, a bond may specified as regular (default), cyclic or exocyclic. the order of the bond type options matches the order of the pivot atom - end atom pairs. For example, if a potential acts for the triad S - C - O, where the C-S bond is a regular bond and the C-O is an exocyclic double bond then the input would look like:

three bond single regular double exocyclic C core S core O core ...... etc

Conditions can also be placed on a potential to check the number of bonds associated with the pivot atom (e.g. if a potential is designed for a 90 degree angle it might be necessary to check that the number of bonds is four or six for square planar or octahedral). The conditions are placed using one or more instances of "nbeq" or "nbne" followed by an integer. Here "nbeq" imples no. of bonds must equal and "nbne" number of bonds must not equal. Hence, "nbeq 4" would require the pivot to have 4 bonds. When multiple terms are used, the logic of nbeq is concatinated with "or", whereas that of nbne is joined by "and". Therefore "nbne 4 nbne 6" would apply to an atom with 3 bonds, but not one with 4 or 6.

NB It's important to specify the first bond as being of "regular" type so that the "exocyclic" attribute is correctly assigned to the second bond since the first term is assumed to apply to the first bond, whereas the second applies to the second bond.

E(three) = K * (1 + isign*cos(n*theta))

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species were 1-2, 2-3 and 1-3 are bonded.

E(three) = K * (C0 + C1*cos(theta) + C2*cos(2*theta))

Here the coefficients are related to theta0 by:

C2 = 1/(2*sin(theta0))**2 C1 = - 4*C2*cos(theta0) C0 = C2*(2*(cos(theta0))**2 + 1)

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species were 1-2, 2-3 and 1-3 are bonded.

E(three) = (2K/n**2) * (1 - cos(n*theta)) + 2K*exp(-beta*(r13 - r0))

Here r13 is the distance between the atoms 1-3 in the triad.

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species were 1-2, 2-3 and 1-3 are bonded.

E(three) = k (1+3*cos(theta1)*cos(theta2)*cos(theta3)) ------------------------------------------- (r12*r13*r23)**3

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species were 1-2, 2-3 and 1-3 are bonded.

E(three) = K.P(Q1,Q2,Q3).exp(-rho.Q1)

P(Q1,Q2,Q3) = c0 + c1.Q1 + c2.Q1**2 + c3(Q2**2 + Q3**2) + c4.Q1**3 + c5.Q1(Q2**2 + Q3**2) + c6(Q3**3 - 3.Q3.Q2**2) + c7.Q1**4 + c8.Q1**2.(Q2**2+Q3**2) + c9.(Q2**2+Q3**2)**2 + c10.Q1(Q3**3 - 3.Q3.Q2**2)

R1 = (r12-r012)/r012 R2 = (r13-r013)/r013 R3 = (r23-r023)/r023

Q1 = (R1+R2+R3)/sqrt(3) Q2 = (R2-R3)/sqrt(2) Q3 = (2*R1-R2-R3)/sqrt(6)

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species were 1-2, 2-3 and 1-3 are bonded.

E(three) = K * exp(-rho1*r12).exp(-rho2*r13).exp(-rho3*r23)

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species where 1-2, 2-3 and 1-3 are bonded.

E = A.exp(rho/(r-rmax)).(B/r**4 - 1)

E = A.exp(rho/(r-rmax)).(B/r**4 - 1)*F(qi,qj)

F(qi,qj) = exp(1/Q)exp(1/(qi + qj - Q)) if qi+qj < Q F(qi,qj) = 0 if qi+qj > Q

See Chem. Eng. Sci., 49, 2991 (2000) for more details.

E(three) = K * exp(rho1/(r12-rmax12) + rho2/(r13-rmax13))(cos-ct0)**2

where cos = cos(theta(jik)) and ct0 = cos(theta0)

intra/inter can also be specified for molecular calculations, as can bond. Note that unlike most potentials, the specification of the cut-off distances is still required when using the "bond" sub-option since the values act as parameters to the model.

E(three) = K*F12*F13*exp(rho1/(r12-rmax12) + rho2/(r13-rmax13))(cos-ct0)**2

where cos = cos(theta(jik)) and ct0 = cos(theta0)

F12(q1,q2) = exp(1/Q)exp(1/(q1 + q2 - Q)) if q1+q2 < Q F12(q1,q2) = 0 if q1+q2 > Q F13(q1,q3) = exp(1/Q)exp(1/(q1 + q3 - Q)) if q1+q3 < Q F13(q1,q3) = 0 if q1+q3 > Q

See Chem. Eng. Sci., 49, 2991 (2000) for more details.

intra/inter can also be specified for molecular calculations, as can bond. Note that unlike most potentials, the specification of the cut-off distances is still required when using the "bond" sub-option since the values act as parameters to the model.

E(three) = K * (1 + b*cos(n.theta)**m) * (r12 - r1) * (r13 - r2)

intra/inter can also be specified for molecular calculations, as can bond. If bond is given no cutoff is required as the potential will only act between species were 1-2 and 1-3 are bonded.