QuantumPackage/scripts/qp_import_trexio.py

#!/usr/bin/env python3
"""
convert TREXIO file to EZFIO

Usage:
    qp_import_trexio [-o EZFIO_DIR] FILE

Options:
    -o --output=EZFIO_DIR    Produced directory
                             by default is FILE.ezfio

"""

import sys
import os
import numpy as np
from functools import reduce
from ezfio import ezfio
from docopt import docopt
import qp_bitmasks

try:
  import trexio
except ImportError:
    print("Error: trexio python module is not found. Try python3 -m pip install trexio")
    sys.exit(1)


try:
    QP_ROOT = os.environ["QP_ROOT"]
    QP_EZFIO = os.environ["QP_EZFIO"]
except KeyError:
    print("Error: QP_ROOT environment variable not found.")
    sys.exit(1)
else:
    sys.path = [QP_EZFIO + "/Python",
                QP_ROOT + "/install/resultsFile",
                QP_ROOT + "/install",
                QP_ROOT + "/scripts"] + sys.path

def uint64_to_int64(u):
    # Check if the most significant bit is set
    if u & (1 << 63):
        # Calculate the two's complement
        result = -int(np.bitwise_not(np.uint64(u))+1)
    else:
        # The number is already positive
        result = u
    return result

def generate_xyz(l):

    def create_z(x,y,z):
       return (x, y, l-(x+y))

    def create_y(accu,x,y,z):
       if y == 0:
          result = [create_z(x,y,z)] + accu
       else:
          result = create_y([create_z(x,y,z)] + accu , x, y-1, z)
       return result

    def create_x(accu,x,y,z):
       if x == 0:
          result = create_y([], x,y,z) + accu
       else:
          xnew = x-1
          ynew = l-xnew
          result = create_x(create_y([],x,y,z) + accu , xnew, ynew, z)
       return result

    result = create_x([], l, 0, 0)
    result.reverse()
    return result



def write_ezfio(trexio_filename, filename):

    try:
        trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_TEXT)
    except:
        trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_HDF5)

    ezfio.set_file(filename)
    ezfio.set_trexio_trexio_file(trexio_filename)

    print("Nuclei\t\t...\t", end=' ')

    charge = [0.]
    if trexio.has_nucleus(trexio_file):
        charge = trexio.read_nucleus_charge(trexio_file)
        ezfio.set_nuclei_nucl_num(len(charge))
        ezfio.set_nuclei_nucl_charge(charge)

        coord = trexio.read_nucleus_coord(trexio_file)
        coord = np.transpose(coord)
        ezfio.set_nuclei_nucl_coord(coord)

        label = trexio.read_nucleus_label(trexio_file)
        nucl_num = trexio.read_nucleus_num(trexio_file)

        # Transformt H1 into H
        import re
        p = re.compile(r'(\d*)$')
        label = [p.sub("", x).capitalize() for x in label]
        ezfio.set_nuclei_nucl_label(label)
        print("OK")

    else:
        ezfio.set_nuclei_nucl_num(1)
        ezfio.set_nuclei_nucl_charge([0.])
        ezfio.set_nuclei_nucl_coord([0.,0.,0.])
        ezfio.set_nuclei_nucl_label(["X"])
        print("None")



    print("Electrons\t...\t", end=' ')

    try:
        num_beta = trexio.read_electron_dn_num(trexio_file)
    except:
        num_beta = int(sum(charge))//2

    try:
        num_alpha = trexio.read_electron_up_num(trexio_file)
    except:
        num_alpha = int(sum(charge)) - num_beta

    if num_alpha == 0:
        print("\n\nError: There are zero electrons in the TREXIO file.\n\n")
        sys.exit(1)
    ezfio.set_electrons_elec_alpha_num(num_alpha)
    ezfio.set_electrons_elec_beta_num(num_beta)

    print(f"{num_alpha} {num_beta}")

    print("Basis\t\t...\t", end=' ')

    shell_num = 0
    try:
        basis_type = trexio.read_basis_type(trexio_file)

        print ("BASIS TYPE: ", basis_type.lower())
        if basis_type.lower() in ["gaussian", "slater"]:
            shell_num   = trexio.read_basis_shell_num(trexio_file)
            prim_num    = trexio.read_basis_prim_num(trexio_file)
            ang_mom     = trexio.read_basis_shell_ang_mom(trexio_file)
            nucl_index  = trexio.read_basis_nucleus_index(trexio_file)
            exponent    = trexio.read_basis_exponent(trexio_file)
            coefficient = trexio.read_basis_coefficient(trexio_file)
            shell_index = trexio.read_basis_shell_index(trexio_file)
            ao_shell    = trexio.read_ao_shell(trexio_file)

            ezfio.set_basis_basis("Read from TREXIO")
            ezfio.set_ao_basis_ao_basis("Read from TREXIO")
            ezfio.set_basis_shell_num(shell_num)
            ezfio.set_basis_prim_num(prim_num)
            ezfio.set_basis_shell_ang_mom(ang_mom)
            ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
            ezfio.set_basis_prim_expo(exponent)
            ezfio.set_basis_prim_coef(coefficient)

            nucl_shell_num = []
            prev = None
            m = 0
            for i in ao_shell:
                if i != prev:
                   m += 1
                   if prev is None or nucl_index[i] != nucl_index[prev]:
                        nucl_shell_num.append(m)
                        m = 0
                prev = i
            assert (len(nucl_shell_num) == nucl_num)

            shell_prim_num = []
            prev = shell_index[0]
            count = 0
            for i in shell_index:
                if i != prev:
                   shell_prim_num.append(count)
                   count = 0
                count += 1
                prev = i
            shell_prim_num.append(count)

            assert (len(shell_prim_num) == shell_num)

            ezfio.set_basis_shell_prim_num(shell_prim_num)
            ezfio.set_basis_shell_index([x+1 for x in shell_index])
            ezfio.set_basis_nucleus_shell_num(nucl_shell_num)


            shell_factor = trexio.read_basis_shell_factor(trexio_file)
            prim_factor  = trexio.read_basis_prim_factor(trexio_file)

        elif basis_type.lower() == "numerical":

            shell_num   = trexio.read_basis_shell_num(trexio_file)
            prim_num    = shell_num
            ang_mom     = trexio.read_basis_shell_ang_mom(trexio_file)
            nucl_index  = trexio.read_basis_nucleus_index(trexio_file)
            exponent    = [1.]*prim_num
            coefficient = [1.]*prim_num
            shell_index = [i for i in range(shell_num)]
            ao_shell    = trexio.read_ao_shell(trexio_file)

            ezfio.set_basis_basis("None")
            ezfio.set_ao_basis_ao_basis("None")
            ezfio.set_basis_shell_num(shell_num)
            ezfio.set_basis_prim_num(prim_num)
            ezfio.set_basis_shell_ang_mom(ang_mom)
            ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
            ezfio.set_basis_prim_expo(exponent)
            ezfio.set_basis_prim_coef(coefficient)

            nucl_shell_num = []
            prev = None
            m = 0
            for i in ao_shell:
                if i != prev:
                   m += 1
                   if prev is None or nucl_index[i] != nucl_index[prev]:
                        nucl_shell_num.append(m)
                        m = 0
                prev = i
            assert (len(nucl_shell_num) == nucl_num)

            shell_prim_num = []
            prev = shell_index[0]
            count = 0
            for i in shell_index:
                if i != prev:
                   shell_prim_num.append(count)
                   count = 0
                count += 1
                prev = i
            shell_prim_num.append(count)

            assert (len(shell_prim_num) == shell_num)

            ezfio.set_basis_shell_prim_num(shell_prim_num)
            ezfio.set_basis_shell_index([x+1 for x in shell_index])
            ezfio.set_basis_nucleus_shell_num(nucl_shell_num)

            shell_factor = trexio.read_basis_shell_factor(trexio_file)
            prim_factor  = [1.]*prim_num
        else:
           raise TypeError

        print(basis_type)
    except:
        print("None")
        ezfio.set_ao_basis_ao_cartesian(True)

    print("AOS\t\t...\t", end=' ')

    try:
        cartesian = trexio.read_ao_cartesian(trexio_file)
    except:
        cartesian = True

    if not cartesian:
        raise TypeError('Only cartesian TREXIO files can be converted')

    ao_num = trexio.read_ao_num(trexio_file)
    ezfio.set_ao_basis_ao_num(ao_num)

    if shell_num > 0:
        ao_shell    = trexio.read_ao_shell(trexio_file)
        at = [ nucl_index[i]+1 for i in ao_shell ]
        ezfio.set_ao_basis_ao_nucl(at)

        num_prim0 = [ 0 for i in range(shell_num) ]
        for i in shell_index:
            num_prim0[i] += 1

        coef = {}
        expo = {}
        for i,c in enumerate(coefficient):
            idx = shell_index[i]
            if idx in coef:
              coef[idx].append(c)
              expo[idx].append(exponent[i])
            else:
              coef[idx] = [c]
              expo[idx] = [exponent[i]]

        coefficient = []
        exponent    = []
        power_x     = []
        power_y     = []
        power_z     = []
        num_prim    = []

        for i in range(shell_num):
            for x,y,z in generate_xyz(ang_mom[i]):
                power_x.append(x)
                power_y.append(y)
                power_z.append(z)
                coefficient.append(coef[i])
                exponent.append(expo[i])
                num_prim.append(num_prim0[i])

        assert (len(coefficient) == ao_num)
        ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
        ezfio.set_ao_basis_ao_prim_num(num_prim)

        prim_num_max = max( [ len(x) for x in coefficient ] )

        for i in range(ao_num):
            coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)]
            exponent   [i] += [0. for j in range(len(exponent[i]), prim_num_max)]

        coefficient = reduce(lambda x, y: x + y, coefficient, [])
        exponent    = reduce(lambda x, y: x + y, exponent   , [])

        coef = []
        expo = []
        for i in range(prim_num_max):
            for j in range(i, len(coefficient), prim_num_max):
                coef.append(coefficient[j])
                expo.append(exponent[j])

#        ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
        ezfio.set_ao_basis_ao_coef(coef)
        ezfio.set_ao_basis_ao_expo(expo)

        print("OK")

    else:
        print("None")


    #                _
    # |\/|  _   _   |_)  _.  _ o  _
    # |  | (_) _>   |_) (_| _> | _>
    #

    print("MOS\t\t...\t", end=' ')

    labels = { "Canonical" : "Canonical",
               "RHF" : "Canonical",
               "BOYS" : "Localized",
               "ROHF" : "Canonical",
               "UHF" : "Canonical",
               "Natural": "Natural" }
    try:
      label = labels[trexio.read_mo_type(trexio_file)]
    except:
      label = "None"
    ezfio.set_mo_basis_mo_label(label)
    ezfio.set_determinants_mo_label(label)

    try:
      clss = trexio.read_mo_class(trexio_file)
      core     = [ i for i in clss if i.lower() == "core" ]
      inactive = [ i for i in clss if i.lower() == "inactive" ]
      active   = [ i for i in clss if i.lower() == "active" ]
      virtual  = [ i for i in clss if i.lower() == "virtual" ]
      deleted  = [ i for i in clss if i.lower() == "deleted" ]
    except trexio.Error:
      pass

    try:
      mo_num = trexio.read_mo_num(trexio_file)
      ezfio.set_mo_basis_mo_num(mo_num)

      MoMatrix = trexio.read_mo_coefficient(trexio_file)
      ezfio.set_mo_basis_mo_coef(MoMatrix)

      mo_occ = [ 0. for i in range(mo_num) ]
      for i in range(num_alpha):
         mo_occ[i] += 1.
      for i in range(num_beta):
         mo_occ[i] += 1.
      ezfio.set_mo_basis_mo_occ(mo_occ)
      print("OK")
    except:
      print("None")



    print("Pseudos\t\t...\t", end=' ')

    ezfio.set_pseudo_do_pseudo(False)

    if trexio.has_ecp_ang_mom(trexio_file):
        ezfio.set_pseudo_do_pseudo(True)
        max_ang_mom_plus_1 = trexio.read_ecp_max_ang_mom_plus_1(trexio_file)
        z_core = trexio.read_ecp_z_core(trexio_file)
        ang_mom = trexio.read_ecp_ang_mom(trexio_file)
        nucleus_index = trexio.read_ecp_nucleus_index(trexio_file)
        exponent = trexio.read_ecp_exponent(trexio_file)
        coefficient = trexio.read_ecp_coefficient(trexio_file)
        power = trexio.read_ecp_power(trexio_file)

        lmax = max( max_ang_mom_plus_1 ) - 1
        ezfio.set_pseudo_pseudo_lmax(lmax)
        ezfio.set_pseudo_nucl_charge_remove(z_core)

        prev_center = None
        ecp = {}
        for i in range(len(ang_mom)):
            center = nucleus_index[i]
            if center != prev_center:
               ecp[center] = { "lmax": max_ang_mom_plus_1[center],
                               "zcore": z_core[center],
                               "contr": {} }
               for j in range(max_ang_mom_plus_1[center]+1):
                    ecp[center]["contr"][j] = []

            ecp[center]["contr"][ang_mom[i]].append( (coefficient[i], power[i], exponent[i]) )
            prev_center = center

        ecp_loc = {}
        ecp_nl  = {}
        kmax    = 0
        klocmax    = 0
        for center in ecp:
            ecp_nl [center] = {}
            for k in ecp[center]["contr"]:
                if k == ecp[center]["lmax"]:
                    ecp_loc[center] = ecp[center]["contr"][k]
                    klocmax = max(len(ecp_loc[center]), klocmax)
                else:
                    ecp_nl [center][k] = ecp[center]["contr"][k]
                    kmax = max(len(ecp_nl [center][k]), kmax)

        ezfio.set_pseudo_pseudo_klocmax(klocmax)
        ezfio.set_pseudo_pseudo_kmax(kmax)

        pseudo_n_k = [[0  for _ in range(nucl_num)] for _ in range(klocmax)]
        pseudo_v_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)]
        pseudo_dz_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)]
        pseudo_n_kl = [[[0  for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
        pseudo_v_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
        pseudo_dz_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
        for center in ecp_loc:
            for k in range( len(ecp_loc[center]) ):
                v, n, dz = ecp_loc[center][k]
                pseudo_n_k[k][center] = n
                pseudo_v_k[k][center] = v
                pseudo_dz_k[k][center] = dz

        ezfio.set_pseudo_pseudo_n_k(pseudo_n_k)
        ezfio.set_pseudo_pseudo_v_k(pseudo_v_k)
        ezfio.set_pseudo_pseudo_dz_k(pseudo_dz_k)

        for center in ecp_nl:
            for l in range( len(ecp_nl[center]) ):
                for k in range( len(ecp_nl[center][l]) ):
                    v, n, dz = ecp_nl[center][l][k]
                    pseudo_n_kl[l][k][center] = n
                    pseudo_v_kl[l][k][center] = v
                    pseudo_dz_kl[l][k][center] = dz

        ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl)
        ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl)
        ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl)
        print("OK")

    else:
        print("None")

    print("Determinant\t...\t", end=' ')
    alpha = [ i for i in range(num_alpha) ]
    beta  = [ i for i in range(num_beta) ]
    if trexio.has_mo_spin(trexio_file):
       spin = trexio.read_mo_spin(trexio_file)
       if max(spin) == 1:
         alpha = [ i for i in range(len(spin)) if spin[i] == 0 ]
         alpha = [ alpha[i] for i in range(num_alpha) ]
         beta  = [ i for i in range(len(spin)) if spin[i] == 1 ]
         beta  = [ beta[i] for i in range(num_beta) ]
         print("Warning -- UHF orbitals --", end=' ')
    alpha_s = ['0']*mo_num
    beta_s  = ['0']*mo_num
    for i in alpha:
      alpha_s[i] = '1'
    for i in beta:
      beta_s[i] = '1'
    alpha_s = ''.join(alpha_s)[::-1]
    beta_s = ''.join(beta_s)[::-1]
    def conv(i):
      try:
        result = np.int64(i)
      except:
        result = np.int64(i-2**63-1)
      return result

    alpha = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(alpha_s) ][::-1]
    beta  = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(beta_s ) ][::-1]
    ezfio.set_determinants_bit_kind(8)
    ezfio.set_determinants_n_int(1+mo_num//64)
    ezfio.set_determinants_n_det(1)
    ezfio.set_determinants_n_states(1)
    ezfio.set_determinants_psi_det(alpha+beta)
    ezfio.set_determinants_psi_coef([[1.0]])
    print("OK")




def get_full_path(file_path):
    file_path = os.path.expanduser(file_path)
    file_path = os.path.expandvars(file_path)
    return file_path


if __name__ == '__main__':
    ARGUMENTS = docopt(__doc__)

    FILE = get_full_path(ARGUMENTS['FILE'])
    trexio_filename = FILE

    if ARGUMENTS["--output"]:
        EZFIO_FILE = get_full_path(ARGUMENTS["--output"])
    else:
        EZFIO_FILE = "{0}.ezfio".format(FILE)

    write_ezfio(trexio_filename, EZFIO_FILE)
    sys.stdout.flush()