scaleLJxml

#!/usr/bin/env python
# scaleLJxml - scale epsilon and sigma LJ parameters in pair.lmp LAMMPS file.
# Agilio Padua <agilio.padua@ens-lyon.fr>
# Kateryna Goloviznina <kateryna.goloviznina@ens-lyon.fr>
# version 2023/11/24

import sys
import argparse
import os
import xml.etree.ElementTree as ET
#from lxml import etree as ET

usage = """
==============================================================================
Scale epsilon and sigma LJ parameters in pair.lmp LAMMPS file
------------------------------------------------------------------------------
Format of file containing specification of Drude oscillators
and polarisability values (alpha.ff):
  # type  dm/u   dq/e  k/(kJ/molA2)  alpha/A3  a_thole
  CT      0.4    -1.0     4184.0      1.016     2.6
  ...
* dm is the mass to place on the Drude particle (taken from its core),
* dq is the charge to place on the Drude particle (taken from its core),
* k is the harmonic force constant of the bond between core and Drude,
* alpha is the polarizability, hyrdogen aroms are not merged,
* a_thole is a parameter of the Thole damping function.
------------------------------------------------------------------------------
Format of file containing monomers and dimers specification (fragment.ff)
  MONOMERS
  # name       q/e       mu/D
  c2c1im       1.0      1.1558
  ..
  DIMERS
  # m1         m2       r_COM/A    k_sapt
  c2c1im       dca       2.935      0.61
  ...
* q is the charge of the monomer,
* mu is the dipole moment of the monomer,
* m1 and m2 are the monomers forming a dimer,
* r_COM distance between the center of mass of the monomers,
* k_sapt is the scaling factor for the epsilon of LJ potential,
obtained by SAPT quantum calculation (optional).
------------------------------------------------------------------------------
Format of file containing fragments list with atomic names (frag.inp)
  #name  atoms
  c2c1im   NA CR CW C1 HCR C1A HCW H1
  c4h10    C2 CS HC CT
  ntf2     CBT F1 SBT NBT OBT
  ...
* atomic names associated with this fragment in field-p.lmp file.
------------------------------------------------------------------------------
Script requires structure files of monomers used in frag.inp file
(.zmat, .xyz, .mol, .pdb)
==============================================================================
"""

class _Const(object):

    sigma_k = 0.985

    # C0 and C1 are coefficients to predict k
    @staticmethod
    def C0():
        return 0.254952

    @staticmethod
    def C1():
        return 0.106906

    @staticmethod
    def isfloat(value):
        try:
            float(value)
            return True
        except ValueError:
            return False

class FieldXML(object):
    '''force field from OpenMM xml file'''

    def __init__(self, infile):
        '''element tree from xml file'''

        self.ftree = ET.parse(infile)
        self.root = self.ftree.getroot()

    def Write(self, outfile):
        '''write force field to xml file'''

        ET.indent(self.root, space=' ')
        self.ftree.write(outfile)

    def Scale(self, ff, syst, sapt, scsig):
        try:
            ## i-i interactions
            atom_lj = self.root.find('LennardJonesForce').findall('Atom')
            for ii_pair in atom_lj:

                i = ii_pair.get('class')
                eps = float(ii_pair.get('epsilon'))
                sig = float(ii_pair.get('sigma'))

                for frp in syst.fragmentpairs:
                    if not sapt:
                        k = frp.k_pred
                    elif frp.k_sapt is not None:
                        k = frp.k_sapt
                    else:
                        raise Exception(f'  error: k_sapt for {frp.fr1.name} {frp.fr2.name} dimer not found in {ff.filename}')
                    if (i in frp.fr1.ind_range and i in frp.fr2.ind_range) and k < 1:
                        eps *= k
                        break

                if (scsig is not None):
                    fi = next((x for x in syst.fragments if (i in x.ind_range)), None)

                    if (fi is not None):
                        if fi.scale_sig:
                            sig *= _Const.sigma_k

                ii_pair.set('epsilon', f'{eps:.5f}')
                ii_pair.set('sigma', f'{sig:.4f}')

            ## i-j interactions
            nb_fix_pair = self.root.find('LennardJonesForce').findall('NBFixPair')

            for ij_pair in nb_fix_pair:
                i = ij_pair.get('class1')
                j = ij_pair.get('class2')
                eps = float(ij_pair.get('epsilon'))
                sig = float(ij_pair.get('sigma'))

                for frp in syst.fragmentpairs:
                    if not sapt:
                        k = frp.k_pred
                    elif frp.k_sapt is not None:
                        k = frp.k_sapt
                    else:
                        raise Exception(f'  error: k_sapt for {frp.fr1.name} {frp.fr2.name} dimer not found in {ff.filename}')
                    if ((i in frp.fr1.ind_range and j in frp.fr2.ind_range) or (j in frp.fr1.ind_range and i in frp.fr2.ind_range)) and k < 1:
                        eps *= k
                        break

                if (scsig is not None):
                    fi = next((x for x in syst.fragments if (i in x.ind_range)), None)
                    fj = next((x for x in syst.fragments if (j in x.ind_range)), None)

                    if (fi is not None and fj is not None):
                        if fi.scale_sig and fj.scale_sig:
                            sig *= _Const.sigma_k

                ij_pair.set('epsilon', f'{eps:.8f}')
                ij_pair.set('sigma', f'{sig:.8f}')

        except Exception as e:
            print(e)
            sys.exit(1)

    def CheckAtoms(self,syst):
        atoms_types_inp = [y for x in syst.fragments for y in x.ind_range]
        atom_types_xml = [x.get('name').split('-')[1] for x in self.root.find('AtomTypes').findall('Type') if x.get('name')[0] != 'D']
        extra_atoms = list(set(atom_types_xml) - set(atoms_types_inp))
        if len(extra_atoms) > 0:
            extra_atoms = ', '.join(extra_atoms)
            print(f'warning: atom(s) {extra_atoms} do/does not belong to any fragment in {syst.filename}')

class Monomer(object):
    '''monomer in fragment.ff'''

    def __init__(self, name, q, mu):
        self.name = name
        self.q = q
        self.mu = mu

    def __str__(self):
        return f'{self.name:10s}  q = {self.q:2d}   mu = {self.mu:6.4f} D'

class Dimer(object):
    ''' dimer in fragment.ff'''

    def __init__(self, m1, m2):
        self.m1 = m1
        self.m2 = m2
        self.r = None
        self.k_sapt = None

    def SetR(self, r):
        self.r = r

    def SetKSAPT(self,  k_sapt):
        self.k_sapt = k_sapt

    def __str__(self):
        res = f'{self.m1.name:10s} {self.m2.name:10s}'
        if self.r is not None:
            res += f'    r = {self.r:6.4f}'
        if self.k_sapt is not None:
            res += f'    k_sapt = {self.k_sapt:4.2f}'
        return res

class FragmentFField(object):
    '''fragment force field from fragment.ff'''

    def __init__(self, filename):
        '''read fragment.ff'''

        self.filename = filename
        self.monomers = []
        self.dimers= []
        try:
            with open(self.filename, 'r') as f:
                for line in f:
                    if line.startswith('#') or line.strip() == '':
                        continue
                    if line.lower().startswith('monomer'):
                        section = 'monomers'
                        continue
                    elif line.lower().startswith('dimer'):
                        section = 'dimers'
                        continue

                    tok = line.strip().split()

                    if section == 'monomers':
                        if tok[0].endswith('+') or tok[0].endswith('-'):
                            tok[0] = tok[0][:-1]
                        name = tok[0].lower()
                        q = float(tok[1])
                        mu  = float (tok[2])

                        if next((x for x in self.monomers if x.name == name), None) is not None:
                            raise Exception(f'  error: monomer {name} is specified twice in {self.filename}')

                        self.monomers.append(Monomer(name, q, mu))

                    elif section == 'dimers':
                        for i in range(0,2):
                            if tok[i].endswith('+') or tok[i].endswith('-'):
                                tok[i] = tok[i][:-1]
                        m1_name = tok[0].lower()
                        m2_name = tok[1].lower()
                        r = float (tok[2])

                        if next((x for x in self.dimers if (x.m1.name == m1_name and x.m2.name == m2_name) or (x.m1.name == m2_name and x.m2.name  == m1_name)), None):
                            raise Exception(f'  error: dimer {m1_name} {m2_name} is specified twice in {self.filename}')

                        d = self.SetDimer(m1_name, m2_name)
                        self.dimers.append(d)
                        d.SetR(r)
                        if len(tok)>3:
                            k_sapt = float (tok[3])
                            d.SetKSAPT(k_sapt)

        except IOError:
            print(f'  error: force field file {self.filename} not found')
            sys.exit(1)
        except IndexError:
            print(f'  error: incorrect force field file line: {line}')
            sys.exit(1)
        except Exception as e:
            print(e)
            sys.exit(1)

    def SetDimer(self, m1_name, m2_name):
        '''create dimer from two monomers'''

        try:
            m1 = next((x for x in self.monomers if x.name == m1_name), None)
            m2 = next((x for x in self.monomers if x.name == m2_name), None)
            if (m1 is not None and m2 is not None):
                d = Dimer(m1, m2)
                return d
            else:
                raise Exception(f'  error: monomer {m1_name} or {m2_name} not described in monomers section of {self.filename}')

        except Exception as e:
            print(e)
            sys.exit(1)

    def __str__(self):
        res = self.filename
        res += '\nMONOMERS'
        for m in self.monomers:
            res += '\n'+str(m)
        res += '\nDIMERS'
        for d in self.dimers:
            res += '\n'+str(d)
        return res

class Fragment(Monomer):
    """
    a fragment in frag.inp
    fragment is based on monomer with index range in current system added
    """

    def __init__(self, m, ind_range, res):
        self.name = m.name
        self.q = m.q
        self.mu = m.mu
        self.ind_range = ind_range
        self.pol_model = res[0]
        self.scale_eps = res[1]
        self.scale_sig = False

    def __str__(self):
            return '{0:10s}  q = {1:5.2f}  mu = {2:6.4f} D    Pol = {3:5s} ScaleEps = {4:5s} ScaleSig = {5:5s} Atoms = {5:s} '.format(self.name, self.q, self.mu,self.pol_model, self.scale_eps, self.scale_sig, self.ind_range)

    @staticmethod
    def PolExclude(m, p):
        '''check if fragment is polarisable'''

        pol_model = False
        scale_eps = True

        if p is not None:
            if m.name in p:
                pol_model = True
                if (m.q == 0 and m.mu == 0):
                    scale_eps = False

        return (pol_model,scale_eps)

    def GetAtomsFromMol(self, pol):
        '''read atoms from .mol for a given fragment'''

        try:
            filename = self.name + '.mol'
            with open(filename, 'r') as f:
                atoms = []
                tok = f.readline().strip().split()
                self.name = tok[0]
                line = f.readline()
                line = f.readline()
                tok = f.readline().strip().split()   # counts line
                natom = int(tok[0])

                for i in range(natom):
                    tok = f.readline().strip().split()
                    at_name = tok[3]
                    print(at_name)
                    atom = next((x for x in pol.atomtypes if x.name == at_name), None)
                    if atom is None:
                        raise Exception(f'  error: atom type {at_name} not found in {pol.filename}')
                    atoms.append(atom)
            self.atoms = atoms

        except IOError:
            print('  error: fragment mol file {filename} not found')
            sys.exit(1)
        except Exception as e:
            print(e)
            sys.exit(1)

    def GetAtomsFromPdb(self, pol):
        '''read atoms from .pdb for a given fragment'''

        try:
            filename = self.name + '.pdb'
            with open(filename, 'r') as f:
                atoms = []
                line = f.readline()
                while 'COMPND' not in line:
                    line = f.readline()
                tok = line.strip().split()
                self.name = tok[1]
                line = f.readline()
                while not (line.startswith('ATOM') or line.startswith('HETATM')):
                    line = f.readline()
                while 'ATOM' in line or 'HETATM' in line:
                    tok = line.strip().split()
                    at_name = tok[2]
                    atom = next((x for x in pol.atomtypes if x.name == at_name), None)
                    if atom is None:
                        raise Exception(f'  error: atom type {at_name} not found in {pol.filename}')
                    atoms.append(atom)
                    line = f.readline()
            self.atoms = atoms

        except IOError:
            print('  error: fragment pdb file {filename} not found')
            sys.exit(1)
        except Exception as e:
            print(e)
            sys.exit(1)


    def GetAtomsFromZmat(self, pol):
        '''read atoms from .zmat for a given fragment'''

        try:
            filename = self.name + '.zmat'
            with open(filename, 'r') as f:
                atoms = []
                line = f.readline()
                while line.strip().startswith('#'):
                    line = f.readline()
                fr_name = line.strip()
                line = f.readline()
                while line.strip().startswith('#') or line.strip() == '':
                    line = f.readline()

                tok = line.strip().split()
                if len(tok) > 1:   # there can be line numbers
                    shift = 1
                else:
                    shift = 0
                while line:
                    tok = line.strip().split()
                    if len(tok) == 0:
                        break
                    at_name = tok[shift]
                    atom = next((x for x in pol.atomtypes if x.name == at_name), None)
                    if atom is None:
                        raise Exception(f'  error: atom type {at_name} not found in {pol.filename}')
                    atoms.append(atom)
                    line = f.readline()
            self.atoms = atoms

        except IOError:
            print('  error: fragment zmat file {filename} not found')
            sys.exit(1)
        except Exception as e:
            print(e)
            sys.exit(1)

    def GetAtomsFromXYZ(self, pol):
        '''read atoms from .xyz for a given fragment'''

        try:
            filename = self.name+'.xyz'
            with open(filename, 'r') as f:
                atoms = []
                line = f.readline()
                while line.strip().startswith('#'):
                    line = f.readline()
                n_atoms = int(line.strip())
                line = f.readline()
                fr_name = line.strip()
                line = f.readline()
                while line.strip().startswith('#') or line.strip() == '':
                    line = f.readline()

                tok = line.strip().split()
                while line:
                    tok = line.strip().split()
                    if len(tok) == 0:
                        break
                    at_name = tok[0]
                    atom = next((x for x in pol.atomtypes if x.name == at_name), None)
                    if atom is None:
                        raise Exception(f'  error: atom type {at_name} not found in {pol.filename}')
                    atoms.append(atom)
                    line = f.readline()
            self.atoms = atoms

        except IOError:
            print(f'  error: fragment xyz file {filename} not found')
            sys.exit(1)
        except Exception as e:
            print(e)
            sys.exit(1)

class FragmentPair(Dimer):
    '''fragment pair based on dimer but using fragments (not monomers) as units'''

    def __init__(self, fr1, fr2, r, k_sapt):
        self.fr1 = fr1
        self.fr2  = fr2
        self.r = r
        self.k_sapt = k_sapt

    def __str__(self):
        res = f'{self.fr1.name:10s} {self.fr2.name:10s}'
        if self.r is not None:
            res += f'    r = {self.r:6.4f}'
        if self.k_sapt is not None:
            res += f'    k_sapt = {self.k_sapt:4.2f}'
        return res

    def PredictK(self):
        """
        predict k factor for a given fragment pair based on:
        total charge, alpha, dipole moment of fragment and interfragment distance
        """

        alpha_fr1 = [sum(x.alpha for x in self.fr1.atoms)][0]
        alpha_fr2 = [sum(x.alpha for x in self.fr2.atoms)][0]
        k_pred = 1.0

        if not self.fr1.pol_model:
            k_pred += _Const.C0()*self.r*self.r*(self.fr2.q*self.fr2.q)/alpha_fr2
            k_pred += _Const.C1()*(self.fr2.mu*self.fr2.mu)/alpha_fr2

        if not self.fr2.pol_model:
            k_pred += _Const.C0()*self.r*self.r*(self.fr1.q*self.fr1.q)/alpha_fr1
            k_pred += _Const.C1()*(self.fr1.mu*self.fr1.mu)/alpha_fr1

        k_pred = 1/k_pred
        self.k_pred = k_pred

class AtomType(object):
    '''atom type with polarisability value'''

    def __init__(self, name, alpha):
        self.name = name
        self.alpha = alpha

    def __str__(self):
        return f'{self.name:s}    {self.alpha:6.3f}'

class Polarisation(object):
    '''polarisability values for all atom types; from alpha.ff'''

    def __init__(self, filename):
        self.filename = filename
        self.atomtypes = []
        try:
            with open(self.filename, "r") as f:
                for line in f:
                    line = line.strip()
                    if line.startswith('#') or len(line) == 0:
                        continue
                    tok = line.split()
                    a_name = tok[0]
                    a_alpha = float(tok[4])
                    if next((x for x in self.atomtypes if x.name == a_name), None) is not None:
                        raise Exception(f'  error: atom type {a_name} is specified twice in {self.filename}')
                    self.atomtypes.append(AtomType(a_name, a_alpha))

        except IOError:
            print(f'  error: polarisation file {self.filename} not found')
            sys.exit(1)
        except Exception as e:
            print(e)
            sys.exit(1)

    def __str__(self):
        res = self.filename
        for at in self.atomtypes:
            res += '\n' + str(at)
        return res

class System(object):
    '''system consisting of fragments from frag.inp'''

    def __init__(self, filename):
        self.fragments = []
        self.fragmentpairs = []
        self.filename = filename

    def GetFragments(self, ff, p):
        '''read fragments from frag.inp'''

        try:
            with open(self.filename, 'r') as f:
                for line in f:
                    if line.startswith('#') or line.strip() == '':
                        continue
                    else:
                        tok = line.strip().split()
                        m_name = tok[0].lower()
                        if next((x for x in self.fragments if x.name == m_name), None) is not None:
                            raise Exception(f'  error: fragment {m_name} is specified twice in {self.filename}')
                        if len(tok)<2:
                            raise Exception(f'  error: no atomic names for {m_name} fragment in {self.filename}')
                        ind_range = tok[1:]
                        if next((x for x in ind_range if _Const.isfloat(x[0])), None) is not None:
                            raise Exception(f'  error: fragment {m_name} requires atomic labels not numbers in {self.filename}')
                    m = next((x for x in ff.monomers if x.name == m_name), None)
                    if m is None:
                        raise Exception(f'  error: fragment {m_name} not found in {ff.filename}')

                    self.fragments.append(Fragment(m,ind_range,Fragment.PolExclude(m,p)))

        except IOError:
            print(f'  error: fragment input file {self.filename} not found')
            sys.exit(1)
        except Exception as e:
            print(e)
            sys.exit(1)

    def GetFragmentPairs(self, ff):
        '''generate fragment pairs based on fragment list; fragment pairs are not formed for ++ and -- combinations'''

        try:
            i = j = 0
            for i in range(0,len(self.fragments)):
                for j in range(i,len(self.fragments)):
                    if ((self.fragments[i].q * self.fragments[j].q <= 0.0) or max(abs(self.fragments[i].q), abs(self.fragments[j].q)) == 0) and (self.fragments[i].scale_eps and self.fragments[j].scale_eps) and (not self.fragments[i].pol_model or not self.fragments[j].pol_model):
                        d = next((x for x in ff.dimers if (x.m1.name == self.fragments[i].name and x.m2.name == self.fragments[j].name) or (x.m1.name == self.fragments[j].name and x.m2.name == self.fragments[i].name)), None)

                        if d is not None:
                            self.fragmentpairs.append(FragmentPair(self.fragments[i], self.fragments[j], d.r, d.k_sapt))
                        else:
                            raise Exception(f'  error: dimer {self.fragments[i].name} {self.fragments[j].name} not found in {ff.filename}')
                    j+=1
                i+=1
        except Exception as e:
            print(e)
            sys.exit(1)

    def ParseScaleSigma(self, scsig):

        if len(scsig) > 0:
            if _Const.isfloat(scsig[0]):
                _Const.sigma_k = float(scsig[0])
                scsig.pop(0)

        if len(scsig) > 0:
            for i in scsig:
                fi = next((x for x in self.fragments if (x.name == i)), None)
                if (fi is None):
                    raise Exception(f'  error: fragment {i} specified with -s option not found in frag.inp ')
                elif fi.pol_model:
                    raise Exception(f'  error: fragment {i} specified with -s option is already polarisable: sigma should not be scaled')
                else:
                    fi.scale_sig = True
        else:
            for f in self.fragments:
                f.scale_sig = True


    def __str__(self):
        res = self.filename
        res += '\nFRAGMENTS'
        for fr in self.fragments:
            res += '\n' + str(fr)
        res+='\nFRAGMENT PAIRS'
        for frp in self.fragmentpairs:
            res += '\n' + str(frp)
        return res

    def GetFragAtoms(self, pol):
        '''get atom list for fragments based on existing atom type'''

        try:
            for fr in self.fragments:
                if os.path.exists('./' + fr.name + '.zmat'):
                    fr.GetAtomsFromZmat(pol)
                elif os.path.exists('./' + fr.name + '.xyz'):
                    fr.GetAtomsFromXYZ(pol)
                elif os.path.exists('./' + fr.name + '.mol'):
                    fr.GetAtomsFromMol(pol)
                elif os.path.exists('./' + fr.name + '.pdb'):
                    fr.GetAtomsFromPdb(pol)
                else:
                    raise Exception(f'  error: structure file (.zmat, .xyz, .pdb or .mol) for fragment {fr.name} not found')

        except Exception as e:
            print(e)
            sys.exit(1)

    def GetKPred(self):
        for frp in self.fragmentpairs:
            frp.PredictK()

def PrintReport(syst, sapt, scsig, polarisable):
    '''print report to screen'''

    report = 'Epsilon LJ coefficients were scaled by '
    if sapt:
        report += 'k_sapt'
    else:
        report += 'k_pred'
    report += ' parameter'
    if polarisable is not None:
        tmp = ', '.join(polarisable)
        report += f'. Fragments {tmp} were already polarisable'

    report += '\nSigma LJ coefficients '
    if scsig is None:
        report += 'were not scaled.'
    else:
        report += f'were scaled by {_Const.sigma_k:5.3f} value.'

        if all([x.scale_sig for x in syst.fragments]):
            report += ' All fragments were scaled.'
        else:
            sig_list = [f.name for f in syst.fragments if f.scale_sig]
            tmp = ', '.join(sig_list)
            report += f' Only {tmp} fragments were scaled.'

    report += '\n------------------------------------------\n'

    report += 'Fragment_i   Fragment_j   k_sapt'
    if not sapt:
        report += '    k_pred'

    for frp in syst.fragmentpairs:
        report += f'\n{frp.fr1.name:10s}   {frp.fr2.name:10s}'
        if sapt:
            report += f'  {frp.k_sapt:6.2f}'
        elif frp.k_sapt is None:
            report += f'        -    {frp.k_pred:6.2f}'
        else:
            report += f'   {frp.k_sapt:6.2f}    {frp.k_pred:6.2f}'

    report += '\n------------------------------------------'
    print(report)

def main():
    parser = argparse.ArgumentParser(description = usage, formatter_class = argparse.RawTextHelpFormatter)
    parser.add_argument('-f', '--ff_file', type=str, default = 'fragment.ff', help = 'fragment force field (default: fragment.ff)')
    parser.add_argument('-a', '--alpha_file', type=str, default = 'alpha.ff', help = 'polarisability values file (default: alpha.ff)')
    parser.add_argument('-i', '--id_file', type=str, default = 'frag.inp', help = 'identify atomic indices of fragments (default: frag.inp)')
    parser.add_argument('-ix', '--inxml', default = 'field-p.xml', help = 'input OpenMM xml file (default: field-p.xml)')
    parser.add_argument('-ox', '--outxml', default = 'field-p-sc.xml', help = 'output OpenMM xml file (default: field-p-sc.xml)')
    parser.add_argument('-q', '--sapt', action = 'store_true', help = 'use sapt calculated k values, default: use predicted k values')
    parser.add_argument('-s', '--scsig', nargs='*', type = str,  help = 'scale sigma if specified; default value: 0.985; \n\
    -s                       - scale all fragments\' sigma by 0.985 \n\
    -s value                 - scale all fragments\' sigma by user-defined value \n\
    -s name1 name2 ...       - scale the specified fragments\' sigma by 0.985 \n\
    -s value name1 name2 ... - scale the specified fragments\' sigma by user-defined value')
    parser.add_argument('-p', '--polarisable', nargs='+', type = str, help = 'already polarisable monomers')

    args = parser.parse_args()

    ff = FragmentFField(args.ff_file)
    syst = System(args.id_file)
    syst.GetFragments(ff,args.polarisable)

    if args.scsig is not None:
        syst.ParseScaleSigma(args.scsig)

    syst.GetFragmentPairs(ff)

    if (not args.sapt):
        pol = Polarisation(args.alpha_file)
        syst.GetFragAtoms(pol)
        syst.GetKPred()

    field = FieldXML(args.inxml)
    field.CheckAtoms(syst)
    field.Scale(ff,syst,args.sapt,args.scsig)
    field.Write(args.outxml)

    PrintReport(syst,args.sapt,args.scsig,args.polarisable)

if __name__ == '__main__':
    main()