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

Reference: JCTC 19, 8481 (2023)
Code: exam_cif2xyz.mw from 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 needed.
  • 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

Reference: JCTC 19, 8481 (2023)
Code: exam_vas.mw from 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

Reference: JCP 159, 024107 (2023)
Code: ITIC4F/dat/main.mw/Crystal (available per request)

• 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 the 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

Reference: JCP 159, 024107 (2023) (initial development)
Code: LocalizeMO.mw from MolMod package, QFT/Models/ECG/ECG.mw (available per request)

• 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

Reference: in preparation, initial development is in Annu Rev Phys Chem 66, 305 (2015)
Code: _exam_Heff.mw (under revision) from MolMod package

• Current implementation is limited to: 1) the same ECG basis for all monomers, 2) nondegenerate electronic states (in crystal symmetry). For π-conjugated systems, condition (1) implies that all electronically active 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). Multimers are uniquely defined by lists of mt=[m,t1,t2,t3], meaning m-th molecule/polymer translated by t, or equivalent string code used in filenames. Polymers can be approximated by long oligomers, see here. For transfer integrals it is usually enough to consider only nearest-neighbor dimers and replace molecules by hydrogen-passivated π-conjugated backbones. For onsite energies at least one shell of nearest neighbors is required. A multimer should be centrosymmetric wrt to the target matrix element. Nearest neighbors can be identified by the contact distance between π-conjugated atoms with a reasonable cutoff of 5 Ao or at the nearest minimum of the corresponding distribution function.
• 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). Note that the computational bottleneck of the entire protocol is in the DFT calculation of the largest multimer.
• 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 achieved by the projection of MOs of the smallest multimers onto reference orbitals (it is important for such projection to preserve the symmetry between all symmetry-identical monomers). Important output parameter here is minocc2 which is the minimum occupation in LocalizeMO output: it should be close to 1 (values below 0.7 require an inspection of the obtained ECG basis). Reference orbitals can be taken either as MOs of the monomer or from ECG of some reference molecule, for which one can take one of the smallest multimers. For nearly degenerate orbitals, a large enough multimer might be needed to obtain the correct reference orbitals (e.g. triphenyltriazine). If monomer=molecule and nof=1, Step 3 can be omitted since HOMO/LUMO can serve as the basis implicitly, but for record-keeping it is better to store the ECG basis explicitly.
• Step 4: Calculate effective ECG Hamiltonian using MultimerH program and save it in this format. 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).
  There is extensive helpful information in several log-files (_mm,_HL,_HS,_HC) generated by the code. Important output parameter here is minW which is the minimum weight in MultimerH output: small values mean that the multimer MOs are not representable in the ECG basis, which is usually caused by bands entanglement (currently used algorithm cannot disentangle bands, possible solution is to expand ECG basis or minimize multimers). In the HS-list, all entries labeled with *! must be inspected (these couplings will be included but their value strongly depends on multimer). In the HC-list, large deviations for couplings near multimer centers indicate nonoptimal choice of multimers or incomplete ECG basis.
• The generated Hamiltonian can be analyzed by QFT/Models/TB.mw code.

Creation and parametrization of tight-binding models

Reference: see here, initial development is in JPCL 12, 4674 (2021)
Code: QFT/Models/TB.mw (available per request)

• 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.