Create an input file for FDMNES from a cif or poscar file

cif2fdmnes.py

#!/usr/bin/env python
#
# Create an input file for FDMNES from a cif or poscar file
# Davide Ceresoli <dceresoli@gmail.com>, 27/02/2024
#

import os, sys, argparse
import ase.io
import ase.data


#======================================================================
# parse CLI arguments
#======================================================================
parser = argparse.ArgumentParser(description='Create input for FDMNES calculations')

parser.add_argument('-a', '--absorber', dest='absorber', action='store', required=True,
                    help='absorber element')

parser.add_argument('-e', '--edge', dest='edge', action='store', default='K',
                    help='absorption egde (default: K)')

parser.add_argument('--rmax', dest='rmax', action='store', default=7.0,
                    help='maximum radius of cluster (default: 7 Å)')

parser.add_argument('--range', dest='range', action='store', default='-20,50',
                    help='energy range (default: -20,50)')

parser.add_argument('--scf', dest='scf', action='store_true', default=True,
                    help='perform SCF calculation (default: True)')

parser.add_argument('-m', '--method', dest='method', action='store', default='MT',
                    help='method: MT or FD (default: MT)')

parser.add_argument('--tddft', dest='tddft', action='store_true', default=False,
                    help='use TDDFT (default: False)')

parser.add_argument('--quadrupole', dest='quadrupole', action='store_true', default=False,
                    help='calculate quadrupole transitions (default: False)')

parser.add_argument('filename', action='store', help='cif or poscar filename')

args = parser.parse_args()
#print(args)


#======================================================================
# read file
#======================================================================
try:
    cry = ase.io.read(args.filename)
except:
    print('error: opening or parsing {0}'.format(args.filename), file=sys.stderr)
    sys.exit(1)

if args.absorber not in cry.get_chemical_symbols():
    print('error: atom {0} not found'.format(args.absorber))
    sys.exit(1)
    

#======================================================================
# create output file
#======================================================================
base = os.path.basename(args.filename)
prefix = os.path.splitext(base)[0]
out = 'FDMNES_' + prefix + '.txt'

with open(out, 'wt') as f:
    print('! FDMNES input file for {0}'.format(args.filename), file=f)
    print('! created by cif2fdnmnes.py', file=f)
    print('', file=f)
    print('Filout\n  ./out\n', file=f)

    emin, emax = args.range.split(',')
    print('Range                       ! Energy range of calculation (eV). Energy of photoelectron relative to Fermi level', file=f)
    print('  {0} 0.1 {1}\n'.format(emin, emax), file=f)

    print('Radius                       ! Radius of the cluster where final state calculation is performed', file=f)
    print('  {0}\n'.format(args.rmax), file=f)

    print('Edge                         ! Threshold type', file=f)
    print('  {0}\n'.format(args.edge), file=f)

    if args.scf:
        print('''SCF                          ! Self consistent solution
Full_atom                    ! Better SCF convergence
Delta_E_conv
  0.010
N_self
  1000
P_self
  0.01
''', file=f)

    if args.method.lower() == 'mt':        
        print('Green                        ! Muffin tin potential', file=f)
    elif args.method.lower() == 'fd':
        print('!Green                        ! Finite differences', file=f)
        raise 'TODO'
    else:
        raise NotImplementedError
    
    if args.quadrupole:
        print('Quadrupole                   ! Allows quadrupolar E1E2 terms', file=f)

    if args.tddft:
        print('TDDFT                        ! Use Timde Dependent DFT', file=f)
        raise 'TODO'

    print('''Relativism
!magnetism                    ! performs magnetic calculations
Density                      ! Outputs the density of states as _sd1.txt
Sphere_all                   ! Outputs the spherical tensors as _sph_.txt
Cartesian                    ! Outputs the cartesian tensors as _car_.txt
energpho                     ! output the energies in real terms
Convolution                  ! Performs the convolution
All_conv

rxs                          ! Resonant x-ray scattering at various peaks
''', file=f)
    
    print('Atom ! s=0,p=1,d=2,f=3, must be neutral, get d states right by moving e to 2s and 2p sites', file=f)
    species = set(cry.get_chemical_symbols())
    for sp in species:
        print('{0} 0 ! {1}'.format(ase.data.atomic_numbers[sp], sp), file=f)
    print('', file=f)
    
    print('Crystal                      ! Periodic material description (unit cell)', file=f)
    a, b, c = cry.cell.lengths()
    alpha, beta, gamma = cry.cell.angles()
    print('  {0:12.6f} {1:12.6f} {2:12.6f} {3:12.6f} {4:12.6f} {5:12.6f}'.format(a, b, c, alpha, beta, gamma), file=f)

    species = list(species)
    
    atoms = cry.get_chemical_symbols()
    frac = cry.get_scaled_positions()
    for i in range(len(cry)):
        ind = species.index(atoms[i])
        if atoms[i] == args.absorber:
            print('{0:3d} {1:20.10f} {2:20.10f} {3:20.10f} ! {4:3d} {5}'.format(ind+1,
                  frac[i][0], frac[i][1], frac[i][2], i, species[ind]), file=f)

    for i in range(len(cry)):
        ind = species.index(atoms[i])
        if atoms[i] != args.absorber:
            print('{0:3d} {1:20.10f} {2:20.10f} {3:20.10f} ! {4:3d} {5}'.format(ind+1,
                  frac[i][0], frac[i][1], frac[i][2], i, species[ind]), file=f)

    print('End\n', file=f)


with open('fdmfile.txt', 'wt') as f:
    print('1', file=f)
    print(out, file=f)
    
print('files {0} and fdmfile.txt written successfully!'.format(out))