JDFTx  1.7.0
Setting up molecule geometries

The previous tutorial showed a water molecule calculation in a cubic box of side 10 bohrs. JDFTx is a plane-wave DFT code, which means that it naturally incorporates periodic boundary conditions. However, for molecules we usually want to calculate the properties of one isolated molecule without incurring interactions between the molecule and its periodic images. JDFTx can eliminate image interactions in general geometries, for molecules, nanowires, slabs etc. using truncated Coulomb potentials. (See command coulomb-interaction.) This tutorial shows you the most convenient and efficient way to set up molecule calculations, starting from molecule geometries in the XYZ format, often used by the chemistry community.

Let's stick to water, but this time fetch the geometry from the NIST CCCBDB database. Picking one of the many calculated geometries of water from that database, we can create water.xyz containing:

3

O    0  0        0.127118
H    0  0.75799 -0.50847
H    0 -0.75799 -0.50847

The first line is the number of atoms, the second line is a comment line that is ignored, and the subsequent lines specify the atom type and cartesian coordinates in Angstroms.

Now, we can use the script xyzToIonposOpt from the JDFTx scripts directory to convert this to JDFTx input (see Getting started to add this directory to your run PATH, or specify full path to script below):

xyzToIonposOpt water.xyz 15

This script outputs ionic positions in water.ionpos and lattice vectors in water.lattice after centering the molecule at the origin and rotating it such that you can use the smallest box with a periodic image distance of 15 bohrs. Using these outputs, we can create the input file water.in:

include water.lattice
include water.ionpos
coords-type cartesian

ion-species GBRV/$ID_pbe.uspp        #Replace $ID by any atom symbol to find a pseudopotential
elec-cutoff 20 100                   #Recommended cutoffs for GBRV pseudopotentials
dump End None                        #Don't output anything for this example

which we can run as

jdftx -i water.in | tee water.out

In this input file, we specified the entire set of GBRV pseudopotentials using the wildcard syntax of command ion-species instead of selecting O and H pseudopotentials specifically. Therefore, we could replace water with another molecule containing other atom types with almost no change to the input file, making it extremely easy to set up similar calculations for a large number of molecules. Now, check the final energy of the water molecule (in Hartrees) using the script:

listEnergy water.out

Next, check how the energy changes upon increasing the box size. Replace 15 with 20, 25 and then 30 in xyzToIonposOpt above, and rerun jdftx and listEnergy for each. You could do this manually, but could easily script all these calculations in the bash shell and retrieve their energies using:

#!/bin/bash
for padding in 15 20 25 30; do
    xyzToIonposOpt water.xyz $padding
    jdftx -i water.in | tee water-$padding.out
done
listEnergy water-*.out

Note that you could directly type the above in the shell, or save them to a file, say run.sh, give execute permissions using chmod +x run.sh and then execute using ./run.sh

Finally, also add the coulomb trunaction commands to the above water.in:

coulomb-interaction Isolated         #Switch from periodic to isolated images mode
coulomb-truncation-embed 0 0 0       #Center of molecule; origin for breaking periodicity

and once again run the calculation for the four box sizes.

I get the following energies with and without Coulomb truncation for the different box sizes:

Padding Without truncation With truncation
15 -17.268078941702573 -17.267860822043808
20 -17.268005076655513 -17.267899023265009
25 -17.267959775753003 -17.267901318302076
30 -17.267937318025449 -17.267901948787543

Note that the energy converges much faster with box size after including truncation. In fact, with truncation, the energy converges exponentially with unit cell size, compared to an inverse-cubic convergence without truncation due to dipole-dipole interactions.

In this example, the energy with 15 bohr spacing with Coulomb truncation is within 5x10-5 Hartrees of the value converged with box size, whereas a 30 bohr spacing is required to achieve that accuracy without truncation. Note that the calculation with the 30 bohr spacing is roughly 8 times slower than that with the 15 bohr spacing, so using optimized unit cells with truncated Coulomb interactions can substantially speed up your calculation.