xcontrol instruction set is inspired by the Turbomole
file syntax. I decided to call it
xcontrol instructions back than,
but here we will just call it (detailed) input for convenience.
The man page can be found here.
The parser implemented is more general and limited by arbitrary choice to this syntax. At some point more common formats like JSON, YAML or XML might become available as alternative input formats.
To read an input file called
> xtb --input xtb.inp coord
In the detailed input you have control about almost every global variable in the program, some instructions even check your input, but most of the time you should know what you are doing. Developed as a feature for developers, this is incredible powerful and naturally way to complicated for the average application. So in most cases you can safely rely on the internal defaults or the shipped global configuration file (should usually be the same).
I will walk you through some selected instructions you might find useful for your application.
Fixing, Constraining and Confining¶
xtb different concepts of constraints are implemented,
so you should know which tool is best for you problem before you
start writing the detailed input.
We will go through this sections using the caffeine molecule
24 C 1.07317 0.04885 -0.07573 N 2.51365 0.01256 -0.07580 C 3.35199 1.09592 -0.07533 N 4.61898 0.73028 -0.07549 C 4.57907 -0.63144 -0.07531 C 3.30131 -1.10256 -0.07524 C 2.98068 -2.48687 -0.07377 O 1.82530 -2.90038 -0.07577 N 4.11440 -3.30433 -0.06936 C 5.45174 -2.85618 -0.07235 O 6.38934 -3.65965 -0.07232 N 5.66240 -1.47682 -0.07487 C 7.00947 -0.93648 -0.07524 C 3.92063 -4.74093 -0.06158 H 0.73398 1.08786 -0.07503 H 0.71239 -0.45698 0.82335 H 0.71240 -0.45580 -0.97549 H 2.99301 2.11762 -0.07478 H 7.76531 -1.72634 -0.07591 H 7.14864 -0.32182 0.81969 H 7.14802 -0.32076 -0.96953 H 2.86501 -5.02316 -0.05833 H 4.40233 -5.15920 0.82837 H 4.40017 -5.16929 -0.94780
In the exact fixing approach the Cartesian position of the selected atom is fixed in space by setting its gradient to zero and the degrees of freedom are removed from the optimization procedure and therefore the atoms stay in place in geometry optimizations.
For dynamics this exact fixing is automatically deactivated, since it usually leads to instabilities in the simulation.
To activate the exact fixing for atoms 1–10 and atom 12 as well as for all oxygen atoms, add
$fix atoms: 1-10,12 elements: O $end
to your detailed input, the atoms keyword refers to the numbering of the individual atoms in your input geometry. With this input the verbose output will show a short summary of the fixed atoms:
------------------------------------------------- | Fixed Atoms | ------------------------------------------------- * 13 fixed atom positions, i.e. in gradient # Z position/Å 1 6 C 1.0731700 0.0488500 -0.0757300 2 7 N 2.5136500 0.0125600 -0.0758000 3 6 C 3.3519900 1.0959200 -0.0753300 4 7 N 4.6189800 0.7302800 -0.0754900 5 6 C 4.5790700 -0.6314400 -0.0753100 6 6 C 3.3013100 -1.1025600 -0.0752400 7 6 C 2.9806800 -2.4868700 -0.0737700 8 8 O 1.8253000 -2.9003800 -0.0757700 9 7 N 4.1144000 -3.3043300 -0.0693600 10 6 C 5.4517400 -2.8561800 -0.0723500 12 7 N 5.6624000 -1.4768200 -0.0748700 8 8 O 1.8253000 -2.9003800 -0.0757700 11 8 O 6.3893400 -3.6596500 -0.0723200
Since version 6.3 the input is sorted and duplicates are removed automatically.
Almost absolute control about anything in your system is archived
by applying constraining potentials. First of all the constraining
potentials offer a weaker version of the exact fixing, which is
invoked by the same syntax in the
$constrain data group as
$constrain atoms: 11 elements: C,N,8 $end
the program will not attempt to hold the Cartesian positions constant, but the distances between all selected atoms, here number 11 and all carbon, nitrogen and oxygen. For each atom pair a harmonic potential is generated to hold the distances at roughly the starting value, this even works without problems in dynamics.
For your caffeine molecule this results in a problem, which can easily be spotted in the verbose output of the constraints summary.
------------------------------------------------- | Constraints | ------------------------------------------------- * 15 constrained atom positions positions referring to input geometry # Z position/Å displ./Å 11 8 O 6.3577347 -3.6327225 -0.0684681 0.0000000 1 6 C 1.0687445 0.0520162 -0.0755782 0.0000000 3 6 C 3.3535252 1.0744217 -0.0774722 0.0000000 5 6 C 4.5969189 -0.6303196 -0.0745855 0.0000000 6 6 C 3.2896462 -1.0950551 -0.0735452 0.0000000 7 6 C 2.9629004 -2.4886091 -0.0702591 0.0000000 10 6 C 5.4425717 -2.8389078 -0.0699228 0.0000000 13 6 C 7.0086271 -0.9538835 -0.0749134 0.0000000 14 6 C 3.9536622 -4.7147069 -0.0641970 0.0000000 2 7 N 2.5030143 0.0336686 -0.0754494 0.0000000 4 7 N 4.6213728 0.7205067 -0.0770658 0.0000000 9 7 N 4.1215924 -3.2704219 -0.0683024 0.0000000 12 7 N 5.6601563 -1.4769082 -0.0732293 0.0000000 8 8 O 1.8493654 -2.9780046 -0.0690630 0.0000000 11 8 O 6.3577347 -3.6327225 -0.0684681 0.0000000 applying 105 atom pairwise harmonic potentials applied force constant per pair: 0.0035714 effective force constant per atom: 0.0500000 constraining energy/grad norm: 0.0000000 0.0000000
Since version 6.3 the input is sorted and duplicates are removed automatically.
Note that in some versions of
xtb this leads to
NaN for the
gradient, therefore double-check the constrained atom list for duplicates.
To constrain the atoms more tightly the force constant can be adjusted
$constrain force constant=1.0 $end
this variable goes directly into the constraining procedure and is given in atomic unites (Hartree/Bohr²), for very high force constants this becomes equivalent to the exact fixing. Note the difference in the syntax as you are required to use an equal-sign instead of a colon, as you are modifying a global variable.
It is also possible to constrain selected internal coordinates, possible are distances, angles and dihedral angles as done here
$constrain distance: 1, 2, 1.4 angle: 5, 7, 8, auto dihedral: 3, 4, 1, 7, 180 $end
This printout is not yet fully integrated in the released versions and
might not work as expected.
For the experimental constraining potentials supporting this features
use the all-caps keyword variants (
Note that those are not available for the scan feature.
Distance constraints are given in Ångström, while angle constraints are given
The distances are defined by two atom number referring to the order in
your coordinate input, angles are defined by three atom numbers and
dihedral angles by four atoms, in any case the atoms do not have to
be connected by bonds. The last argument is always the value which should
be used in the constraining potential as reference, if you decide to
use the current value
auto can be passed. The constraints will be
printed to the screen (the newer implementation may require the verbose mode,
to trigger the printout of the constraint summary), we check this setup
for the caffeine molecule
------------------------------------------------- | Constraints | ------------------------------------------------- * 1 constrained distance # Z # Z value/Å actual/Å 1 6 C 2 7 N 1.4000000 1.4409371 constraining potential exponent: 2.0000000 applied force constant per dist.: 0.0500000 effective force constant per atom: 0.0250000 constraining energy/grad norm: 0.0002992 0.0109403 * 1 constrained angle # Z # Z # Z value/° actual/° 5 6 C 7 6 C 8 8 O 150.4357763 150.4357763 applied force constant per angle: 0.0500000 effective force constant per atom: 0.0166667 constraining energy/grad norm: 0.0000000 0.0000000 * 1 constrained dihedral angle # Z # Z # Z # Z value/° actual/° 3 6 C 4 7 N 1 6 C 7 6 C 180.0000000 -179.9396548 applied force constant per angle: 0.0500000 effective force constant per atom: 0.0125000 constraining energy/grad norm: 0.0000000 0.0000629 total constraint energy/grad norm: 0.0002993 0.0110032
You can find the constraint energy and gradient at the end of the summary, check for unphysical high values of the energy and gradient here to verify your constraint setup otherwise you might encounter strange behaviour in the following optimization or dynamics to adhere this constraints.
If you are not quite sure which distances or angles you want to constrain, run
> cat geosum.inp $write distances=true angles=true torsions=true $end > xtb --define --verbose --input geosum.inp caffeine.xyz
and have a look at the geometry summary for your molecule. The
data group toggles the printout in the property section and also some
printouts in the input section.
For caffeine the geometry summary including only the distances looks like this
------------------------------------------------- | Geometry Summary | ------------------------------------------------- molecular mass/u : 194.1926000 center of mass at/Å : 4.0569420 -1.6298957 -0.0733327 moments of inertia/u·Å² : 4.7317175E+02 7.1109348E+02 1.1745947E+03 rotational constants/cm⁻¹ : 3.5626878E-02 2.3706633E-02 1.4351872E-02 * 25 selected distances # Z # Z value/Å 1 6 C 2 7 N 1.4409371 2 7 N 3 6 C 1.3698478 3 6 C 4 7 N 1.3186949 4 7 N 5 6 C 1.3623047 2 7 N 6 6 C 1.3652477 5 6 C 6 6 C 1.3618461 6 6 C 7 6 C 1.4209574 7 6 C 8 8 O 1.2271501 7 6 C 9 7 N 1.3977057 9 7 N 10 6 C 1.4104346 10 6 C 11 8 O 1.2347703 5 6 C 12 7 N 1.3741439 10 6 C 12 7 N 1.3953559 12 7 N 13 6 C 1.4514011 (max) 9 7 N 14 6 C 1.4496299 1 6 C 15 1 H 1.0929740 1 6 C 16 1 H 1.0928728 1 6 C 17 1 H 1.0928837 3 6 C 18 1 H 1.0829302 (min) 13 6 C 19 1 H 1.0932399 13 6 C 20 1 H 1.0945661 13 6 C 21 1 H 1.0945601 14 6 C 22 1 H 1.0927021 14 6 C 23 1 H 1.0949866 14 6 C 24 1 H 1.0949141 * 4 distinct bonds (by element types) Z Z # av. dist./Å max./Å min./Å 1 H 6 C 10 1.0926630 1.0949866 1.0829302 6 C 6 C 2 1.3914017 1.4209574 1.3618461 6 C 7 N 11 1.3941548 1.4514011 1.3186949 6 C 8 O 2 1.2309602 1.2347703 1.2271501
There is no electronic structure information used at this point but a simple geometric model to select distances, which can get too few or too many bonds or angles in this printout.
Confining in a Cavity¶
If you are running dynamics for systems that are non-covalently bound, you may encounter dissociation in the dynamics. If you want to study the bound complex, you can try to confine the simulation in a little sphere, which keeps the molecules from escaping. The detailed input looks like
$wall potential=logfermi sphere: auto, all $end
You can be more precise on the radius by giving the value in Bohr instead
auto. The automatically determined radius is based on the largest
interatomic distance in the structure plus some offset.
The logfermi potential is best suited for confinements, but not yet the
When using this input with the caffeine molecule the automatically determined radius is about 5.6 Å, which should be large enough to contain a molecule of its size. At first it might be surprising to find that the confining energy is about +84 kcal/mol, but we did we did not account for the placement of the molecule, relative to the center of potential within our chosen input. Currently, the center point of the spherical logfermi potential is set at the origin (0,0,0) and the center of mass of the caffeine molecule is about 4.4 Å away from it, so our molecule is stuck halfway in the wall we just created.
To cope with this we should put the center of mass of the caffeine molecule
at the origin, this can be done by adding the
$cma instruction to the input
file, which shifts the coordinates with the center of mass and aligns
the molecule to its principal axes of inertia.
The center point for all wall potentials is always placed at the origin, which cannot be changed with the currently available input options. Therefore, we resort to modifying our input coordinates here.
Different Potential Shapes¶
Currently two different potential shapes are implemented and can be
selected with the
The logfermi potential shape is given by the expression
kB is the Boltzmann constant,
T is formally the temperature but can be used to scale
the strength of the potential (adjustable by
temp=<real>, within the
β is the steepness of the potential (adjustable by
RA are the cartesian coordinates of atom A,
O is the origin (0,0,0) and
Rsphere is the radius of the sphere used for confining.
The default potential shape is a simple polynomial to the power α
The formula that is evaluated in the program is
The main (dis)advantage of this shape is that the radius of the sphere
is a relative quantity compared to the size of the molecule.
auto generator of the sphere radius takes this into account,
by rescaling the largest distance in the molecule instead of adding
a constant shift.
A clear disadvantage of this potential shape it that the gradient does
not vanish inside the sphere and can compress a molecule artificially.
For some molecules an isotropic spherical cavity is not suitable for confinement, since the molecule might have a rod-like or oblate shape. Instead of sphere we can use an ellipsoid to construct an anisotropic cavity, there is no limitation for the potential shape since we use a simple rescaling to introduce anisotropy.
The input file for an anisotropic potential would look like
$wall potential=logfermi ellipsoid: 13.5,11.1,8.6,all # values in Bohr $end
As for the isotropic one can use the
auto keyword to replace any of the
the three radii with an automatically determined value.
The automatic determined value is the automatic isotropic sphere radius,
so letting all three values be autodetermined results in an isotropic potential.
As before, we have to deal with the issue that the center of mass of our caffeine molecule and the origin do not coincident, this time we use a Python interpreter with ASE support for this job
from ase.io import read, write mol = read('caffeine.xyz') mol.set_positions(mol.get_positions() - mol.get_center_of_mass()) write('caffeine_shifted.xyz', mol)
Finally we can check
xtb with the new coordinates and the above input file and
we find that the confining energy is zero in the initial geometry.
Using Multiple Potentials¶
Since version 6.0 an arbitrary number of wall potentials is supported.
Similar to the constraint keywords one could create multiple wall potentials
ellipsoid instructions like
$wall potential=logfermi sphere: auto, all ellipsoid: 13.5,11.1,8.6,all # values in Bohr $end
This could be used to confine different fragments in different sized spheres. The only restriction is that the potential shape is global.
As I promised you can control almost everything, the
xcontrol(7) man page
is a good starting point to get acquainted with the detailed input. This
poses the usual hindrance of actually reading the documentation
(since you are here, you are already above average, thumbs up).
A practical alternative is to just dump the complete internal settings of the program to an input file and start playing around with it. To do so, run
> xtb --input default.inp --define --copy coord
default.inp has not to be present when starting the program
--copy mode, since the
default.inp will be generated for you.
--define flags makes sure that the program only checks your setup
and does not perform any calculation on the input coordinates.
Have a look at the first lines of
$cmd xtb --input default.inp --define --copy coord $date 2019/03/05 at 08:50:26.651 $chrg 0 $spin 0 ...
This is actually the command you used in the first place to invoke the
program, next you find the timestamp when the program was started and
then system specific information about charge and spinstate of your system,
this is what I understand as a self-documenting program run.
$date are cosmetic features and will never influence
any calculation if included in the detailed input, but I figured that
they might become handy if you look back into your calculations when
putting together the manuscript or taking over a project from your,
now graduated, fellow coworker.
The rest of the file represent every accessible variable documented
xcontrol(7) man page with its current setting, this should be
quite a lot. So lets focus say on the
... $wall potential=polynomial alpha=30 beta=6.000000000000000 temp=300.0000000000000 autoscale=1.000000000000000 axisshift=3.500000000000000 ...
The default potential is a
polynomial one, you want to change this to
alpha is only needed for the
potential, we use
temp in our potential.
The steepness of our potential can be adjusted by modifying the value
beta, since our potential is multiplied with the thermal energy
we can scale it by increasing it temperature in
autoscale is a factor the automatic determined sphere axes are
multiplied with, a default of 1.0 seems reasonable here, but sometimes
we need more space or want to squeeze everything a bit together.
We can also adjust the constant shift value used in the generation
of the automatic axes, but on a second thought this value might be
just fine, so we do not modify
This is an awful lot of information in a small block and quite essential
for your calculation using a confining potential, all details on this
can be found in
xcontrol(7) man page at the group instruction
If you are happy with all this setting you can just use this file as
.xtbrc and place it somewhere in your
Global Configuration File¶
The global configuration file called
.xtbrc has to be around somewhere
xtb is able to find it and uses the very same
syntax as the detailed input. Every instruction (
key=value) you can
use in your detailed input file can be present in your global configuration
file. System specific instructions (
key: value) will not work, of course.
To check which
.xtbrc is read, start the program in verbose mode and
check the Calculation Setup section in the output.
Rules for Control¶
This section is intended to briefly explain the currently applied rules to parse the detailed input file.
Every instruction is started by a flag (
$) and terminated by the next flag.
An instruction is only valid if the flag is in the first letter, the
instruction name is the rest of the register (line).
A valid instruction opens its block with its own options, every option
is a key-value pair. Invalid instructions are ignored without further warning.
There are two kind of instructions, logical and groups. Logical instructions toggle a specific operation and cannot contain a option block while group instructions only open the option block without any further actions.
Groups with the same name can occur multiple times and are merged before parsing.
There are two kind of options,
key=value pairs set global variables and
can only be used once, they are locked at the first encounter of the key,
regardless of the value, in case an invalid value is given the default is used
as fallback and cannot be modified by subsequent options with the same key.
Options of the kind
key: values,... can be present multiple times and
are handled differently depending on the context they are used in.
For example the
atoms: instruction usually appends atoms to an list,
$constrain applies a quadratic potential to
the atom pair specified.