GULP

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

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)

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 EEM 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.

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.

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))

Note that the cut-offs are set by the manybody potential

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))

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 however 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

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.

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.

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.

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

Note that in the case of the Voter taper form there is no taper range since it applies over the whole 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>

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>

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.

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.

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.

E = De.((1-exp(-a(r-r0)))**2 - 1.0) -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.

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.

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. 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.

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.

E(three) = K * (isign*cos(n*theta) + 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) = 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.

E(three) = K * (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.

E(three) = [K1*(r12 - r1) + K2*(r13 - r2)].(theta-theta0)

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.

E(three) = 1/2 * K * (r - r0)**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.

title N (where N=no of title lines) <line 1> <line 2> : : <line N>

OR

title <line 1> <line 2> : : <line N> end

E = 1/2 K2*(r-r0)**2 + 1/6 K3*(r-r0)**3 + 1/24 K4*(r-r0)**4 -coul*qi*qj/r

k3 and k4 terms can be included as follows:

harm k3 atom1 atom2 k2 k3 r0 <coul> <rmin> rmax <3*flags>

harm k4 atom1 atom2 k2 k4 r0 <coul> <rmin> rmax <3*flags>

harm k3 k4 atom1 atom2 k2 k3 k4 r0 <coul> <rmin> rmax <4*flags>

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 = 1/4 K2*(r**2 - r0**2)**2 - 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.

Example :

bornq 1 1 xx 1.95

Note: Born effective charges cannot currently be calculated with electronegativity equalisation.

Example :

monopoleq 1 1 1.95

E = c0 + c1*r + c2*r**2 + ... + cn*r**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.

E = A.exp(-r/rho)/r**m - C/r**n <-E(rmax) <+Grad correction term>>

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 = -(C6/r**6)*f6(r) -(C8/r**8)*f8(r) -(C10/r**10)*f10(r)

where f2n(r) is given by:

f2n(r) = 1 - { sum(k=0->2n) [(b*r)**k]/k!} * exp(-b*r)

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. 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.

Note that this potential cannot yet be used with the c6 keyword.

For geometric mean fit: constrain <fit> <no. of constraints> nvari1 nvari2 mean nfix <coefficient>

Geometric constraints:

Variable numbers are coded as strains = 1-6, internal variable n of atom i = 3*(i-1)+n+6, n=1,3 => x,y,z and radius of atom i = 6+3*nasym+i, where nasym = no. of atoms in asymmetric unit. Now nvari and nfix can be entered as: atom1 x/y/z atom2 x/y/z coefficient offset e.g. 3 x 6 z 1.0 0.0 Radii can be entered as: atom1 r atom2 r coefficient offset e.g. 3 r 6 r 1.0 0.2

Additive use of geometric constraints is allowed so that an atom may be fixed at the centroid of two other atoms. When this is the case for 3-D systems, there is an ambiguity due to the presence of periodic replications. GULP will use the two nearest images for the constraints.

Fitting constraints:

Variable numbers refer to the number of the parameter as included in the variables table, which in general is the same as the order of specification in the input. If nvar1=nvar2 for geometric mean constraints then the square root of just the one parameter is taken. Note that the variable numbers for the constraint command in fitting may be changed in the restart file relative to those in the input. This is to take account of the reordering of potentials.

To input multiple constraints the number can be added after the constraint command.

Note: The ability to raise the variable to power is currently only available for fitting constraint, though it shouldn't be needed for geometric variables anyway.

marvin - generates a file suitable for input to Marvin for surface calculations. Index of surface plane needs to be added first, though this may be passed through GULP to the marvin input but using the "marvin" option. If the file extension ".mvn" is not given by the user it will be added by the program. e.g. output marvin alumina produces alumina.mvn If a file for Marvin2 is required, then use the sub-option marvin2 instead. e.g. output marvin2 alumina

thbrel - converts the final bulk structure into a file format suitable for the THBREL suite of programs. Not all features are readily transferable between the programs so no guarantee is made that the input file is perfect and for similar reasons both programs may not always yield the same results unless the user is careful to make sure the files are equivalent. If a phonon run has been requested then the output is modified for THBPHON instead of THBREL.

xtl - this is only applicable to crystal structures with a single structure per file. Produces a .xtl file for input into the BIOSYM software, though use of the archive file is better unless symmetry is present. If the filename input is not already terminated with ".xtl" then this will be added by the program. e.g. output xtl alumina produces alumina.xtl

xr - this will output a modified CSSR file suitable for input into the Oxford Materials graphical interface program Crysalis. The file will have the extension ".xr" added. At the moment this is only applicable to 3D systems. e.g. output xr alumina produces alumina.xr

cssr - this will output a CSSR file suitable for input into the MSI graphical interface Cerius2. The file will have the extension ".cssr" added. At the moment this is only applicable to 3D systems. e.g. output cssr alumina produces alumina.cssr

arc - alternatively known as a ".car" file. This option produces archive files suitable for input into the BIOSYM Insight software and will handle bulk, cluster and defect calculations, all with multiple structures. The username should input a root name, e.g. output arc alumina The program will then produce archive file names with either "_3D", "_defe", or "_0D" appended to distinguish the files resulting from a particular section of the run. If multiple structures are present then an underscore followed by the structure number will be added, all followed by ".arc". e.g. if the above input for alumina contained two bulk structures then the files produced would be "alumina_bulk_1.arc" and "alumina_bulk_2.arc". If the word "movie" is specified then all structures during an optimisation will be included in the arcfile, rather than just the final one, so that the optimisation may be viewed. e.g. output movie arc alumina Note : for MD runs the name of the archive file is set by "mdarchive" to avoid overwriting the optimisation archive file if present. Also in the case of MD the file is only written during the production phase.

xyz - this will output an xyz file suitable for input into programs such as Molden and with slight modification Unichem. The file will have the extension ".xyz" At the moment this is only applicable to non-periodic systems. When called with a periodic case just the Cartesian coordinates will be output, without the unit cell. e.g. output xyz cluster produces cluster.xyz If the word "movie" is specified then all structures during an optimisation will be included in the xyz file, rather than just the final one, so that the optimisation may be viewed. e.g. output movie xyz cluster

trajectory - this is a binary file which stores the coordinates and velocities from a molecular dynamics run, as well as some of the system properties. This file is used by the Cerius interface to generate a .trj file for analysis in Cerius. Note that this file is only written during the production phase by default. e.g. output trajectory alumina For this sub-option the further option "ascii" may be specified in order to obtain an ASCII format trajectory file. Note that this is less efficient for disk space usage! e.g. output trajectory ascii alumina Furthermore, the sub-option "equil" may be added to force writing during the equilibration phase of the simulation.

history - this is a text file in the DLPOLY HISTORY file format containing the structure and velocities sampled from a MD run. This file can be used for post-MD analysis using the same programs as for DLPOLY such as "After".

fdf - this is a text file in the Flexible Data Format of Alberto Garcia and Jose M. Soler. This file can then form the basis of an input deck for the program SIESTA.

drv - this is a text file containing the energy and appropriate derivatives calculated by GULP at the last function evaluation. This can be used for QM/MM schemes.

frc - this is a text file containing the energy and appropriate force constants for cores only calculated by GULP during a phonon calculation. This can be used for QM/MM schemes, such as in the program QMPOT. A phonon calculation must be specified otherwise no file will be output.

cif - this is only applicable to crystal structures with a single structure per file. Produces a .cif file for input into a variety of programs. If the filename input is not already terminated with ".cif" then this will be added by the program. Note that at present the files are output in P1 regardless of symmetry for simplicity. e.g. output cif alumina produces alumina.cif

str - this is a file in the CRYSTAL98 format, suitable for reading by DLV. If the filename input is not already terminated with ".str" then this will be added. If multiple structures are specified then the filename will be modified so that a separate file can be written for each structure.

phonon - this option leads to the density of states from a phonon calculation and any dispersion curve points being output to the files <filename>.dens and <filename>.disp, respectively, in a format suitable for input into a curve plotting program.

frequency - this option produces a list of frequencies over all k points for use in an analysis program. A typical application might be in the calculation of the Gruneisen parameter by differencing of the frequencies with respect to the cell volume. The file can either be output in binary format (to maintain precision - the default): e.g. output freq binary <filename> or as ordinary text: e.g. output freq text <filename>

pressure - this option outputs a file containing the pressure tensor during a molecular dynamics run, both instanteously and as an average. File is output as <filename>.pre using units of GPa. osc - this is a text file containing the oscillator strengths for the phonon modes so that the frequency dependent data can be handled in post processing.

or if esff option is specified

atom1 atom2 atom3 atom4 k1 k2 <+/->n rmax(1-2) rmax(2-3) rmax(3-4) rmax(4-1) <2 x flags>

E = k*(1+isign*cos(n*phi-phi0))

ESFF :

E = k1*(sin1**2)*(sin2**2)+isign*k2*(sin1**n)*(sin2**n)*cos(n*phi)

intra/inter can also be specified for molecular calculations, as can bond. If bond is given then no cutoffs are required as the potential will act only when i-j, j-k and k-l are bonded.

if the option esff is specified after torsion then the potential takes the ESFF form. Here the potential becomes dependent on the two angles involved in the torsional angle so that it smoothly goes to zero when the torsional angle becomes indeterminate. Strictly the two coefficients k1 and k2 should be related by :

k1 = K / (sin10**2.sin20**2) k2 = K / (sin10**n.sin20**n)

where sin10 and sin20 are the sines of the equilibrium values of angles 1 and 2 respectively. Here they are just input as coefficients for generality.

or if esff option is specified

atom1 atom2 atom3 atom4 k1 k2 <+/->nphase rtaper rmax(1-2) rmax(2-3) rmax(3-4) rmax(4-1) <2 x flags>

E = k*(1+isign*cos(nphase*phi-phi0))*f(r12)*f(r23)*f(r34)

ESFF :

E = [k1*(sin1**2)*(sin2**2)+isign*k2*(sin1**n)*(sin2**n)*cos(n*phi)] *f(r12)*f(r23)*f(r34)

where f(r) is a cosine tapering function.

or if esff option is specified

atom1 atom2 atom3 atom4 k1 k2 <+/->nphase rho12 rho23 rho34 rmax(1-2) rmax(2-3) rmax(3-4) rmax(4-1) <5 x flags>

E = k*(1+isign*cos(nphase*phi-phi0))*exp(-r12/rho12)*exp(-r23/rho23)*exp(-r34/rho34)

ESFF :

E = [k1*(sin1**2)*(sin2**2)+isign*k2*(sin1**n)*(sin2**n)*cos(n*phi)]* *exp(-r12/rho12)*exp(-r23/rho23)*exp(-r34/rho34)

E = 1/2 K(phi - phi0)**2

E = K.cos(phi).(theta - theta0).(theta' - theta0')

where theta is the angle between 1-2-3 & theta' is between 2-3-4

E = k.d**2 + k4.d**4, where d=distance to the plane