Before you start

You need to get Sally or Louise to request an account for you on ib-server. Frank will give you a password.

You will also need to be added to the castep19 usergroup. You are covered by Sally's license, but academic licenses are free anyway.

You will find example .param files for testing convergence and running geometry optimizations, and utilities for extracting convergence test results and energies in /home/lsp/for_Archer2/ You may find it convenient to copy these to your local folder.

Every day:

module add intel/2015.6 castep/intel-2015.6/19.11

You will need to type this in most days you are using CASTEP, as the environment is needed by a lot of the commands you will use.

Before you start

There are some good instructions at https://docs.archer2.ac.uk/quick-start/quickstart-users/

You will need:

• An account on safe. Follow the instructions on https://safe.epcc.ed.ac.uk/
• An account on ARCHER2. Log in to safe and choose "Request login account" under "Login accounts".
• Your project to be included in the current funding cycle. Speak to Louise and Sally about this.

You will also benefit from having the following in your work space:

• cif2cell
  • Copy the whole folder /home/lsp/for_Archer2/cif2cell-1.2.10/ from ib-server to Archer2
  • From within that folder, type

    python setup.py build

  • You might also want to add the path to your .bashrc file, by adding

    PATH="/home/e05/e05/uccalsp/cif2cell-1.2.10/:$PATH"

• Example .param files for testing convergence and running geometry optimizations, and utilities for extracting convergence test results and energies. These are all found in /home/lsp/for_Archer2/ on ib-server

Converting your .cif file

You should start from an optimized .cif file, whenever possible. If you have to start from an experimental crystal structure, consider incremental optimization to make things more efficient.

If you are starting from a .res file, be aware that converting to .cif in Mercury will truncate the cell parameters to the errors, so check that the values in the ZERR line are apprpriate. If you have VERY old DMACRYS output, the errors are written as 0.0, and need to be deleted in order to use Mercury for the structure conversion.

CASTEP uses .cell format as input files. cif2cell should be used to convert .cif format to .cell format. The command is:

cif2cell -p castep STEM.cif

where STEM is the structure name.

If you have any problems with this, it is worth checking that the file is in unix format. Use

file STEM.cif

which should NOT report that it has CRLF deliminators. If it does, run dos2unix on the .cif file.

Testing convergence criteria to be used

Before running a minimization in CASTEP, it is necessary to choose what basis set and what k-point spacing to use. You may also want to change the pseudopotential from the default.

When you are starting out, it is simplest to use the example .param file as described here, and scan through a range of energy cutoffs with a certain k-point spacing. However, with experience (and noting that these two quantities are almost completely independent) you may be able to run just one energy cutoff scan with your chosen k-point spacing and a couple of single point calculations with different k-point spacings. You need to be sure that the energy difference between a pair of polymorphs is stable, so run these calculations for two different structures. If you only have 1 structure, then the convergence will require a higher accuracy.

Basis set

The basis set used by CASTEP is a plane-wave basis set, which can be imagined as being made up of Sine and Cosine waves of increasing frequency, and hence accuracy. The highest accuracy will increase the time taken to perform the calculation, and will not necessarily give any better results, hence it is desirable to test what cutoff energy should be selected for this.

In your folder where you have your .cif file, copy an example .param file with the same file seed. Start with the conv.param file in the folder of examples. Within this, you will see:

• "sedc" is the dispersion correction, which is commented out for the convergence testing
• "task" is a singlepoint energy calculation
• "xc_functions" is the DFT functional, usually PBE for our work
• "cut_off_energy" is 1500 eV.  The "finite_basis_npoints" is 9, and the "finite_basis_spacing" is 100 (finite_basis_corr is 2, but finite_basis_corr set to 0 switches off the scan).  This means that 9 single point energy evaluations will be carried out, with a maximum cutoff energy for the final point of 1500 eV, with the earlier points having the maximum cutoff energy increasing by 100 eV per step.  In practice, the first task will be a single point energy calculation with 700 eV as the cutoff energy for the basis set generation.  The resulting model will be used as the starting point for the next step, which will increase the cutoff energy for the basis set generation to 800 eV.  This will continue until the final (9th) step uses a cutoff energy of 1500 eV.
• "grid_scale" and "fine_grid_scale" tell CASTEP how to set up the k-point space. We have found that 4.00 for the fine-grid_scale is best, but many people will say that this is very high.
• From this point down to "elec_energy_tol" is about how the SCF calculation is run.

k-point spacing

This is set at the end of the .cell file, and is of the format

kpoints_mp_spacing 0.10

This is a k-point spacing of 0.10 A^-1. Just add this single line to your .cell file.

pseudopotential

By default, CASTEP calculates ultrasoft-pseudopotentials on-the-fly for each calculation you run. This takes a matter of seconds (about 20 s for C, H, N, O containing compounds), so we would always calculate them in this way. The pseudopotentials are used to calculate the energy of the core electrons for each of the atoms, but are also used somehow in the energy evaluation of the entire crystal structure. For more notes, please see the separate page on pseudopotentials.

Now you are ready to run CASTEP for the first time.

castep
Seed?: IBPRAC_dfEe1_Ee1
Number of cores (default=2): 48
Version (6.11, 7.0.3, 8.0, 18.1, 18.1-intel, 19.11, 19.11-intel): 19.11-intel
Selected version is 18.1-intel
Jobname: JEKNOC12
Your job 471307 ("JEKNOC12") has been submitted

Our group can use up to 96 cores between us in the all.q.  Keep communicating with Rui and Louise about who may use how many cores.  Always use 19.11-intel as the version.  (There is an old note that Louise was having problems using this version - needs testing again...)

If you find that you need to run on more than 96 cores, you can run on up to 160. Run castep as above, requesting up to 160 cores. This will not run as there are not 160 cores available in the all.q for us. So kill the job and edit the .sh file that was generated, so the line starting "#$ -q" says "#$ -q all.q,common.q" Of course now you are competing with other groups for the available cores, and you can't just ask Louise and Rui what to do! You might still have to wait a while...

You can copy and paste a runcastep.slurm file from https://docs.archer2.ac.uk/research-software/castep/castep/ The key points are:

#!/bin/bash

# Request 2 nodes with 128 MPI tasks per node for 20 minutes
#SBATCH --job-name=CASTEP

• Replace CASTEP with anything that you will recognize in the job queue.

#SBATCH --nodes=2
#SBATCH --tasks-per-node=128
#SBATCH --cpus-per-task=1
#SBATCH --time=00:20:00

# Replace [budget code] below with your project code (e.g. t01)
#SBATCH --account=[budget code]

• Our budget code is e05-biosoft-sal

#SBATCH --partition=standard
#SBATCH --qos=standard

# Load the CASTEP module, avoid any unintentional OpenMP threading by
# setting OMP_NUM_THREADS, and launch the code.
module load castep
export OMP_NUM_THREADS=1
srun --distribution=block:block --hint=nomultithread castep.mpi test_calc

• Replace test_calc with the STEM of your filenames.

sbatch runcastep.slurm

squeue --me (or squeue to see the whole queue)

scancel

Once CASTEP has been running for a minute or so, you should check the output file.

cif2cell bug

There is a bug in cif2cell, which means that some fractional coordinates which are close to symmetry positions are snapped to the symmetry sites.  This is because CASTEP was written for Inorganic structures, and this is usually desirable.  In the .castep output file, if the number in the

Maximum deviation from symmetry

field is of the order of <1e-10, then you are fine, but if this around 1e-3 or higher then you have fallen foul of this bug.  Louise was confident that the version of cif2cell installed on ib-server and the one in her folder on ib-server which you will have installed on Archer2 have the bug corrected by Rui. However, Louise has just encountered the bug. It is actually relatively simple to fix. You can identify quickly which atoms are affected - they will have a different number of significant figures from the other atoms. The whole atom block can be copied and pasted into Excel, and the sum or difference of atomic positions calculated so you know what fix needs to be made. In Louise's case, two oxygen atoms had x positions of 0.333333333 and 0.666666667. Creating a table of all oxygen positions, and calculating the sum of the x positions, made it simple to see what the relationships should have been and calculate the correct values for these two atoms. Then the .cell file could be fixed manually.

Memory

The memory required per core should be less than 1000 MB (the number under "Disk"). If it's getting close to 1000, then you might need to kill the job and resubmit with more cores.

This is complicatd, and Louise doesn't understand it yet. She thinks that less than 256 MB is okay, and that you can just increase nodes until you get below this. But still there are Out Of Memory errors in some runs.

Energies

Then you will se the SCF convergence for each energy cutoff. Three energies are printed out (they will all be the same), and it is the "Final energy, E" that you need. The next step, with a higher energy cutoff for the basis set will have a better energy. There will also be a section "Total energy corrected for finite basis set" with gives you the energy extrapolated to an infinite basis set. This is only written out for a single point energy.

With a Hartree-Fock calculation, a better description of the basis set will give you a more accurate energy, which will always be lower than the approximation you have. With DFT methods, this is not the case, and the energy extrapolated to the infinite basis set can be slightly higher than with the finite basis set you have selected.

Analysis

Once CASTEP has been run with this k-point spacing, alternative spacings of, say, 0.05 and 0.06 should be used, with a more limited range of cutoff energies (or just 700 eV and finite_basis_corr set to 0 which switches off the next two things) for the basis set.

In the first file, you have scanned through a number of energy values. Generate a plot of cutoff energy in eV vs. Final energy, and see where the plot levels off. If you have more than one polymorph, calculate the energy difference between the polymorphs at each cutoff energy, and see where that converges.

There is a script, cutoffconv that you can use to extract this using:

./cutoffconv *.castep

This gives a file called conv.dat which can be imported into Excel. You need to divide the energy (in eV) by the number of molecules in the unit cell and convert it to kJ mol^-1 (1eV = 96.487 kJ mol^-1). For the analysis of the k-point spacing convergence, you can manually extract the numbers or you can use this script.

If you have more than one crystal structure, you want the difference in energy between the two structures to not change by more than 1 kJ mol^-1 between one step and the next. If you only have one crystal structure, then the difference between one step and the next must be below this threshold. Louise isn't sure whether this threshold changes for bigger molecules.

As the accuracy of the basis set required for a molecule depends on the elements contained within the molecule, it should be safe to only test convergence of the cutoff energy required once for each combination of elements. The following table is a list of the element sets that Louise has tested with the PBE method and the TS dispersion correction (data can be downloaded from here).

Similarly, you should check at what k-point spacing your energy has converged. This does not depend on elements, and should be checked for each chemical system that you run.

Running a cell optimization

You can use the same .cell file, changing the k-point spacing to the value you have selected fromyour convergence tests. You will also need a different .param file.

You can copy /home/rui/CASTEPinput_template/VURYOW_eme-opt.param to your folder

Edit this .param file:

• The "sedc" line is not commented out, so the dispersion correction will be applied. We have found that TS is better for geometry optimization that G06, but it was mentioned recently that TS gives a reproduction closer to room temperature structures than should be expected (0 K is what we are notionally assuming!). MBD* is too costly for geometry optimization in most cases.

• The task is "geometryoptimization"
• Change the cutoff energy to the value you have determined from yout testing
• "finite_basis_corr" is needed for CASTEP to extrapolate to infinite convergence limit, but leave it as 2. There is no stepsize, so CASTEP uses 5. This needs checking, as we don't think that this is calculated for geometry optimizations at the moment.
• "geom_force_tol" is the convergence threshold and 0.001 eV/ang is required if you are going to follow up with phonon calculations. Looser convergence will obviously be faster, but stick with 0.001 until you know what you are doing.
• Keep the "scf_cycles" at 100
• The backup interval is how frequently the wavefunction is backed up in units of seconds (21600 s = 6 h). This is more applicable on Archer2, where jobs are limited to 24 hours.
• "write_cell_structure" and "calculate_stress" need to be "true" so it writes out the final cell and calculates the stress. Stress should be very low.

This will take a lot longer to run than the convergence testing, so think about the number of processors you want to ask for.

On ARCHER2, we get charged by the number of cores we are using (we should probably check this, as the charging has changed and it seems a bit more complicated). The way CASTEP parallelises is that each of the k-points is set to the same number of processors. So if you have a 2*2*3 grid, you should request a multiple of 12 processors for maximum efficiency.

After a geometry optimization

Extract the energy for the final point. One of the output files is STEM-out.cell.

Use:

cell2shx < STEM-out.cell > STEM/whatever.cif

to convert this to a .res file, and then use PLATON to find the correct symmetry.

Single point energy with MBD* dispersion correction

You can use the output .cell file as the starting point for this calculation. You just need to change the final line which specifies the "kpoint_mp_grid" back to the "kpoint_mp_spacing" you had before. Edit the .param file as necessary. When you submit the job, keep an eye on it until the k-point grid is written out. If the grid is not the same as that in the geometry optimization, then stop the job as soon as you notice. You will need to run a geometry optimization starting from the end point (with the same .param file as previously), so that you are using the same k-point grid for both the geometry optimization and subsequent single point energy steps.

It is possible, although very rare, that you will have a structure that fluctuates between two k-point grids when you use the same k-point spacing. In this case, you should probably run with the k-point grid definition, to keep a roughly consistent k-point spacing with the other structures of this molecule. But you MUST make a note that you have done this.

Restarting a timed out job

On ARCHER2, the jobs are limited to 24 hours. As such, you will often find that a larger job will not complete within the time limit and needs to be restarted.

Restarting from the .check file

Louise has done this once or twice in the past, but the default with our .param file is to write out a .check file every 6 hours, and so what you usually have is something 6 hours out of date. There is also a .check_bak, which is the previous .check file, and is about 12 hours out of date. Maybe, if you want this to be an option for you, you should change the backup_interval in your .param file to something more sensible like 3600.

Manually restarting

You can get all the information on the latest geometry from the STEM.castep file, and extract in manually. There are two sections that you need.

Unit cell parameters

Search for the last incidence of "Unit Cell" in your STEM.castep file, and copy the three lines underneath it to the clipboard. Paste it into your STEM.cell file, in the section between %BLOCK LATTICE_CART and %ENDBLOCK LATTICE_CART. Ensure you only have the lattice and not the reciprocal lattice as well, removing the last three columns.

Fractional coordinates

Just below the Unit Cell parameters in your STEM.castep file you will have the Fractional coordinates of the atoms. Copy this block of text to the clipboard. Paste it into your STEM.cell file, in the section between %BLOCK POSITIONS_FRAC and %ENDBLOCK POSITIONS_FRAC. Ensure you only have the atomic type and three numbers, removing the atom number and x characters from the section.

For both these sections, the precision of the numbers in the output STEM.castep file is much lower than what you originally used in your STEM.cell file, but most of these extra values are zeros (in the fractional coordinates section at least), and so this is as good a way as any. Then you just need your original runcastep.slurm and STEM.param files, and you are ready to resubmit the job.

Dryrun calculations

A quick way to check all your input files, before you submit a long job that might not start quickly, is to do a dryrun. This will calculate the ultrasoft pseupotentials, write out the k-point grid, write out the maximum deviation from symmetry (so you can check whether there was a problem with cif2cell snapping an atom to a symmetry site) and give an estimate of the memory requirements.

Do not use the castep command as normal to set up the job, but run the serial version of castep on the headnode:

• castep -d STEM

Edit your runcastep.slurm file so that the final line includes the dryrun flag -d, and looks something like:

• srun --distribution=block:block --hint=nomultithread castep.mpi -d STEM

CASTEP Incremental Optimization

All the instructions in this document are based on information in the CASTEP keywords webpage.

Motivation

If you are working with an experimental crystal structure file, there are a lot of variables to be optimized at once. This can cause problems.

Historically, the periodic electronic structure calculations would be carried out in a number of steps. In most cases this would save time (I think). I don't know whether it has been checked that the resulting crystal structure and energy is the same regardless of route.

Job steps

• Optimize hydrogen atoms
• Optimize all atoms
• Optimize complete cell

What you actually do is constrain the cell and heavy atom positions in the first run, constrain just the unit cell in the second run, and remove all constraints in the final run.

Constraining the unit cell

In your .cell file, there is a block which looks like:

%BLOCK LATTICE_CART
      ang    # angstrom units
      9.048000000000000   0.000000000000000   0.000000000000000
      0.000000000000000   4.870000000000000   0.000000000000000
      0.000000000000000   0.000000000000000  18.262000000000000
%ENDBLOCK LATTICE_CART

After this, add a block that looks like:

%BLOCK CELL_CONSTRAINTS
      0   0   0
      0   0   0
%ENDBLOCK CELL_CONSTRAINTS

The first three numbers of this row of this address the cell lengths, and in this case they will all be fixed. If you want to allow some to vary, replce the "0 0 0" with "1 2 3" to allow them all to optimize freely, or link some variables, such as "1 2 1" would mean that the a and c lengths must stay the same but the b length may vary. The last three numbers deal with the cell angles similarly. Make sure that these do not contain numbers that are also in the cell lengths line (apart from 0), or you will find that your beta angle is constrained to be the same as your c cell length!

Constraining atom positions

In your .cell file, there is a block which looks like:

%BLOCK POSITIONS_FRAC
      O    0.288000000000000   0.553900000000000   0.098000000000000
      O    0.712000000000000   0.446100000000000   0.598000000000000
      O    0.212000000000000   0.053900000000000   0.598000000000000
%ENDBLOCK POSITIONS_FRAC

This contains the fractional coordinates of all the ions in the unit cell. Don't worry about charges - a neutral atom is just an ion with no charge. CASTEP calls them all ions. I'm not sure when this is related to Cartesian coordinates.

After this, add a block that looks like:

%BLOCK IONIC_CONSTRAINTS
      1   O   1   1   0   0
      2   O   1   0   1   0
      3   O   1   0   0   1
      4   O   2   1   0   0
      5   O   2   0   1   0
      6   O   2   0   0   1
      7   O   3   1   0   0
      8   O   3   0   1   0
      9   O   3   0   0   1
      10  C   1   1   0   0
      11  C   1   0   1   0
      12  C   1   0   0   1
      13  C   2   1   0   0
      14  C   2   0   1   0
      15  C   2   0   0   1
      16  C   3   1   0   0
      17  C   3   0   1   0
      18  C   3   0   0   1
%ENDBLOCK IONIC_CONSTRAINTS

The first column is the number of the constraint, the second column is the species of the ion, the third column is the ion number of that type from the previous list, and the fourth, fifth and sixth columns are the fractional coordinates that you want to constrain. As you can see, you need to constrain each atom in each of the three directions separately. This is very tedious, and I do it in Excel.

In the first step of your workflow, constrain the cell and all the non-H atoms, and run an optimization calculation in CASTEP. In the second step, take the output *-out.cell file, and remove constraints on the atoms, but retain the unit cell constraints, and run another optimization. In the final step, take your new output *-out.cell, and remove the unit cell constrains, and run another optimization.

It is possible to optimize the unit cell parameters, but keep the Cartesian coordinates fixed. In this case, it's not clear what would happen, and it's likely that it is actually the fractional coordinates that stay fixed. This would lead to a distortion of the molecule, so the unit cell would probably not change at all. Don't do it.

Considerations of conformation

There are two reasons why you might want to consider the conformational energy of the molecule. Firstly, if you want to calculate the Lattice Energy of a crystal structure, that is relative to the infinitely separated molecules in their lowest energy conformation. Hence you would need to know the energy of the molecule and the energy of the molecule in its lowest energy conformation. Secondly, there is the issue reported by Beran of the inaccuracies of DFT conformational energies. See this paper for more details.

Calculating Lattice Energy

There are a few components to consider here. Firstly, you need to know the intermolecular contribution to the energy. You can get this by subtracting the energy of the isolated molecule in the same conformation as in the crystal from the total crystal energy. To get this, you run a single point energy CASTEP calculation on just a single molecule in a large enough box that the periodicity does not give any molecular interactions. Rui wrote a document on how to do this.

On ib-server, Rui has written a fortran code to assess a molecule in a mol.xyz file for longest interatomic distances, and then put this in an orthorhombic box with a specified distance between molecules in adjacent boxes, writing out a mol.cell file with most of the things you would want in it. The only .f90 Louise can find is /home/rui/temp/singlemol_cellmaker.f90 which has also been copied to /home/lsp/CASTEP/ for safe keeping. The compiled utility is /home/rui/succinic/beta/ncp/singlemol/singlemol_cellmaker/singlemol_cellmaker although alternative versions are available in /home/rui/fluorinePot/ but it is not documented what the differences are.

Rui's documentation suggests that 20 Angstroms is about right, but whatever the limit is that increasing this separation by 2 Angstroms gives a change in absolute energy of less than 0.1 kJ mol^-1 is what you should use. Therefore, run this utility on your mol.xyz file a number of times with different separations and run single point energy calculations to test the convergence.

This utility has not been tested on Archer2.

Secondly, you need to consider how far this molecular conformation differs from the gas phase optimized conformation. Rui and Louise have discussed this. Since you would need to run an optimization calculation in CASTEP, you cannot just put the molecule in a large orthorhombic box and hope it stays there sufficiently that there are no intermolecular interactions from the periodicity. You can fix atomic positions, but if you fix more than one, you risk not fully optimizing the molecule.

Since the box has no symmetry, constraining the molecule to near its middle is unnecessary. The only issue would arise if a molecule rotated within the box. Louise's first optimization barely moved in over 70 steps, so this might not be an issue. Please tell Louise if you want any more notes added here.

An alternative method of calculating the intramolecular energy penalty would be to run a geometry optimization in GAUSSIAN with a "large enough" basis set (as Rui said), and see how much the energy drops in that. But what constitutes a "large enough" basis set? Some of the comments in the paper by Beran may be of use.

Correcting monomer energies

Inaccurate Conformational Energies Still Hinder Crystal Structure Prediction in Flexible Organic Molecules by Greenwell and Beran noted that the conformational energies calculated by DFT with the PBE functional were not accurate enough, and caused the incorrect energy ranking of (particularly conformational) polymorphs. The PBE functional has a tendency to overestimate the conformational stability gained by extended pi conjugation. Their recommendation is to correct the monomer energy, by subtracting the DFT-PBE calculated monomer energy from the crystal energy and adding on the MP2D monomer energy calculated separately.

Unfortunately, MP2D is not available in GAUSSIAN at the moment.

From the SI of the paper by Beran:

DFT monomer energies: To evaluate the monomer correction, B86bPBE-XDM single point energies were computed for each monomer in the asymmetric unit with planewave B86bPBE-XDM in a large periodic box with 20 Angstroms of vacuum spacing in all directions. As discussed in Section 1.2, this is typically sufficient to converge the relative monomer conformational energies to within about 0.1-0.2 kJ/mol (which is in good agreement with what Rui found). The monomer geometries were extracted directly from the crystal structures. These calculations were performed at the gamma point with the same planewave cutoff and PAW potentials as for the crystal geometry optimizations.

MP2D monomer energies: Complete basis set MP2D6 monomer energies were computed on the same monomer geometries using PSI4. 7 Density fitting was employed for both the Hartree-Fock and MP2 portions, and core electrons were frozen in the MP2 portion of the calculation. Extrapolation of the correlation energy to the complete basis set (CBS) limit was performed from the aug-cc-pVTZ and aug-cc-pVQZ basis sets using the standard two-point extrapolation approach.

We do not have PSI4, and do not have MP2D in GAUSSIAN. Running a simpler MP2 calculation in GAUSSIAN with one of these augmented basis sets could take a very long time, even though it is just a single point energy being calculated.

Another interesting paper on this subject is Overcoming the difficulties of predicting conformational polymorph energetics in molecular crystals via correlated wavefunction methods by Greenwell et al.