# Ab initio random structure searching This example illustrates crystal structure prediction with ab initio random structure searching ([AIRSS](https://www.mtg.msm.cam.ac.uk/Codes/AIRSS)). The first step is to generate random structures, having sensible initial volumes and interatomic distances but otherwise unconstrained (at least in this example). Then, one relaxes each to a local minimum in the potential energy landscape and analyzes the results. The procedure is simple but effective for predicting crystal structures entirely from first principles. ## Seed file The AIRSS package generates structures based on a seed file, assumed here to be named ``Mg.cell``. For a magnesium search over structures with 2-4 atoms, a reasonable seed would be #SPECIES=Mg%NUM=2-4 #TARGVOL=2.5-3.0 #MINSEP=1.5 The `TARGVOL` keyword sets the range of volumes (per atom) permitted in the random structures. Any structures with atoms closer than `MINSEP` are rejected. ## Structure relaxations To perform the search we require a `relax.py` script, which accepts a random structure as input and optimizes the geometry with ``profess``. An example script is provided here. import numpy as np import os import sys sys.path.append(os.path.join(os.path.dirname( os.path.abspath(__file__)), '../../../../external/ase/')) from ase import Atoms import ase.io import profess import ase_tools to_eV_per_A3 = 1.0 / 160.2177 # from GPa to_Ha_per_B3 = 1.0 / 29421.02648 # from GPa #------------------------------------------------ # Collect command line input and load random structure pressure = float(sys.argv[1]) seed = sys.argv[2] element = seed.split('-')[0][2:] atoms = ase.io.read(seed + '.xyz') box_vectors = atoms.cell[...] xyz_coords = atoms.get_positions() # Relax in stages to refresh the grid periodically stage = 0 num_converged = 0 while stage < 100: stage += 1 print('---------- STAGE {} ----------'.format(stage)) steps_per_stage = 20 if stage <= 3: # use short stages at first steps_per_stage = 3 # Build system and relax geometry energy_cutoff = 600 potential = '../potentials/' + element.lower() + '.gga.recpot' system = ( profess.System.create(box_vectors, energy_cutoff, ['a','ev']) .add_ions(potential, xyz_coords, 'a') .add_electrons() .add_generic_nonlocal_a_functional( # kinetic energy (5+np.sqrt(5))/6, (5-np.sqrt(5))/6, lambda x: np.exp(x), lambda x: np.exp(x)) .add_hartree_functional() # electrostatic energies .add_ion_electron_functional() .add_ion_ion_interaction() .add_perdew_burke_ernzerhof_functional() # exchange-correlation ).minimize_energy() converged = ase_tools.minimize_forces_stress( system, 'BFGSLineSearch', tol=2e-4, steps=steps_per_stage, scalar_pressure=pressure*to_Ha_per_B3) box_vectors = np.array(system.box_vectors('a')) xyz_coords = np.array(system.ions_xyz_coords('a')) # If converged twice, output final structure and fake castep file if converged: num_converged += 1 if num_converged == 2: atoms = Atoms(element+str(xyz_coords.shape[0]), positions=xyz_coords, cell=box_vectors), ase.io.write(seed + '-out.xyz', atoms, 'extxyz') volume = system.volume('a3') pressure = system.pressure('gpa') enthalpy = system.enthalpy('ev') with open(seed + '.castep', 'w') as f: f.write("Current cell volume = {:25.15f} A**3\n".format(volume)) f.write("* Pressure: {:25.15f}\n".format(pressure)) f.write("Python: Final Enthalpy = {:25.15f} eV\n".format(enthalpy)) break ## Performing the search One initiates a search via the command line airss.pl -python -max 10 -seed Mg where `-max` sets the number of initial random structures. ## Analyzing the results of a search Other AIRSS tools analyze the results. The command ca -r might produce (without the column labels) identifier pressure volume energy num space group Mg-12790-5343-2 -0.00 23.099 -24.221 2 Mg P63/mmc 1 Mg-12790-5343-3 0.00 23.096 0.000 4 Mg P63/mmc 1 Mg-12790-5343-1 0.00 23.094 0.001 4 Mg P63/mmc 1 Mg-12790-5343-6 0.01 23.090 0.002 2 Mg P63/mmc 1 Mg-12790-5343-5 0.00 23.122 0.002 3 Mg R-3m 1 Mg-12790-5343-7 -0.00 23.214 0.009 3 Mg I4/mmm 1 Mg-12790-5343-8 -0.00 23.215 0.009 3 Mg I4/mmm 1 Mg-12790-5343-4 0.00 23.215 0.009 2 Mg I4/mmm 1 Mg-12790-5343-9 0.00 23.212 0.009 4 Mg I4/mmm 1 Mg-12790-5343-10 0.01 23.210 0.010 3 Mg I4/mmm 1 Each structure appears once, ranked by energy. The lowest-energy structure, having space group P63/mmc, is the HCP structure investigated in another example. Alternatively, for these same results, the command ca -u 0.01 -r would produce Mg-12790-5343-2 -0.00 23.099 -24.221 2 Mg P63/mmc 4 Mg-12790-5343-5 0.00 23.122 0.002 3 Mg R-3m 1 Mg-12790-5343-7 -0.00 23.214 0.009 3 Mg I4/mmm 5 Similar structures (based on a tolerance value 0.01) are now united into a single row. The rightmost column indicates frequency---in this case, the HCP structure was discovered 4 times, and R-3m and I4/mmm structures were found 1 and 5 times, respectively. The same `ca -u 0.01 -r` command, after a search of 100 structures having between 2-8 atoms, might yield Mg-19149-2215-4 -0.01 23.100 -24.222 2 Mg P63/mmc 26 Mg-19149-2215-8 -0.01 23.128 0.002 3 Mg R-3m 12 Mg-19149-2215-35 -0.01 23.218 0.009 2 Mg I4/mmm 56 Mg-19149-2215-77 -0.00 23.241 0.010 3 Mg Fmmm 2 Mg-19149-2215-47 -0.00 23.243 0.010 2 Mg Fmmm 2 Mg-19149-2215-9 -0.00 23.268 0.024 8 Mg C2/m 2 The larger search will require several hours on a laptop.