GULP

GULP, the General Utility Lattice Program, is a forcefield code.

Instructions

Make sure that the following variables are set correctly:

  • $GULP_LIB : path to the folder containing the potential files

  • $ASE_GULP_COMMAND="/path/to/gulp < PREFIX.gin > PREFIX.got"

The latter defaults to gulp < PREFIX.gin > PREFIX.got, assuming that the gulp command is in $PATH.

GULP uses several variables to write the input file. The most important ones are:

  • keywords: string of space-separated keywords, e.g.: 'angle compare comp conp'

  • options: list of strings with options, e.g.: ['xtol opt 5.0']

  • library: filename of the potential library file (e.g. 'reaxff.lib'). Recall that the search path for this library is given by $GULP_LIB.

If the potential uses different atom types, one must define an instance of conditions with the atom object to label them and add rules to rename atoms. The only rule suported right now is min_distance_rule which selects atoms of one type close to atoms of another type and renames them. For example:

c = Conditions(atoms)
c.min_distance_rule('O', 'H', ifcloselabel1='O2', ifcloselabel2='H',
                    elselabel1='O1')
calc = GULP(conditions=c)

This will assign, for each H in the system, a corresponding O atom as O2, selecting those O which are the closest to H. The rest of oxygens will be labeled O1.

Finally, if the potential file requires the use of shells, the variable shel must be used:

GULP(..., shel=['O'])

Example

Here is an example of how to use the GULP calculator.

# flake8: noqa
import numpy as np

from ase import Atoms
from ase.calculators.gulp import GULP, Conditions

cluster = Atoms(symbols='O4SiOSiO2SiO2SiO2SiOSiO2SiO3SiO3H8',
                pbc=np.array([False, False, False], dtype=bool),
                cell=np.array(
                    [[0., 0., 0.],
                     [0., 0., 0.],
                        [0., 0., 0.]]),
                positions=np.array(
                    [[-1.444348, -0.43209, -2.054785],
                     [-0.236947, 2.98731, 1.200025],
                        [3.060238, -1.05911, 0.579909],
                        [2.958277, -3.289076, 2.027579],
                        [-0.522747, 0.847624, -2.47521],
                        [-2.830486, -2.7236, -2.020633],
                        [-0.764328, -1.251141, 1.402431],
                        [3.334801, 0.041643, -4.168601],
                        [-1.35204, -2.009562, 0.075892],
                        [-1.454655, -1.985635, -1.554533],
                        [0.85504, 0.298129, -3.159972],
                        [1.75833, 1.256026, 0.690171],
                        [2.376446, -0.239522, -2.881245],
                        [1.806515, -4.484208, -2.686456],
                        [-0.144193, -2.74503, -2.177778],
                        [0.167583, 1.582976, 0.47998],
                        [-1.30716, 1.796853, -3.542121],
                        [1.441364, -3.072993, -1.958788],
                        [-1.694171, -1.558913, 2.704219],
                        [4.417516, 1.263796, 0.563573],
                        [3.066366, 0.49743, 0.071898],
                        [-0.704497, 0.351869, 1.102318],
                        [2.958884, 0.51505, -1.556651],
                        [1.73983, -3.161794, -0.356577],
                        [2.131519, -2.336982, 0.996026],
                        [0.752313, -1.788039, 1.687183],
                        [-0.142347, 1.685301, -1.12086],
                        [2.32407, -1.845905, -2.588202],
                        [-2.571557, -1.937877, 2.604727],
                        [2.556369, -4.551103, -3.2836],
                        [3.032586, 0.591698, -4.896276],
                        [-1.67818, 2.640745, -3.27092],
                        [5.145483, 0.775188, 0.95687],
                        [-2.81059, -3.4492, -2.650319],
                        [2.558023, -3.594544, 2.845928],
                        [0.400993, 3.469148, 1.733289]]))


c = Conditions(cluster)
c.min_distance_rule('O', 'H', ifcloselabel1='O2',
                    ifcloselabel2='H', elselabel1='O1')
calc = GULP(keywords='conp', shel=['O1', 'O2'], conditions=c)

# Single point calculation
cluster.calc = calc
print(cluster.get_potential_energy())

# Optimization using the internal optimizer of GULP
calc.set(keywords='conp opti')
opt = calc.get_optimizer(cluster)
opt.run(fmax=0.05)
print(cluster.get_potential_energy())

The script performs a single-point calculation and then an optimization of the H8Si8O20 double ring structure within ASE using the ffsioh potential. The method get_optimizer() returns an object which works like an ASE optimizer, but actually triggers an optimization using GULP’s internal optimizer. This works on GULP calculators that have the 'opti' keyword.

_GULP: http://gulp.curtin.edu.au/gulp/