The current status and capabilities of the atomistic simulation code GULP are described.  In particular, the differences between versions 1.3.2 and 3.0 are detailed, as well as a concise pointer to applications in computational crystallography

Gale J.

Rohl A.

Woodley S.

English

Crossref

GULP: Capabilities and prospects

Gale, Julian

espace@Curtin

GULP: Capabilities and prospectsJulian D. Gale*Nanochemistry Research Institute, Department of Applied Chemistry, Curtin University of Technology, P.O. Box U1987, Perth 6845,Western AustraliaReceived July 9, 2004; accepted October 3, 2004Atomistic simulation / Lattice dynamics / Phonons /Electrostatics / GULP Computer program /Computational crystallographyAbstract. The current status and capabilities of the ato-mistic simulation code GULP are described. In particular,the differences between versions 1.3.2 and 3.0 are de-tailed, as well as a concise pointer to applications in com-putational crystallography.1. IntroductionThe program GULP [1–3] is designed as a general forcefieldsimulation package with a strong emphasis on condensedphases. In the present description of the code the functional-ity that will be available in the next release, 3.0, is described.However, the core capabilities for 3-D solids and finite clus-ters are already included in the present release, 1.3.2.While forcefield molecular dynamics codes are numer-ous, instead GULP is designed primarily for the field oflattice dynamics. In short, the program aims to providemost properties that can readily be determined based onthe energy and analytic first, second and third derivativesof this quantity with respect to atomic coordinates andlattice strains. Furthermore, the code is designed to handlesystems of 0-D (molecules, clusters and embedded de-fects), 1-D (polymers and dislocations), 2-D (slabs, sur-faces and grain boundaries), and 3-D solids. In the latercase, the use of crystallographic symmetry is availablesuch that structures can be built using the space group.This also has the benefit of accelerating any structural op-timisation through reducing the number of variables andonly evaluating energy terms for the asymmetric unit.2. CapabilitiesMost calculations in GULP consist of the optimisation ofa trial structure to the local energy minimum, under givenconditions of pressure and temperature, followed by thecalculation of the physical properties that depend on thecurvature about that point. However, there has been con-siderable work by Woodley et al. [4] to extend this to glo-bal optimization through the use of genetic algorithms. InTable 1 a concise list of the key properties that can becalculated with GULP is presented, as well as the genericcapabilities.In order to calculate the energy and its derivatives inGULP it is necessary to define the functional form of theinteractions between particles and any parameters that arerequired. Normally the energy can be regarded as being ex-panded as a sum over terms that depend on increasing num-bers of atomic centres. Within the code at present two-, three-and four-body interatomic potentials are allowed for, as wellas certain specific many-body functional forms. The typesof forcefield implemented span those typically required tosimulate most classes of material, including ionic solids,semiconductors, metals and organic molecular crystals.Aside from the interatomic potentials, a key contribu-tion to the energy for many systems is the electrostaticenergy. GULP has the capability to perform convergentelectrostatic lattice sums for all dimensionalities, includingthe Ewald sum for 3-D [5] and the Parry sum for 2-D [6,7]. The charges can either be fixed via the input to eitherformal or partial values, or can be determined as a func-tion of structure according to electronegativity equalisa-tion. Coulomb subtraction and scaling within molecules isalso available. Furthermore, dipolar ion polarization canbe included via the shell model [8], and radial compressi-bility via the breathing shell model [9].A further capability that can be useful in a crystallo-graphic context is the ability to incorporate partial occu-pancies for sites. This is referred to as either the mean-field or virtual crystal approximation, in which the interac-tions between species are weighted according to their siteoccupancies. Where different chemical species exist withpartial occupancies on the same site, they are automaticallyconstrained to move as a single entity within the program.Such an approach is useful for simulating crystallographicaverages for disordered materials, though thermodynamicinformation, including excess quantities, typically requiresexplicit supercell orderings for an accurate description.3. AvailabilityGULP is available to academics free of charge subject to theconditions that the program is not further distributed with-out permission or used for commercial research. More infor-552 Z. Kristallogr. 220 (2005) 552–554# by Oldenbourg Wissenschaftsverlag, München* e-mail: julian@power.curtin.edu.aumation can found at the web site http://gulp.curtin.edu.au/.While version 1.3.2, and those prior to this, were writtenlargely in f77, with non-standard extensions to enablesome limited use of dynamic memory for arrays that areproportional to the square of the size of the problem, ver-sion 3.0 is written in f90 and has fully dynamic memoryallocation. Support is provided for Unix/Linux-based oper-ating systems, including Mac OSX. For problems that re-quire only first derivatives, such as molecular dynamics orconjugate gradient optimization, then parallel executioncan be performed using MPI. At present this is based onreplicated data, with data distribution over processorseither via atom number or by domain decomposition.4. ApplicationsGULP has been widely applied in many aspects of solidstate chemistry, physics, mineralogy and crystallography.Both the functionality of the code and the literature con-cerning its application have recently been reviewed [3].Hence, we will present just a brief illustration here.While there are many programs based on different Ha-miltonians that make it possible to optimize the structureof a material, the capability of GULP that is more rarelyfound is the ability to minimize under symmetry con-straints and to rapidly analyze the phonon behaviour ofthe final structure. It is often crucial to determine the pho-non modes at the Brillouin zone centre, especially follow-ing a symmetry constrained optimization, in order to ver-ify that the system is stable within the given space group.In order for a system to be at an energy minimum it isnecessary that the first three phonon modes (i.e. theacoustic branch) be equal to zero, as a consequence oftransitional invariance. Furthermore, the remaining vibra-tional frequencies must all be real, instead of imaginary.When this is not the case, this indicates that the imposedsymmetry is too high and that the space group must bechanged. This information can only be obtained from thephonon spectrum and not from the Hessian matrix, whichwill show no instabilities under the space group con-straint.During studies of silica and aluminophosphate poly-morphs [10, 11], many of which are microporous phases,there have been a number of instances of theoretical pre-diction of symmetry lowering based on lattice dynamicalinformation. One of the best known examples is that of thematerial silicalite, which exhibits an orthorhombic structureat elevated temperatures, such as those employed when itis used as a heterogeneous catalyst. However, energy mini-mization followed by a phonon calculation indicates that atlower temperatures, as implied by the neglect of thermaleffects in the simulation, the system should distort to elevi-ate the imaginary mode that is found when studying theorthorhombic unit cell. Indeed, it is now well-known thatsilicalite undergoes a monoclinic-orthorhombic transitionwith temperature. However, the calculations demonstratethat the lower symmetry is an intrinsic property of the frame-work structure, rather than being a consequence of the pre-sence of impurities in the material. Although we shouldpoint out that the original calculations on silicalite [12]were not actually performed with GULP, the result is re-produced with the present code. More recently, the abilityto determine analytical derivatives of vibrational free energywithin GULP [13] has been used to examine the tempera-ture dependent variation of the monoclinic distortion in si-licalite. Although the quasiharmonic approximation breaksdown before phase transition is reached, the results ofGrau-Crespo et al. [14] indicate that the monoclinic angledoes indeed tend towards a right angle as the temperaturerises. It should also be noted that the use of a shell-modelis crucial for the proper description of such distortionssince they are driven by polarization.The existance of imaginary modes in the phonon spec-trum is not confined to the zone centre. Occasionally, themaximum in the modulus of the imaginary mode may beat some other special point in the Brillouin zone. Here theinverse of the fractional coordinates in reciprocal spaceindicates the dimension of the supercell that the systemwishes to adopt. This situation occured in a recent studyof the (10–14) surface of the mineral calcite (CaCO3)[15]. Examination of the phonons for the relaxed surfaceshowed the presence of an imaginary mode at (0.5, 0.0) inthe Brillouin zone. After the creation of a (2 1) super-GULP: Capabilities and prospects 553Capability Version Property VersionOptimization 1.3.2 Elastic constants 1.3.2Free energy minimization 1.3.2 Bulk moduli 1.3.2Molecular dynamics 1.3.2 Shear moduli 3.0Monte Carlo 3.0 Young’s moduli 1.3.2Defect calculations 1.3.2 Piezoelectric constants 1.3.2Forcefield fitting 1.3.2 Dielectric constants 1.3.2Transition state location 1.3.2 Refractive indices 1.3.2Genetic algorithms 1.3.2 Phonon frequencies 1.3.2Phonon free energy 1.3.2 Born effective charges 3.0Non-analytic frequency correction at G 3.0Frequency dependent properties 3.0Electrostatic site potentials 1.3.2Electric field gradient tensor 1.3.2Surface and attachment energy 3.0Table 1. Summary of the main capabilities ofGULP and the properties that can be deter-mined directly by the program. In each casethe applicable version is indicated, with allfunctionality of version 1.3.2 being present inversion 3.0.cell, followed by perturbation along the direction of theimaginary mode eigenvector and re-optimization, the finalsurface structure was determined. Alternate carbonategroups were found to be rotated differently, leading to asuperlattice that is consist with observations from experi-mental LEED work [16].5. Future directionsWhile the program encompasses a wide variety of force-fields and simulation types, there are still many develop-ments of GULP in progress. In order to broaden the rangeof applicability of the code, and to improve the accuracy,new model Hamiltonians are always being added. Presently,the many different forms of bond order potentials are beingintroduced, with those due to Tersoff [17] and Brenner et al.[18] already implemented. These are particularly useful formany covalent systems of current interest, such as siliconand carbon nanotubes. The trend towards including modelsthat are increasingly influenced by quantum mechanical ormany-body effects will also continue. Examples, of some ofthe development work includes the incorporation of ligandfield effects into the description of transition metal systems[19], especially those that exhibit a Jahn-Teller distortion,and use of periodic MNDO-based Hamiltonians to achievea self-consistent dipolar and quadrupolar description ofatomic polarization in the solid state [20].Acknowledgments. We would like to thank the Government of Wes-tern Australia for its support through the Premier’s Research Fellow-ship programme, as well as funding from Curtin University of Tech-nology. The contribution of the ARC to recent developments throughthe provision of an IREX grant is also appreciated.References[1] Gale, J.: GULP: A computer program for the symmetry-adaptedsimulation of solids. J. Chem. Soc., Faraday Trans. 1 (1997)629–637.[2] Gale, J.: Empirical potential derivation for ionic materials. Phil.Mag. B73 (1996) 3–19.[3] Gale, J.; Rohl, A.: The general utility lattice program. Mol. Si-mul. 29 (2003) 291–341.[4] Woodley, S.; Battle, P.; Gale, J.; Catlow, C.: The prediction ofinorganic crystal structures using a genetic algorithm and energyminimisation. Phys. Chem. Chem. Phys. 1 (1999) 2535–2542.[5] Ewald, P.: Die Berechnung optischer und elektrostatischer Git-terpotentiale. Ann. Phys. 64 (1921) 253–287.[6] Parry, D.: The electrostatic potential in the surface region of anionic crystal. Surf. Sci. 49 (1975) 433–440.[7] Parry, D.: Errata: The electrostatic potential in the surface regionof an ionic crystal. Surf. Sci. 54 (1976) 195–195.[8] Dick, B.; Overhauser, A.: Theory of the dielectric constants ofalkali halide crystals. Phys. Rev. 112 (1958) 90–103.[9] Schroeder, U.: A new model for lattice dynamics (breathingshell model). Solid State Commun. 4 (1966) 347–349.[10] Henson, N.; Cheetham, A.; Gale, J.: Computational studies ofaluminium phosphate polymorphs. Chem. Mater. 8 (1996) 664–670.[11] Henson, N.; Cheetham, A.; Gale, J.: Theoretical calculations onsilica frameworks and their correlation with experiment. Chem.Mater. 6 (1994) 1647–1650.[12] Bell, R.; Jackson, R.; Catlow, C.: Computer simulation of themonoclinic distortion in silicalite. JCS Chem. Comm. (1990)782–783.[13] Gale, J.: Analytical free energy minimization of silica poly-morphs. J. Phys. Chem. B102 (1998) 5423–5431.[14] Grau-Crespo, R.; Acuay, E.; Ruiz-Salvador, A.: A free energyminimization study of the monoclinic-orthorhombic transition inmfi zeolite. JCS Chem. Commun. (2002) 2544–2545.[15] Rohl, A.; Wright, K.; Gale, J.: Evidence from surface phononsfor the (2 1) reconstruction of the (10–14) surface of calcitefrom computer simulation. Am. Miner. 88 (2003) 921–925.[16] Stipp, S.: Toward a conceptual model of the calcite surface: Hy-dration, hydrolysis, and surface potential. Geochimica Cosmo-chimica Acta 63 (1999) 3121–3131.[17] Tersoff, J.: New empirical-model for the structural-properties ofsilicon. Phys. Rev. Lett. 56 (1986) 632–635.[18] Brenner, D.; Shenderova, O.; Harrison, J.; Stuart, S.; Ni, B.;Sinnott, S.: A second-generation reactive empirical bond order(rebo) potential energy expression for hydrocarbons. J. Phys.:Condens. Matter 14 (2002) 783– 802.[19] Woodley, S.; Battle, P.; Catlow, C.; Gale, J.: Development of anew interatomic potential for the modeling of ligand field ef-fects. J. Phys. Chem. B105 (2001) 6824–6830.[20] Gale, J.: Semi-empirical methods as a tool in solid state chemis-try. Faraday Discuss. 106 (1997) 219–232.554 Julian D. Gale

https://espace.curtin.edu.au/bitstream/20.500.11937/35300/2/19420_downloaded_stream_512.pdf

GULP: Capabilities and prospects

Abstract

Similar works

Full text

Available Versions

Crossref

espace@Curtin