Computational materials science protocols

• Molecular crystal CIF refinement
• Molecular crystal geometry optimization
• Determination of molecular crystal geometry from MD snapshot of a large supercell
• Electronic coarse-graining of insulators
• Calculation of coarse-grained electronic Hamiltonian for molecular crystals
• Creation and parametrization of tight-binding models

Molecular crystal CIF refinement

• Use MolMod/exam_cif2xyz.mw, see MolMod package.
• Create CIF draft according to this template.
• Preprocess CIF draft:
  • Generate entire molecule and move it (to the origin if possible).
  • Change symmetry settings.
  • If there are multiple symmetry-inequivalent molecules in the unit cell, create separate CIF for each of them and process separately. If molecules are naturally grouped (solute+solvent, charge-transfer pair), they can be processed as a single meta-molecule. All molecules must have identical and separated numbering.
  • Make molecules connected.
  • Reorder atoms for consistency with the rest of data. In complex cases extract molecule, reorder it with MolMod/exam_Geometry.mw/Molecules: Partition molecule and identify fragments code, and use the result as reference for reordering.
  • If symmetry-equivalent conformers are present, chose the conformer consistent with the rest of data, i.e. it should differ from a reference monomer by a translation only, and, if needed, a proper rotation.
  • Symmetrize molecule at fixed digits.
  • Add missing or replace ill-positioned hydrogens with PM7-optmized ones using MolMod/exam_Geometry.mw/Crystals: add missing hydrogens.
• Provide generators g11,g12,...,g21,g22,... explicitly and modify them in CIF to get the optimally packed unit cell (gm1 generates m-th molecule, it should be a proper rotation if possible).
• Save XYZ (example) and keep also thus created CIF (example).
• Store the full symmetry information in XYZ-file including SG,nm,nas,orbs defined here.

Molecular crystal geometry optimization

• Use MolMod/exam_vas.mw, see MolMod package.
• Start with refined CIF (example) converted to XYZ (example).
• Each disorder group requires separate optimization.
• Use PBE-D3/PAW as optimal trade-off between accuracy and complexity.
• Use fast SCF iterators like ALGO=F or VF.
• If initial geometry is problematic preoptimize at fixed cell.
• Optimize cell+coordinates with k-spacing≈0.2 Å^-1 (expected numerical accuracy is 0.1 meV/atom), ECUT=900 eV, and PREC=Accurate.
• Fix all used k-grids in fopts.ini (example) according to the following guidelines:
  • Note that for small-molecule organic semiconductors typical values of K₀ and V₀ are 10 GPa and 10 ų resulting in 0.3 meV/atom energy difference at 3% volume change (example).
  • k-spacing≈0.2 Å^-1 is usually enough for smooth EoS with numerical accuracy better than 0.1 meV/atom.
  • Use comparable-size k-grids for different polymorphs.
  • Use the same number of k-points along topologically equivalent directions.
  • MBD, highly dispersive bands, and flat PES might require k-spacing as low as 0.15 iAo.
• Gaussian smearing with default SIGMA=0.2 eV might be inaccurate for narrow-gap systems like hexacene with PBE functional or charge-transfer compounds.
• Set convergence criteria according to the following guidelines:
  • Default EDIFF=1e-4 and tightened EDIFFG=1e-4 (0.1 meV) should give geometry with 1% accuracy and final gradients of the order of tens meV/Ao.
  • For more robust convergence, option tol=-1 sets EDIFF=1e-5, EDIFFG=1e-4.
  • Set EDIFFG=1e-5 to improve accuracy and get gradients below 10 meV/Ao (option tol=-2 sets EDIFF=1e-6, EDIFFG=1e-5).
  • For numerically precise geometry optimization set EDIFF=1e-6 and EDIFFG=-1e-3 (1 meV/Ao).
• For large unit cells, accuracy can be increased gradually in the following order:
  1. Increase ECUT from 400.
  2. Increase PREC from normal to accurate.
  3. Increase k-grid from minimal reasonable (such as Γ-point only, adding more points for highly dispersive directions).
• If feasible, get final geometry by 5-point volume scan with 3% volume change step at ECUT=600 (example, if number of geo iterations exceed 30-40, additional run is needed to converge all volume scan points with precision of fractions of meV/atom). Otherwise reoptimize atomic positions at ECUT=400 (minimum for organics).
• Decrease k-spacing to 0.06-0.1 Å^-1 to calculate final total energy and eDOS.
• Check consistency of all calculations by volume scan plot and EoS fit (especially K'₀).
• Symmetrize or at least check symmetry at each run. Check that CIF is not corrupted. Check supersymmetry and symmetry breaking.
• If feasible calculate Γ-point phonons and elastic tensor to check the stability (example).

Determination of molecular crystal geometry from MD snapshot of a large supercell

• Use ITIC4F/dat/main.mw/Crystal
• It is assumed that all molecules are connected.
• Group molecules into translationally invariant sets.
• Determine cell translations using a reference atom and then fold cells.
• Superimpose molecules to inspect conformational disorder. It is assumed that the conjugated backbone is not conformationally disordered.
• Determine disorder groups for side chains: 1) Determine disordered part of side chains and type of disorder using both cartesian and internal coordinates. 2) Analyze disorder by histograms in proper coordinates. 3) Use cluster analysis (density-based such as DBSCAN and Gaussian distribution-based) to determine disorder groups. 4) Discard nonessential groups.
• Fix a disorder group and calculate average molecule. Identify atoms whose average position is meaningless and determine their mean positions in internal coordinates. In particular, aliphatic side chains often can be defined only by fixing C-C-C-C dihedrals to their most probable value consistent with average position of the chain in cartesian coordinates.
• Check intramolecular symmetry and symmetrize.
• Determine unit cell origin, check crystal symmetry, and symmetrize.
• Optimize geometry with any light but accurate method to redetermine the space group.

Electronic coarse-graining of insulators

• Use QFT/Models/ECG/ECG.mw or MolMod/LocalizeMO.mw, see MolMod package.
• PBC are supported only at Gamma point (supercell approach).
• Fragmentize system, identify repeating fragments (monomers), and save to file following quantities: 1) fragments as F=list(list(atom#)); 2) bonds between fragments as B=list([atom#,atom#]); 3) terminal bonds as T=list([atom#,atom#]) with the 1st atom belonging to a fragment; 4) fragments to monomers map as f2m=list(monomer#); 5) monomers to fragments map as m2f=list(list(fragment#)); 6) reference fragment for each monomer as m2fref=list(fragment#); 7) optionally monomer names as MN=list(string); 8) optionally symmetry information (translation period, asymmetric unit).
• Calculate MOs and save out,evl,evc,S (or evr).
• Select subspace of MOs for coarse-graining as accurate as possible, though small number of irrelevant MOs can be discarded later during orbital localization.
• Prelocalize by geometry or other methods such as NBO or spatial localization, so that you will get a list of pre-LMOs.
• Localize by projection onto pre-LMOs, filter out redundant MOs and pre-LMOs, iterate until convergence (usually the second iteration is either final or meaningless). Correct LMO signs by a sign function, or by reference orbitals, or manually, or applying a combination of these methods. Save the obtained LMOs as evc,evr,h1e (optionally s1e) binary files and plain-text records containing: 1) copy of fragmentation information (see above); 2) basis set fragmentation as K=list(list(basisfunction#)); 3) LMOs to fragments map as o2f=list(fragment#); 4) fragments to LMOs map as f2o=list(list(LMO#)); 5) list of coarse-grained MOs MOs=list(MO#); 6) localization method as method=string.
• There is a code to print the Hamiltonian and overlaps. Such files are labeled with "_HS". The output contains list of LMOs for which the matrix elements are printed. The first row and column of the printed table contains LMO energies in eV. The second and third columns give fragment and orbital numbers. The rest of the table are matrix elements (in eV for Hamiltonian).
• If needed, create reference LMOs for a fragment by cutting the geometry and basis. Save out,evc,evr,h1e.

Calculation of coarse-grained electronic Hamiltonian for molecular crystals

• Use MolMod/_exam_Heff.mw, see MolMod package.
• Current implementation is limited to: 1) the same ECG basis for all monomers, 2) nondegenerate electronic states. For π-conjugated systems, condition (1) implies that all monomers have the same π-conjugated backbone and are interconnected by single bonds.
• Load crystal structure, monomerize it, then rationalize site positions and initialize FGS subpackage of FiniteGroups package. Important input parameters:
  • SG space group,
  • M matrix of translation vectors augmented to 3D,
  • nm number of molecules/polymers in the unit cell,
  • F fragmentation of each molecule/polymer.
• Step 1: Generate symmetry-unique multimers or supercell here or elsewhere (e.g. Polymers.mw) as a set of xyz-files. In the simplest approximation, it is enough to consider only 1) nearest-neighbor dimers (up to 5~Ao between π-conjugated atoms), 2) hydrogen-passivated π-conjugated backbones instead of entire molecules.
• Step 2: Calculate MOs of multimers/supercells and all translationally unique molecules. For each calculation save out,evc,evr,evl files for frontier orbitals (keep enough number of MOs per monomer or save entire evc and s1e files).
• Step 3: Generate translationally-unique ECG-projection basis for each translationally unique monomer and save it as out,evc,evr,evl files. This can be done by ECG of all translationally-unique molecules using some reference orbitals. Reference orbitals can be taken either from MO of the monomer or from ECG of an oligomer. If monomer=molecule, Step 3 can be omitted (in this case the molecule's HOMO/LUMO serves as the basis).
• Step 4: Calculate effective ECG Hamiltonian using MultimerH program and save it in this format. There is an extensive helpful information in several log-files (_mm,_HL,_HS,_HC) generated by the code. Important input parameters:
  • typ,nof electron/hole and number of orbitals per monomer to define the projection basis,
  • rmax max distance (in lattice coordinates) from the multimer center to the center of two fragments for which coupling is calculated, reasonable values are from 0 to 1, the default is 0.55,
  • nf3m number of fragments per each end of each molecule to ignore in calculation of Hamiltonian elements, reasonable values are from 1 to 3,
  • CDmin,tmin minimal contact distance and value of couplings to keep matrix elements (couplings are discarded if CD>CDmin and t<tmin), to strictly obey symmetry set CDmin<0,
  • Hsortf,Hpickid sorting function for symmetry-identical matrix elements (the most relevant one goes first) and how to pick the best value (currently two options are supported: first (default) or mean).
• Step 5: Print transfer integrals, calculate hopping amplitudes and band structure. The report is saved to txt-file.

Creation and parametrization of tight-binding models

• Use QFT/Models/TB.mw
• Create lattice geometry initialization file in this format supplemented with CIF.
• Create TB Hamiltonian initialization file in this format, including list of symmetry-unique transfer integrals, visualized with CIF and picture.
• Save symbolic Hamiltonian in this format.
• Parameterize this TB Hamiltonian from supercell or cluster calculations, saving the result in the TB Hamiltonian file.
• Derive eigenvectors and eigenvalues at high-symmetry k-points and k-paths.
• Derive a minimal TB model and plot its band structure and DOS.
• If relevant, parameterize the TB model from high-symmetry k-points.
• Derive hopping amplitudes tensor and inverse mass tensors at selected eigenstates.
• Study convergence of TB models wrt number of parameters and supercell/cluster size.