berrypi

#!/usr/bin/env python
'''
The main program that should be run from a command line

  $ berrypi (options)

See

  $ berrypi -h

for details.
'''

from __future__ import print_function # Python 2 & 3 compatible print function
import os, os.path
import sys, shutil
import argparse # parse line arguments
import subprocess
import numpy # math library
import re
import functools # needed for functools.reduce()

import errorCheck as b_PyError
import submoduleProcess as b_PySubProcess
import calculations as b_PyCalc

import config; from config import BERRY_DEFAULT_CONSOLE_PREFIX as DEFAULT_PREFIX
import mmn2pathphase
import testerror; from testerror import testerror
import rmerror; from rmerror import rmerror
from os import environ


class Logger(object): # redirect stdout to case.outputberry file and console
    def __init__(self):
        self.terminal = sys.stdout
        # case mane is taken from structure_name
        outputFileName = structure_name + '.outputberry'
        self.log = open(outputFileName, "w")
        print("Output will also be written in " + outputFileName)

    def write(self, message):
        self.terminal.write(message)
        self.log.write(message)  

    def flush(self):
        #this flush method is needed for python 3 compatibility.
        #this handles the flush command by doing nothing.
        #you might want to specify some extra behavior here.
        pass    

def changeToWorkingDirectory(path):
    os.chdir(path)

def postProcBerryPhase(configFile, phasesRaw, spCalc):
    '''
    Pass a config file to it in order to perform the final calculation
    on the simulation. Returns the main calculation as the class object
    '''
    #set up path, case name, file extensions
    structurePath = configFile['Structure Path']
    changeToWorkingDirectory(configFile['Structure Path'])
    structName = configFile['Structure Name']
    print(DEFAULT_PREFIX + "Performing main calculation on " + structName)
    structFileExtensions = {
        'struct' : '.struct',
        'scf2' : '.scf2',
        'inc' : '.inc',
        }

    #perform the main calculation
    mainCalculation = b_PyCalc.MainCalculationContainer(
        phases = phasesRaw,
        sp = spCalc,
        orb = orbCalc,
        so = soCalc,
        kmesh = configFile['K-Mesh Divisions'],
        file_struct = structName + structFileExtensions['struct'],
        file_inc = structName + structFileExtensions['inc'],
        file_scf = structName + structFileExtensions['scf2'],
        )
    finalPolarizationValue = mainCalculation()
    return mainCalculation
        

def rawBerryPhase(configFile):
    '''
    Runs the automation process. When a command fails, returns False,
    otherwise it returns True. This allows for several automation
    processes to be batched.
    '''
    print(DEFAULT_PREFIX + \
      "Starting BerryPI for", configFile['Structure Name'] )

    # get calculation parameters and options
    runLAPW = configFile['Perform LAPW']

    # get bands
    bands = configFile['Bloch Band']
    if bands == None: # create an empty array if band range is not set
        bands = []
    
    #get the starting directory and switch to it
    newWorkingDirectory = configFile['Structure Path']
    changeToWorkingDirectory(newWorkingDirectory)

    #ensure background env is also in the right directory
    os.system('cd ' + newWorkingDirectory)
    print(DEFAULT_PREFIX + 'New working directory: ' + newWorkingDirectory)

    # setup wien2k, berrypi path and python path
    w2kpath = os.getenv('WIENROOT')
    if w2kpath == None:
        print('Check the existance by executing "echo $WIENROOT" '
                'or set it up by "export WIENROOT=/path/to/wien2k"')
        raise RuntimeError('Unable to find the WIENROOT enviroment variable')
    else:
        print(DEFAULT_PREFIX + 'w2kpath = ' + w2kpath)

    pypath = sys.executable # get that to python
    print(DEFAULT_PREFIX + 'pypath = ' + pypath)
    # check python version
    if sys.version_info[0] <= 2: # only major python version 3 is supported
        raise RuntimeError('Python 2 is not supported')
        
    print(DEFAULT_PREFIX, 'python major version =', sys.version_info[0])
    print(DEFAULT_PREFIX, 'python minor version =', sys.version_info[1])
    print(DEFAULT_PREFIX, 'python micro version =', sys.version_info[2])
        
    bppath = w2kpath + "/SRC_BerryPI/BerryPI"
    if not os.path.isdir(bppath): # check if BerryPI directory exists
        print(DEFAULT_PREFIX, 'bppath =', bppath)
        raise FileNotFoundError('BerryPI directory does not exist')
    else:
        print(DEFAULT_PREFIX, 'bppath =', bppath)


    # print version
    with open(bppath+'/version.info', 'r') as versionFile: # open version file
        print(DEFAULT_PREFIX + versionFile.read())

    print(DEFAULT_PREFIX + "Numpy version: " + numpy.version.version)
    
    # Make sure certain files (prerequisites) exist
    structureName = configFile['Structure Name']
    print(DEFAULT_PREFIX + "Checking prerequisite files for " + structureName)
    if not wCalc: # any calculation except for Weyl points/Wilson loop
        structFileExtensions = {
            'struct' : '.struct',
            'inc' : '.inc'
            }
    else: # case.inc files are not needed for Weyl points/Wilson loop
        structFileExtensions = {
            'struct' : '.struct',
            'klist' : '.klist'
            }

    for extension in structFileExtensions.values():
        tempFilename = structureName + extension
        print(DEFAULT_PREFIX + "Checking existance of",tempFilename)
        if not b_PyError.fileExists(tempFilename):
            print( DEFAULT_PREFIX +\
                "{} does not exist, cannot finish calculation".\
                format(tempFilename))
            exit(2)
        else:
            print(DEFAULT_PREFIX + '-- OK')
    if pCalc: # parallel calc.
        tempFilename = '.machines'
        print(DEFAULT_PREFIX + "Checking existance of",tempFilename)
        if not b_PyError.fileExists(tempFilename):
            print( DEFAULT_PREFIX +\
                "{} does not exist, cannot finish calculation".\
                format(tempFilename))
            exit(2)
        else:
            print(DEFAULT_PREFIX + '-- OK')

    # Check SCRATCH folder in parallel calculations to make sure that
    # it is accessible from all nodes
    if pCalc and (environ.get('SCRATCH') is not None): # parallel calc. and SCRATCH is set
        print(DEFAULT_PREFIX + "You run parallel calculation with SCRATCH variable set")
        scratch = environ.get('SCRATCH')
        print(DEFAULT_PREFIX + f"SCRATCH = {scratch}")
        abspathscratch = os.path.abspath(scratch)
        print(DEFAULT_PREFIX + f"The absolute path to SCRATCH is {abspathscratch}")
        print(DEFAULT_PREFIX + f"The absolute path to case folder is {newWorkingDirectory}")
        if not(abspathscratch == newWorkingDirectory):
            print(DEFAULT_PREFIX + "Checking the number of distinct "+\
                    "hosts in the .machines file")
            cmd = "grep : .machines | grep -Ev "+\
                    "'#|omp_|lapw0|lapw2_|dstart|nlvdw|nmr|granular|"+\
                    "residue|extrafine|balance|lapw2_vector' "+\
                    "| sed 's/: /:/' | uniq |cut -d: -f2"
            print(DEFAULT_PREFIX + f"Trying to execute: {cmd}")
            cmdresult = subprocess.run([cmd], capture_output=True, text=True, shell=True)
            msgwarn = "WARNING: "+\
                    "Unable to determine the list of distinct "+\
                    "hosts in the .machines file. The calculation will proceed, "+\
                    "but it can fail. If it does, please set SCRATCH to the case "+\
                    "folder."
            hostlist = cmdresult.stdout.splitlines() # create list of unique hosts
            print(DEFAULT_PREFIX +\
                    f'Your list of unique hosts from .machines file is {hostlist}')
            numdisthosts = len(hostlist) # number of fosts
            if cmdresult.returncode == 0:
                if numdisthosts == 0:
                    print(DEFAULT_PREFIX + msgwarn)
                elif numdisthosts == 1:
                    msg = "The SCRATCH does not point to the case folder. "+\
                            "But, we should not expect any problems with "+\
                            "locating parallel vector files for w2w since "+\
                            "you run calculations on one machine."
                    print(DEFAULT_PREFIX + msg)
                elif numdisthosts > 1:
                    msg = "It seems that you run a parralel calculation using "+\
                            f"{numdisthosts} distinct hosts {hostlist} "+\
                            "and your SCRATCH variable does not point to the case "+\
                            "folder. In this case there will be problems with "+\
                            "locating parallel vector files for w2w."+\
                            "Please set SCRATCH to the case folder or "+\
                            "modify .machines to run all parallel calculations "+\
                            "on one machine and then rerun BerryPI."
                    raise RuntimeError(msg)
            else:
                print(DEFAULT_PREFIX + msgwarn)
        else:
            msg = "The SCRATCH points to the case folder. "+\
                    "We should not expect any problems with "+\
                    "locating parallel vector files for w2w."
            print(DEFAULT_PREFIX + msg)


    ######################################
    # copy structure file                #
    # cp $filename.struct $filename.ksym #
    ######################################
    fileNameSrc = structureName + '.struct'
    fileNameDest = structureName + '.ksym'
    shutil.copy(fileNameSrc, fileNameDest)
    print(DEFAULT_PREFIX + "Copied " + fileNameSrc + " to " + fileNameDest)


    ###########################################
    # copy klist file                         #
    # cp $filename.klist $filename.klist_orig #
    ###########################################
    if runLAPW and not wCalc:
        fileNameSrc = structureName + '.klist'
        fileNameDest = structureName + '.klist_orig'
        print(DEFAULT_PREFIX + "Trying to copy " + fileNameSrc + " to " + fileNameDest)
        if os.path.isfile(fileNameDest): # check if $filename.klist_orig already exists
            msg = "file " + fileNameDest + " already exists. Copy aborted."
            print(DEFAULT_PREFIX + msg)
            resoreKlist = False
        else:
            resoreKlist = True
            shutil.copy(fileNameSrc, fileNameDest)
            print(DEFAULT_PREFIX + "-- OK")


    ################
    # execute kgen #
    ################
    if runLAPW and not wCalc:
        #get our input parameters from the configuration file
        numKpoints = configFile['Number of K-Points']
        numKDivisions = configFile['K-Mesh Divisions']
        kShift = configFile['K-Mesh Shift']
        #input is a string per line representation of each input for that particular file
        inputListing = [str(numKpoints),
                        functools.reduce(lambda i,j: str(i) + ' ' + str(j), numKDivisions),
                        kShift,
                        ]
        shell_x_kgen = b_PySubProcess.VirtualShellInstance(
            'x', 'kgen', '-fbz',
            input=inputListing,
            ) # -fbz forces full BZ mesh
        try:
            shell_x_kgen.run()
        except subprocess.CalledProcessError as err:
            print(DEFAULT_PREFIX + "ERROR: in automation of " + str(structureName))
            print(DEFAULT_PREFIX + "ERROR --> " + shell_x_kgen.getCommandString())
            print(str(err))
            return False
    elif not runLAPW:
        print(DEFAULT_PREFIX + "Skip k-mesh generation, since you "+ \
            "opted for not performing LAPW.\n Proceed assuming that " + \
            "the proper k-mesh generated in the FULL BZ and vector "+ \
            "files were created in the previous run")
        print(DEFAULT_PREFIX + "k-mesh will be read from case.klist file")
        f = open( str(structureName)+'.klist' , 'r' )
        for line in f:
            kmesh = re.findall( r'\(.*?\)' , line )
            kmesh2 = re.findall( r'\d+' , kmesh[0] )
            numKDivisions = map(int, kmesh2)
            print(DEFAULT_PREFIX + "The following k-mesh will be used:", \
                numKDivisions)
            break
        f.close()


    ####################################
    # get length of k-path (-w option) #
    ####################################
    if wCalc:
        numKDivisions = [0]*3 # creale list [0,0,0]
        print(DEFAULT_PREFIX + "k-path will be read from case.klist file")
        f = open( str(structureName)+'.klist' , 'r' )
        word = ("END", "end")
        numKDivisions[0] = 0
        numKDivisions[1] = 1
        numKDivisions[2] = 1
        for line in f:
            if any(s in line for s in word): # check if end of list
                break
            numKDivisions[0] += 1
        f.close()
        print(DEFAULT_PREFIX + "The following k-path length will be used:", \
            numKDivisions)


    ##########################################
    # cp $filename.klist $filename.klist_w90 #
    ##########################################
    shell_pathphase_y_copy = b_PySubProcess.VirtualShellInstance(
        'cp',
        str(structureName)+'.klist',
        str(structureName)+'.klist_w90',
        )
    try:
        shell_pathphase_y_copy()
    except subprocess.CalledProcessError as err:
        print(DEFAULT_PREFIX + "ERROR: in automation of " + \
            str(structureName))
        print(DEFAULT_PREFIX + "ERROR --> " + \
            shell_pathphase_y_copy.getCommandString())
        print(str(err))
        return False


    #############################################################
    # Fake spin-polarized calculation if spin-orbit coupling is #
    # activated without spin polarization                       #
    # cp $filename.vsp $filename.vsp[up/dn]                     #
    # cp $filename.vns $filename.vns[up/dn]                     #
    #############################################################
    if soCalc and not spCalc:
        print(DEFAULT_PREFIX + "Faking spin-polarized calculation...")
        for ext in ['.vsp', '.vns']:
            for spin in ['up', 'dn']:
                shell_pathphase_y_copy = \
                    b_PySubProcess.VirtualShellInstance(
                    'cp',
                    str(structureName) + ext,
                    str(structureName) + ext + spin,
                    )
                try:
                    shell_pathphase_y_copy()
                except subprocess.CalledProcessError as err:
                    print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                        str(structureName))
                    print(DEFAULT_PREFIX + "ERROR --> " + \
                        shell_pathphase_y_copy.getCommandString())
                    print(str(err))
                    return False


    ##################
    # Parallel mode? #
    ##################
    if pCalc:
        pOpt = '-p'
    else:
        pOpt = ''

    ##############################
    # run lapw1 (-up/-dn) (-orb) #
    ##############################
    if runLAPW:
        rmerror('lapw1') # rm error files
        if spCalc or soCalc or orbCalc: # spin polarized mode
            # SOC calculation needs a faked SP (for w2w reasons)
            spOptions = ['-up', '-dn']
            if sp_cCalc: # constrained (non-magneric SP up=dn)
                spOptions = ['-up']
        else: # regular (non-spin polarized calc.)
            spOptions = ' ' # need space here!
        for suffix in spOptions:
            suffix = suffix + ' ' + pOpt # take care of parallel mode
            if orbCalc and not(soCalc): # use [-orb] for LAPW1 without SOC only
                suffix = suffix + ' ' + '-orb' # activate OP switch

            shell_x_lapw1 = b_PySubProcess.VirtualShellInstance(
                'x', 'lapw1', suffix
                )

            try:
                shell_x_lapw1()
            except subprocess.CalledProcessError as err:
                print(DEFAULT_PREFIX + "ERROR: in automation of " + str(structureName))
                print(DEFAULT_PREFIX + "ERROR --> " + shell_x_lapw1.getCommandString())
                print(str(err))
                return False
        testerror('lapw1') # test for non-empty error files
    else:
        print(DEFAULT_PREFIX + "Skipping 'x lapw1'")


    ###########################
    # run lapwso (-up) (-orb) #
    ###########################
    if soCalc and runLAPW: # spin-orbit coupling ONLY
        rmerror('lapwso') # rm error files
        suffix = '-up' # SOC calculation needs a faked SP (for w2w reasons)
        suffix = suffix + ' ' + pOpt # take care of parallel mode
        if orbCalc: # OP option
            suffix = suffix + ' ' + '-orb'

        shell_x_lapwso = b_PySubProcess.VirtualShellInstance(
            'x', 'lapwso', suffix
            )

        try:
            shell_x_lapwso()
        except subprocess.CalledProcessError as err:
            print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                str(structureName))
            print(DEFAULT_PREFIX + "ERROR --> " + \
                shell_x_lapwso.getCommandString())
            print(str(err))
            return False
        testerror('lapwso') # test for non-empty error files

    #############################################
    # run lapw2 -fermi -in1orig (-so) (-up/-dn) #
    #############################################
    # the purpose is to generate case.scf2 file that contains band occupancies
    # (-orb) is not required here as per Laurence Marks suggestion (Sep 2020)
    # (-in1orig) is added to prevent LAPW2 from midifying case.in1(c) file.
    #            Otherwise the electronic polarization results will be slightly
    #            different for consequitive runs of BerryPI as per Laurence Marks
    #            and Peter Blaha suggestion (Sep 2020)
    if not wCalc and bands == []:
        # do not run LAPW2 with pre-set band range
        # also not in the case of Weyl point analysys/Wilson loop
        rmerror('lapw2') # rm error files
        if spCalc or soCalc or orbCalc: # spin polarized mode
            # SOC calculation needs a faked SP (for w2w reasons)
            spOptions = ['-up', '-dn']
            if sp_cCalc: # constrained (non-magneric SP up=dn)
                spOptions = ['-up']
                # cp case.energyup* case.energydn* (needed for LAPW2 -up)
                # for f in lambda1.energyup*; do cp "$f" "$(sed 's/\(.*\)up/\1dn/' <<< "$f")"; done
                shell_copy = \
                    b_PySubProcess.VirtualShellInstance(
                    "for f in",
                    str(structureName) + '.energyup*;',
                    'do cp "' + "$f" + '" "$(sed ' + "'s/\(.*\)up/\\1dn/' " +\
                         '<<< "$f")"; done',
                    )
                try:
                    shell_copy()
                except subprocess.CalledProcessError as err:
                    print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                        str(structureName))
                    print(DEFAULT_PREFIX + "ERROR --> " + \
                        shell_copy.getCommandString())
                    print(str(err))
                    return False
        else: # regular (non-spin polarized calc. or SOC)
            spOptions = ' ' # need space here!
        for suffix in spOptions:
            suffix = suffix + ' ' + pOpt # take care of parallel mode
            if soCalc: # SOC
                suffix = suffix + ' ' + '-so'
            shell_x_lapw1 = b_PySubProcess.VirtualShellInstance(
                'x', 'lapw2 -fermi -in1orig', suffix
                )

            try:
                shell_x_lapw1()
            except subprocess.CalledProcessError as err:
                print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                    str(structureName))
                print(DEFAULT_PREFIX + "ERROR --> " + \
                    shell_x_lapw1.getCommandString())
                print(str(err))
                return False
        testerror('lapw2') # test for non-empty error files


    #############################################
    # Calculation to get the occupied blochBand #
    #############################################
    if bands == []: # band range is not specified
        if spCalc or soCalc or orbCalc: # spin polarized mode
            print(DEFAULT_PREFIX + "Determine number of bloch bands " + \
                "in spin-polarized mode based on *.scf2(up/dn)")
            spinList = ['up', 'dn']
            if sp_cCalc: # constrained (non-magneric SP up=dn)
                # it is not a typo (there is no case.scfdn file in this case)
                spinList = ['up', 'up']
            blochBandSettings = [] # allocate array for the range of bands
            for spin in spinList: # loop over spins
                print(DEFAULT_PREFIX + " "*2 + "spin = " + spin)
                try:
                    blochBandCalculation = \
                        b_PyCalc.CalculateNumberOfBands(structureName + \
                        '.scf2' + spin)
                except b_PyError.ParseError as err:
                    print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                        str(structureName))
                    print(DEFAULT_PREFIX + str(err))
                    print(DEFAULT_PREFIX + "ERROR: missing tags: " + \
                        str(err.errorTags))
                    return False
                blochBandRange = [ 1 , \
                    int(blochBandCalculation.getNumberOfBands(spCalc,soCalc, \
                    orbCalc,wCalc)) ]
                print(DEFAULT_PREFIX + " "*4 + \
                    "Number of bloch bands is {}".format(blochBandRange))
                blochBandSettings.append(blochBandRange)
        else: # regular (non-spin polarized calc.)
            print(DEFAULT_PREFIX + \
                "Determine number of bloch bands in non-sp mode")
            try:
                blochBandCalculation = \
                    b_PyCalc.CalculateNumberOfBands(structureName + '.scf2')
                # it is better to take *.scf2 file instead of
            except b_PyError.ParseError as err:
                print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                    str(structureName))
                print(DEFAULT_PREFIX + str(err))
                print(DEFAULT_PREFIX + "ERROR: missing tags: " + \
                    str(err.errorTags))
                return False
            blochBandSettings = [[ 1, \
                int(blochBandCalculation.getNumberOfBands(spCalc,soCalc, \
                orbCalc,wCalc)) ]]
            print(DEFAULT_PREFIX + " "*4 + \
                "Number of bloch bands is {}".format(blochBandSettings))
    else: # the band range is specified
        blochBandSettings = []
        blochBandSettings.append(configFile['Bloch Band'])
        if spCalc or soCalc or orbCalc: # spin polarized mode
            # append the same band range again for the 2nd spin chanel
            blochBandSettings.append(configFile['Bloch Band'])

   
    ###########################################
    # allocate an array for storing the pases #
    ###########################################
    if spCalc and not soCalc: # inludes -sp; -sp -orb, but not -so, -sp -so
        nkpt = numKDivisions[1]*numKDivisions[2]
        phaseX = numpy.zeros((2,nkpt,2)) # spin (up/dn), kpoint-path-start, phase
        nkpt = numKDivisions[0]*numKDivisions[2]
        phaseY = numpy.zeros((2,nkpt,2))
        nkpt = numKDivisions[0]*numKDivisions[1]
        phaseZ = numpy.zeros((2,nkpt,2))
    else: # non-spin polarized calc. or -so or -sp -so
        nkpt = numKDivisions[1]*numKDivisions[2]
        phaseX = numpy.zeros((1,nkpt,2)) # spin degenerate(1), kpoint-path-start, phase
        nkpt = numKDivisions[0]*numKDivisions[2]
        phaseY = numpy.zeros((1,nkpt,2))
        nkpt = numKDivisions[0]*numKDivisions[1]
        phaseZ = numpy.zeros((1,nkpt,2))

    ##############################################
    # THIS IS THE MAIN LOOP OVER SPINS AND BANDS #
    ##############################################
    # Prepare spin lables even for non-sp calculation
    if spCalc or soCalc or orbCalc: # spin-polar. calculation or -so
        spinIndexes = [0, 1] # 0 - up, 1 - dn
        spinLables = ["up", "dn"]
        spinOptions = ["-up", "-dn"]
        spOption = "-sp"
        if sp_cCalc: # constrained (non-magneric SP up=dn)
            spinIndexes = [0] # 0 - up (dn=up)
            spinLables = ["up"]
            spinOptions = ["-up"]
            spOption = "-sp"
    else:# regular (non-spin polarized calc.)
        spinIndexes = [0]
        spinLables = [""] # use empty; brackets [] are for zip()
        spinOptions = [""]
        spOption = ""
    args = zip(spinIndexes, spinLables, spinOptions)
    for spinIndex, spinLable, spinOption in args: # LOOP OVER SPINS

        ####################################
        # create input list for write_inwf #
        ####################################
        if bands != []: # range of bands is specified (typicaly for Wilson loop)
            inputListing = '-mode MMN -bands ' + str(bands[0]) + ' ' +\
                 str(bands[1])
        else: # no range of bands is specified (default)
            nb = blochBandSettings[spinIndex][1] # last occupied band
            inputListing = '-mode MMN -bands 1 ' + str(nb)


        ######################
        # execute write_inwf #
        ######################
        exe = 'write_inwf'
        shell_write_inwf = b_PySubProcess.VirtualShellInstance(
            pypath,
            w2kpath + '/' + exe,
            inputListing,
            )
        try:
            shell_write_inwf.run()
        except subprocess.CalledProcessError as err:
            print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                str(structureName))
            print(DEFAULT_PREFIX + "ERROR --> " + \
                shell_write_inwf.getCommandString())
            print(str(err))
            return False

        ##################################
        # execute write_win -band nofile #
        ##################################
        exe = 'write_win -band nofile'
        shell_write_win = b_PySubProcess.VirtualShellInstance(
            exe,
            )
        case = str(structureName)
        fnameold = case + '.win'
        if os.path.isfile(fnameold): # remove old case.win if present
            print(DEFAULT_PREFIX, fnameold, 'is present and will be removed')
            os.remove(fnameold)
            if not os.path.isfile(fnameold):
                print('             ... done')
            else:
                print(DEFAULT_PREFIX, 'unable to remove', fnameold)
                print('             will continue at your own risk, w2w may crash')
        try:
            shell_write_win()
        except subprocess.CalledProcessError as err:
            print(DEFAULT_PREFIX + "ERROR: in automation of " + str(structureName))
            print(DEFAULT_PREFIX + "ERROR --> " + shell_write_win.getCommandString())
            print(str(err))
            return False
        if not spinLable == "": # for spin-polarized calc.
            fnamenew = fnameold + spinLable
            shutil.copy(fnameold, fnamenew) # copy case.win to case.win[up/dn]

        ######################################
        # move case.inwf -> case.inwf[up/dn] #
        ######################################
        if not spinLable == "": # for spin-polarized calc.
            case = str(structureName)
            fnameold = case + '.inwf'
            fnamenew = fnameold + spinLable
            print(DEFAULT_PREFIX + "Move " + fnameold + " into " + fnamenew)
            shutil.move(fnameold, fnamenew) # move case.inwf to case.inwf[up/dn] 


        #######################
        # execute win2nnkp.py #
        #######################
        if wCalc:
            wOption = '-w' + ' ' + str(wTranslDir)
        else:
            wOption = ''
        shell_win2nnkp = b_PySubProcess.VirtualShellInstance(
                pypath,
                bppath + '/' + 'win2nnkp.py',
                str(structureName), wOption,
                )
        try:
            shell_win2nnkp()
        except subprocess.CalledProcessError as err:
            print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                  str(structureName) + suffix)
            print(DEFAULT_PREFIX + "ERROR --> " + \
                  shell_win2nnkp.getCommandString())
            print(str(err))
            return False


        #####################
        # create case.fermi #
        #####################
        case = str(structureName)
        fname = case + '.fermi' + spinLable # write Fermi energy to case.fermi file
        if os.path.isfile(fname): # remove old case.fermi if present
            print(DEFAULT_PREFIX, fname, 'is present and will be removed')
            os.remove(fname)
            if not os.path.isfile(fname):
                print('             ... done')
            else:
                print(DEFAULT_PREFIX, 'unable to remove', fname)
                print('             will continue at your own risk')
        file = open(fname, "w")
        # write Ef = 0. The value makes no difference for w2w (only file with 
        # a value should be present). The value is important only for a Wannier
        # interpolated band structure plot.
        file.write("0.0" + '\n') # add end of line
        file.close()
        print(DEFAULT_PREFIX + 'Fermi energy file', fname, 'created')


        ###############
        # execute w2w #
        ###############
        suffix = pOpt # take care of parallel mode
        rmerror('w2w') # rm error files
        if soCalc: # for spin-orbit call executables directly
            shell_w2w = b_PySubProcess.VirtualShellInstance( \
                'x w2w -so', spinOption, suffix )
        else:
            shell_w2w = b_PySubProcess.VirtualShellInstance( \
                'x w2w', spinOption, suffix )
        try:
            shell_w2w()
        except subprocess.CalledProcessError as err:
            print(DEFAULT_PREFIX + "ERROR: in automation of " + str(structureName))
            print(DEFAULT_PREFIX + "ERROR --> " + shell_w2w.getCommandString())
            print(str(err))
            return False
        testerror('w2w') # test for non-empty error files


        ####################
        # execute w2waddsp #
        ####################
        if soCalc and spinLable=='dn': # spin-orbit last passage
            shell_csf = b_PySubProcess.VirtualShellInstance( \
                'x w2waddsp' )
            try:
                shell_csf()
            except subprocess.CalledProcessError as err:
                print(DEFAULT_PREFIX + "ERROR: in automation of " + str(structureName))
                print(DEFAULT_PREFIX + "ERROR --> " + shell_csf.getCommandString())
                print(str(err))
                return False


        ######################################################
        # Compute Berry phase:                               #
        # python ./mmn2pathphase.py case x/y/z (-up/-dn/-so) #
        ######################################################
        if soCalc and spinLable=='up':
            pass # spin-orbit skip first passage
        else:
            if soCalc:
                spinIndex = 0 # reset counter for proper writing below
            directions = ['x', 'y', 'z']
            for direction in directions:
                if not soCalc:
                    opt = spinOption
                else:
                    opt = '-so' # work around for spin-orbit
                if wCalc:
                    opt = opt + ' ' + '-w'
                args = [str(structureName), direction, opt]
                pathPhases = mmn2pathphase.main(args)
                if wCalc:
                    break
                if   direction == 'x':
                    phaseX[spinIndex,:,:] = pathPhases
                elif direction == 'y':
                    phaseY[spinIndex,:,:] = pathPhases
                elif direction == 'z':
                    phaseZ[spinIndex,:,:] = pathPhases
                if sp_cCalc: # constrained non-magnetic SP (up=dn)
                    # copy results up=dn
                    phaseX[spinIndex+1,:,:] = phaseX[spinIndex,:,:]
                    phaseY[spinIndex+1,:,:] = phaseY[spinIndex,:,:]
                    phaseZ[spinIndex+1,:,:] = phaseZ[spinIndex,:,:]

    # END OF THE MAIN LOOP FOR BERRY PHASE CALCULATION+++++++++++++++++++
    
    results = [phaseX,phaseY,phaseZ]

    ####################################################################
    # Clean up files intended to fake SP calculation if SO coupling is #
    # activated without spin polarization                              #
    # rm $filename.vsp[up/dn]                                          #
    # rm $filename.vns[up/dn]                                          #
    #                                                                  #
    # In case of constrained spin-polarized calculation                #
    # rm case.energydn*                                                #
    ####################################################################
    if soCalc and not spCalc:
        print(DEFAULT_PREFIX + "Removing duplicated files...")
        for ext in ['.vsp', '.vns']:
            for spin in ['up', 'dn']:
                shell_pathphase_y_rm = \
                    b_PySubProcess.VirtualShellInstance(
                    'rm',
                    str(structureName) + ext + spin,
                    )
                try:
                    shell_pathphase_y_rm()
                except subprocess.CalledProcessError as err:
                    print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                        str(structureName))
                    print(DEFAULT_PREFIX + "ERROR --> " + \
                        shell_pathphase_y_rm.getCommandString())
                    print(str(err))
                    return False
    elif sp_cCalc: # rm case.energydn (created for LAPW2 -up)
        print(DEFAULT_PREFIX + "Removing duplicated files...")
        shell_copy = \
            b_PySubProcess.VirtualShellInstance(
            'rm',
            str(structureName) + '.energydn*',
            )
        try:
            shell_copy()
        except subprocess.CalledProcessError as err:
            print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                str(structureName))
            print(DEFAULT_PREFIX + "ERROR --> " + \
                shell_copy.getCommandString())
            print(str(err))
            return False


    ###########################################
    # Restore klist file                      #
    # mv $filename.klist_orig $filename.klist #
    ###########################################
    if 'resoreKlist' in locals(): # check if the logical variable was set
        if resoreKlist:
            fileNameSrc = structureName + '.klist_orig'
            fileNameDest = structureName + '.klist'
            print(DEFAULT_PREFIX + "Trying to move " + fileNameSrc + " to " + fileNameDest)
            if os.path.isfile(fileNameSrc): # check if $filename.klist_orig is in place
                shutil.move(fileNameSrc, fileNameDest)
                print(DEFAULT_PREFIX + "-- OK")
            else:
                msg = "file " + fileNameSrc + " does not exist. Move aborted."
                print(DEFAULT_PREFIX + msg)
    

    print(DEFAULT_PREFIX + "Finished Berry phase computation for " + \
        structureName)

    return [True,results]
# END rawBerryPhase


if __name__ == "__main__":
    # Set up parser for line arguments
    parser = argparse.ArgumentParser()
    parser.add_argument("-k",\
        help="k mesh for sampling the BZ (default 4 4 4)",\
        nargs=3,\
        type=int,\
        default=[4, 4, 4])
    parser.add_argument("-sp",\
        help="spin polarized calculation",\
        action="store_true")
    parser.add_argument("-sp_c",\
        help="constrained spin polarization (no net magnetic moment, up=dn)",\
        action="store_true")
    parser.add_argument("-orb",\
        help="calculation with an additional orbit potential (e.g., LSDA+U)",\
        action="store_true")
    parser.add_argument("-so",\
        help="enable spin-orbit coupling",\
        action="store_true")
    parser.add_argument("--skip-lapw",\
        help="skip lapw1 (and lapwso if -so is active) "+\
            "and proceed directly to w2w",\
        action="store_true")
    parser.add_argument("-p",\
        help='run "x lapw1 -p", "x lapwso -p" and "x w2w -p" in parallel '+\
            'mode (needs .machines file)',\
        action="store_true")
    parser.add_argument("-w",\
        help="compute Berry phase along a specific (closed loop) k-path given "+\
            "in the case.klist file (used for topological Weyl semimetals "+\
            "and Wannier charge centers). "+\
            "First k(1) and last k(n) points in the case.klist file will be "+\
            "joined to form a closed loop k(1) -> k(2) -> ... -> k(n) -> k(1). "+\
            "The argument (optional) specifies when a periodic image of k(1) "+\
            "should be used. For example '-w 2' means the following path: "+\
            "k(1) -> k(2) -> ... -> k(n) -> k(1)+G[2], where G[2] is the second "+\
            "reciprocal lattice vector. By default we assume '-w 0', which implies "+\
            "that no translation is added at the end of the loop.",\
        nargs='?',\
        const=0,\
        choices=[0, 1, 2, 3],\
        type=int)
    parser.add_argument("-b",\
        help="range of Bloch bands used to construct projections "+\
            "(used for Weil semimetals only in conjunction with -w option)",\
        nargs=2,\
        type=int)
    parser.add_argument("-v",\
        help="increase output verbosity",\
        action="store_true")
    args = parser.parse_args()
    if args.v:
        print("verbosity turned on")
        print("args=", args)

    #boolean for whether to check config for
    #configuration dicts
    checkConfig = True

    #boolean for whether to automate the entire process
    checkAuto = True

    # Assign line arguments parsed by "argparse"
    structure_kmesh = args.k # k mesh
    bands = args.b # range of bands
    runLAPW = not(args.skip_lapw) # skip LAPW
    VERBOSE = args.v # verbose output
    pCalc = args.p   # calculation is serial by defaul
    spCalc = args.sp  # no spin polarization by default
    orbCalc = args.orb # no additional orbital potential by default
    soCalc = args.so  # no spin-orbit coupling by default
    wTranslDir = args.w # no k-path specified by default
    sp_cCalc = args.sp_c # no constrained spin-polarized by default
    
    # Set default case name and path based on where berrypi is executed
    structure_path = os.getcwd()
    structure_name = os.path.split(os.getcwd())[-1]
    
    # redirect stdout to case.outputberry file and console
    # case mane is taken from structure_name
    sys.stdout = Logger()

    #echo command line arguments and switch implied options
    if spCalc:
        print(DEFAULT_PREFIX+'Spin polarization is activated')
    if sp_cCalc:
        print(DEFAULT_PREFIX+'Constrained (non-magnetic) ' +\
            'spin polarization is activated')
        spCalc = True
    if orbCalc:
        print(DEFAULT_PREFIX+'Calculation with an additional '+ \
            'orbital potential is activated')
        print(DEFAULT_PREFIX+'Spin polarization is activated '+ \
            'automatically with adding orbital potential')
        spCalc = True
    if soCalc:
        print(DEFAULT_PREFIX+'Calculation with spin-orbit '+ \
            'coupling is activated')
    if not(runLAPW):
        print(DEFAULT_PREFIX+'Calculations will skip LAPW')
    if pCalc:
        print(DEFAULT_PREFIX+"Parallel calculation is activated")
        print("            it is assumed that you have .machines file ready")
    if bands == None:
        print(DEFAULT_PREFIX+'Range of bands [1, last occupied] will be '+\
            'determined automatically')
    else:
        print(DEFAULT_PREFIX+"Only selected range of bands will be used ")
        print(DEFAULT_PREFIX, "bands=", bands)
    if wTranslDir == None: # -w option not set
        wCalc = False
    else:
        wCalc = True
        print(DEFAULT_PREFIX+"Enable Berry phase calculation along a "+\
            "specific k-path given in the case.klist file. ")
        print(DEFAULT_PREFIX+"Any input of the k-mesh (-k ... option) "+\
            "will be ignored")
        if wTranslDir == 0:
            print(DEFAULT_PREFIX+"The first and last k point in "+\
                "the case.klist file "+\
                "will be connected to form a closed loop")
        else:
            print(DEFAULT_PREFIX+"The first and last k point in "+\
                "the case.klist file "+\
                f"will be connected by adding G[{wTranslDir}] reciprocal "+\
                "lattice vector to form a closed loop.")
        if not(soCalc):
            soCalc = True
            print(DEFAULT_PREFIX+'Spin-orbit coupling is activated '+ \
                'automatically with -w option')
        

    # Check input consistency
    if len(structure_kmesh)!=3:
        #Makes sure that the k-mesh length of the tuple is 3
        print(f'You entered: {structure_kmesh}')
        raise ValueError('Please check that you have 3 integers in the kmesh')
    elif 0 in structure_kmesh:
        #Makes sure that there are no values of 0 in the kmesh
        print(f'You entered: {structure_kmesh}')
        raise ValueError('Cannot have a dimension of 0 in the kmesh')
    elif min(structure_kmesh)<0:
        #Sets the lower bound of the dimensions
        print(f'You entered: {structure_kmesh}')
        raise ValueError('You cannot have negative integers as kmesh values')
    elif wCalc and bands == None:
        #Band range is not specified for Wilson loop calculation
        raise ValueError('Wilson loop calculation requires the band range '
            'to be specified. Check "berrypi -h"')

    if bands != None and len(bands)!=2:
        #Makes sure that the band range length of the tuple is 2
        print(f'You entered: {bands}')
        raise ValueError('Please check that you have 2 integers in the band range')
    elif bands != None and 0 in bands:
        #Makes sure that there are no values of 0 in the bands
        print(f'You entered: {bands}')
        raise ValueError('Cannot have a dimension of 0 in the bands range')
    elif bands != None and min(bands)<0:
        #Sets the lower bound of the dimensions
        print(f'You entered: {bands}')
        raise ValueError('You cannot have negative integers as bands range values')
    elif bands != None and bands[0] > bands[1]:
        #The range of bands should start with a smaller value
        print(f'You entered: {bands}')
        raise ValueError('The range of bands should start with a smaller value')
            
    #Determines whether or not W2WANNIER is installed
    wienroot = os.getenv('WIENROOT') # get path to Wien2k installation
    try:
        with open(wienroot + '/w2w'): pass
    except IOError:
        print('Please make sure that W2W is installed before trying to run BerryPI')
        
    #if configuration values are provided, create config file to pass
    #to the automation and calculations
    if structure_name or structure_path:
        if structure_name and structure_path:
            #follows same format as config.py configurations
            configFile = {
                'Structure Name' : structure_name,
                'Structure Path' : structure_path,
                'Number of K-Points' : config.DEFAULT_NUMBER_OF_KPOINTS,
                'K-Mesh Divisions' : structure_kmesh,
                'K-Mesh Shift' : config.DEFAULT_KMESH_SHIFT, 
                'Bloch Band' : bands,
                'Number of Wannier Functions' : config.DEFAULT_WANNIER_FUNCTIONS,
                'Center Atom and Character' : config.DEFAULT_CENTER_ATOM_AND_CHARACTER,
                'Perform LAPW' : runLAPW,
                'Parallel mode' : pCalc,
                'Spin polarized calculation' : spCalc,
                'Add orbital potential' : orbCalc,
                'Spin-orbit coupling' : soCalc,
                }
            checkConfig = False
        else:
            raise RuntimeError("structure name and structure path must be provided")

    #check config (based on if structure name or
    #structure path has been passed to berryPyRunAutomation.py)
    bAutomationErrorChk = True
    if checkConfig:
        ConfigurationFiles = config.berryPyConfigAutomationList
        for berryConfig in ConfigurationFiles:
            #check whether to automate
            if checkAuto:
                bAutomationErrorChk = rawBerryPhase(berryConfig)
            
            if bAutomationErrorChk and not wCalc:
                try:
                    theCalc = postProcBerryPhase(berryConfig)
                except b_PyError.ParseError as err:
                    print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                        berryConfig['Structure Name'])
                    print(DEFAULT_PREFIX + str(err))
                    print(DEFAULT_PREFIX + "ERROR: missing tags: " + \
                        str(err.errorTags))
                    continue
                
                raise RuntimeError('There is an inssue with config variables. '
                        'Check the code after "if checkConfig:"')
    else:
        #assuming the configFile was passed correctly
        if checkAuto:
            [bAutomationErrorChk,phasesRaw] = rawBerryPhase(configFile)
            
        if bAutomationErrorChk and not wCalc:
            try:
                theCalc = postProcBerryPhase(configFile,phasesRaw,spCalc)
            except b_PyError.ParseError as err:
                print(DEFAULT_PREFIX + "ERROR: in automation of " + \
                    configFile['Structure Name'])
                print(DEFAULT_PREFIX + str(err))
                print(DEFAULT_PREFIX + "ERROR: missing tags: " + \
                    str(err.errorTags))
    
    
    # Final remarks
    w2kpath = os.getenv('WIENROOT') # get path to WIEN2k
    # open version file
    with open(w2kpath+'/SRC_BerryPI/BerryPI/version.info', 'r') as versionFile:
        print("\nCompleted using BerryPI version:", versionFile.read())
    print('''
Suggested reference:
[1] S.J.Ahmed, J.Kivinen, B.Zaporzan, L.Curiel, S.Pichardo and O.Rubel
    "BerryPI: A software for studying polarization of crystalline solids with 
    WIEN2k density functional all-electron package"
    Comp. Phys. Commun. 184, 647 (2013)
    https://doi.org/10.1016/j.cpc.2012.10.028

Questions and comments are to be communicated via the WIEN2k mailing list
(see http://susi.theochem.tuwien.ac.at/reg_user/mailing_list)''')