qp_plugins_scemama/devel/svdwf/buildpsi_diagSVDit.py

import sys, os
QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")

from ezfio import ezfio
from datetime import datetime
import time
import numpy as np
import subprocess
from scipy.linalg import eig, eigh


# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
def get_energy():
    buffer = subprocess.check_output( ['qmcchem', 'result', '-e', 'e_loc', filename]
                                      , encoding='UTF-8' )
    if buffer.strip() != "":
        buffer = buffer.splitlines()[-1]
        _, energy, error = [float(x) for x in buffer.split()]
        return energy, error
    else:
      return None, None
# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def run_qmc():
    return subprocess.check_output(['qmcchem', 'run', filename])

def stop_qmc():
    subprocess.check_output(['qmcchem', 'stop', filename])

def set_vmc_params():
    #subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
    #                                 '-m', 'VMC',
    #                                 '-l', str(20),
    #                                 '--time-step=0.3',
    #                                 '--stop-time=36000',
    #                                 '--norm=1.e-5',
    #                                 '-w', '10',
    #                                 filename])
    subprocess.check_output(['qmcchem', 'edit', '-c'
                                      , '-j', 'None'
                                      , '-l', str(block_time)
                                      , '-t', str(total_time)
                                      , filename])

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !





# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def get_Aref():

    Aref = np.zeros( (n_alpha, n_beta) )
    for k in range(n_det):
        i         = A_rows[k] - 1
        j         = A_cols[k] - 1
        Aref[i,j] = A_vals[0][k]
    return( Aref )

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def check_svd_error( A, U, S, V ):

    # || A - U x S x transpose(V) ||

    # to normalize
    norm_A = np.linalg.norm(A, ord='fro')
    
    # vector S ==> matrix S
    _, na  = U.shape
    _, nb  = V.shape
    S_mat  = np.zeros( (na,nb) )
    for i in range( min(na,nb)  ):
        S_mat[i,i] = S[i]

    #A_SVD   = np.linalg.multi_dot([ U, np.diag(S), Vt ])
    A_SVD   = np.linalg.multi_dot([ U, S_mat, V.T ])
    err_SVD = 100. * np.linalg.norm( A - A_SVD, ord="fro") / norm_A

    print(' error between A_SVD and Aref = {} %\n'.format(err_SVD) )

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !




# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def FSVD_save_EZFIO():

    U_toEZFIO = np.zeros( ( 1, U_FSVD.shape[1], U_FSVD.shape[0] ) )
    V_toEZFIO = np.zeros( ( 1, V_FSVD.shape[1], V_FSVD.shape[0] ) )
    U_toEZFIO[0,:,:] = U_FSVD.T
    V_toEZFIO[0,:,:] = V_FSVD.T

    ezfio.set_spindeterminants_psi_svd_alpha_unique( U_toEZFIO )
    ezfio.set_spindeterminants_psi_svd_beta_unique ( V_toEZFIO )
    ezfio.set_spindeterminants_psi_svd_coefs_unique( S_FSVD )

    print(' Full SVD unique vectors & coeff are saved to EZFIO ')

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def numerote_selected_vectors():
    
    numalpha_selected = []
    numbeta_selected  = []

    for i in range(n_TSVD): 
        for j in range(n_TSVD):
            numalpha_selected.append(j+1)
            numbeta_selected.append(i+1)

    if( (len(numalpha_selected)!=n_selected) or (len(numbeta_selected)!=n_selected) ) :
        print(' error in numerating selectod vectors')
        print(' {} != {} != {}'.format(n_selected,len(numalpha_selected),len(numbeta_selected)) )
    else:
        ezfio.set_spindeterminants_psi_svd_alpha_numselected(numalpha_selected)
        ezfio.set_spindeterminants_psi_svd_beta_numselected (numbeta_selected)
        print(' selected vectors are numeroted in EZFIO')

    return( numalpha_selected,numbeta_selected  )

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def numerote_toselect_vectors( choix  ):
    
    numalpha_toselect = []
    numbeta_toselect  = []

    if( choix == 3 ):
        # nondiagonal blocs
        for i in range(n_TSVD): 
            for j in range(n_TSVD,n_beta):
                numalpha_toselect.append(i+1)
                numbeta_toselect.append(j+1)
        for i in range(n_TSVD,n_alpha):
            for j in range(n_TSVD): 
                numalpha_toselect.append(i+1)
                numbeta_toselect.append(j+1)
        # diagonal bloc
        for i in range(n_TSVD,n_alpha):
            for j in range(n_TSVD,n_beta):
                numalpha_toselect.append(i+1)
                numbeta_toselect.append(j+1)
    elif( choix == 2 ):
        # nondiagonal blocs
        for i in range(n_TSVD): 
            for j in range(n_TSVD,n_beta):
                numalpha_toselect.append(i+1)
                numbeta_toselect.append(j+1)
        for i in range(n_TSVD,n_alpha):
            for j in range(n_TSVD): 
                numalpha_toselect.append(i+1)
                numbeta_toselect.append(j+1)
    else:
        print(' choix = 2 ou 3' )
        exit()

    if( (len(numalpha_toselect)!=n_toselect) or (len(numbeta_toselect)!=n_toselect) ) :
        print(' error in numerating vectors to select')
        print(' {} != {} != {}'.format(n_toselect,len(numalpha_toselect),len(numbeta_toselect)) )
    else:
        ezfio.set_spindeterminants_psi_svd_alpha_numtoselect(numalpha_toselect)
        ezfio.set_spindeterminants_psi_svd_beta_numtoselect (numbeta_toselect)
        print(' vectors to select are numeroted in EZFIO')

    return( numalpha_toselect,numbeta_toselect  )

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def normalize_S_TSVD():

    global S_TSVD

    norm_S_TSVD = np.linalg.norm(S_TSVD, ord='fro')
    S_TSVD = S_TSVD / norm_S_TSVD

    return()

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def get_h_selected_matrix():

    h_selected_matrix = np.zeros( (n_selected,n_selected) )
    h_selected_stater = np.zeros( (n_selected,n_selected) )

    beg_h_selected_matrix = results.find('h_selected_matrix : [ ') + len('h_selected_matrix : [ ')
    end_h_selected_matrix = len(results)
    h_selected_matrix_buf = results[beg_h_selected_matrix:end_h_selected_matrix]
    h_selected_matrix_buf = h_selected_matrix_buf.split( '\n' )

    for i in range(1,n_selected+1):
        ii0 = (i-1) * n_selected
        for j in range(1,n_selected+1):
            iline = ii0 + j

            line = h_selected_matrix_buf[iline].split()
            indc = int  ( line[0] )

            if( indc != iline ):
                print('Error in get_h_selected_matrix')
                exit()
            else:
                h_selected_matrix[i-1][j-1] = float( line[2] )
                h_selected_stater[i-1][j-1] = float( line[4] )

    return(h_selected_matrix,h_selected_stater)

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def get_o_selected_matrix():

    o_selected_matrix = np.zeros( (n_selected,n_selected) )
    o_selected_stater = np.zeros( (n_selected,n_selected) )

    beg_o_selected_matrix = results.find('overlop_selected_matrix : [ ') + len('overlop_selected_matrix : [ ')
    end_o_selected_matrix = len(results)
    o_selected_matrix_buf = results[beg_o_selected_matrix:end_o_selected_matrix]
    o_selected_matrix_buf = o_selected_matrix_buf.split( '\n' )

    for i in range(1,n_selected+1):
        ii0 = (i-1) * n_selected
        for j in range(1,n_selected+1):
            iline = ii0 + j

            line = o_selected_matrix_buf[iline].split()
            indc = int  ( line[0] )

            if( indc != iline ):
                print('Error in get_o_selected_matrix')
                exit()
            else:
                o_selected_matrix[i-1][j-1] = float( line[2] )
                o_selected_stater[i-1][j-1] = float( line[4] )

    return(o_selected_matrix,o_selected_stater)

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def get_Epostsvd():

    # symmetrise and diagonalise
    aa = h_selected_matrix
    aa = 0.5*( aa + aa.T )
    bb = o_selected_matrix
    eigvals_postsvd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True,
            check_finite=True, homogeneous_eigvals=False)

    d_postsvd = np.diagflat(S_TSVD)
    d_postsvd = d_postsvd.reshape( (1,n_selected*n_selected) )
    recouvre_postsvd = np.abs(d_postsvd @ vr)
    ind_gspostsvd = np.argmax(recouvre_postsvd)

    E_postsvd = eigvals_postsvd[ind_gspostsvd]

    return( E_postsvd,  vr[:,ind_gspostsvd] )

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def Hqmc_svd_diag():

    # read CI_SVD matrices
    Ci_h_matrix_svd       = get_Ci_h_matrix_svd()
    Ci_overlap_matrix_svd = get_Ci_overlap_matrix_svd()

    # symmetrise
    aa = Ci_h_matrix_svd
    aa = 0.5*( aa + aa.T )

    bb = Ci_overlap_matrix_svd

    # diagonalise
    eigvals_svd, vr = eig( aa, bb, left=False, right=True
                         , overwrite_a=True, overwrite_b=True
                         , check_finite=True, homogeneous_eigvals=False)

    recouvre_svd = np.abs( np.dot( S_FSVD, vr) )
    ind_gssvd    = np.argmax(recouvre_svd)
    E_svd        = eigvals_svd[ind_gssvd] + nuc_energy

    return( E_svd, vr[:,ind_gssvd] )

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def SVD_postsvd():

    # reshape sigma_postsvd
    sigma_postsvd_mat = np.zeros( (n_selected,n_selected) )
    for i in range(n_TSVD):
        ii = i*n_TSVD
        for j in range(n_TSVD):
            jj = i*n_TSVD + j
            sigma_postsvd_mat[i][j] = sigma_postsvd[jj]

    # construct the new matrix Y & perform a new SVD
    Y = np.dot( U_TSVD , np.dot(sigma_postsvd_mat , Vt_TSVD) )
    U, S, Vt = np.linalg.svd(Y, full_matrices=True)

    return(U, S, Vt)

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def get_hij_fm():

    hij_fm        = np.zeros( (n_toselect) )
    hij_fm_stater = np.zeros( (n_toselect) )

    beg_hij = results.find('hij_fm : [ ') + len('hij_fm : [ ')
    end_hij = len(results)
    hij_buf = results[beg_hij:end_hij]
    hij_buf = hij_buf.split( '\n' )

    for iline in range(1,n_toselect+1):
        line = hij_buf[iline].split()
        indc = int( line[0] )
        if( indc != iline ):
            print('Error in get_hij_fm')
            exit()
        else:
            hij_fm       [iline-1] = float( line[2] )
            hij_fm_stater[iline-1] = float( line[4] )

    return(hij_fm, hij_fm_stater)

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def get_hij_sm():

    hij_sm        = np.zeros( (n_toselect) )
    hij_sm_stater = np.zeros( (n_toselect) )

    beg_hij = results.find('hij_sm : [ ') + len('hij_sm : [ ')
    end_hij = len(results)
    hij_buf = results[beg_hij:end_hij]
    hij_buf = hij_buf.split( '\n' )

    for iline in range(1,n_toselect+1):
        line = hij_buf[iline].split()
        indc = int( line[0] )
        if( indc != iline ):
            print('Error in get_hij_sm')
            exit()
        else:
            hij_sm       [iline-1] = float( line[2] )
            hij_sm_stater[iline-1] = float( line[4] )

    return(hij_sm, hij_sm_stater)

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !



# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !

def get_xij_diag():

    xij_diag        = np.zeros( (n_toselect) )
    xij_diag_stater = np.zeros( (n_toselect) )

    beg_xij = results.find('xij_diag : [ ') + len('xij_diag : [ ')
    end_xij = len(results)
    xij_buf = results[beg_xij:end_xij]
    xij_buf = xij_buf.split( '\n' )

    for iline in range(1,n_toselect+1):
        line = xij_buf[iline].split()
        indc = int( line[0] )
        if( indc != iline ):
            print('Error in get_xij_diag')
            exit()
        else:
            xij_diag       [iline-1] = float( line[2] )
            xij_diag_stater[iline-1] = float( line[4] )

    return(xij_diag, xij_diag_stater)

# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !





if __name__ == '__main__':

    t0 = time.time()
    print("Today's date:", datetime.now() )

    # EZFIO file
    filename = "/home/aammar/qp2/src/svdwf/2h2_work/2h2_cisd"
    ezfio.set_file(filename)
    print("filename = {}".format(filename))

    # parameters
    energ_nuc = 1.711353545183182	# for 2h2 

    # get spindeterminant data from EZFIO
    n_alpha = ezfio.get_spindeterminants_n_det_alpha()
    n_beta  = ezfio.get_spindeterminants_n_det_beta()
    A_rows  = np.array(ezfio.get_spindeterminants_psi_coef_matrix_rows()   )
    A_cols  = np.array(ezfio.get_spindeterminants_psi_coef_matrix_columns())
    A_vals  = np.array(ezfio.get_spindeterminants_psi_coef_matrix_values() )
    n_det   = A_rows.shape[0]

    print('~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~')
    print(' matrix: {} x {}'.format(n_alpha,n_beta))
    print(' n_det = {}'.format(n_det))
    print('~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~')

    # construct intial dense matrix
    Aref = get_Aref()

    # perform initial Full SVD
    print('')
    print(' ----- Performing Full SVD ----- ')
    U_FSVD, S_FSVD, Vt_FSVD  = np.linalg.svd(Aref, full_matrices=True)
    V_FSVD = Vt_FSVD.T

    # check Full SVD error
    check_svd_error( Aref, U_FSVD, S_FSVD, V_FSVD )

    # save SVD vectors & coefficients in EZFIO 
    ezfio.set_spindeterminants_n_svd_coefs_unique(min(n_alpha,n_beta))
    FSVD_save_EZFIO()

    # truncated SVD
    n_TSVD = 15
    ezfio.set_spindeterminants_n_svd_coefs(n_TSVD) 
    U_TSVD = U_FSVD[:,:n_TSVD]
    V_TSVD = V_FSVD[:,:n_TSVD]
    S_TSVD = S_FSVD[:n_TSVD]

    # check truncataed SVD error
    check_svd_error( Aref, U_TSVD, S_TSVD, V_TSVD )

    # numerote selected vectors & save to EZFIO
    n_selected = n_TSVD * n_TSVD
    print(' n_selected = {}'.format(n_selected))
    ezfio.set_spindeterminants_n_svd_selected(n_selected)
    numalpha_selected,numbeta_selected = numerote_selected_vectors()

    # numerote vectors to select & save to EZFIO
    n_toselect = n_alpha * n_beta - n_selected 
    print(' n_toselect = {}'.format(n_toselect))
    ezfio.set_spindeterminants_n_svd_toselect(n_toselect)
    numalpha_toselect,numbeta_toselect = numerote_toselect_vectors(choix=3)


    #_________________________________________________________________________________________
    #
    #                          loop over SVD iterations
    #_________________________________________________________________________________________
    
    it_svd     = 0
    it_svd_max = 1

    while( it_svd < it_svd_max ):

        it_svd = it_svd + 1

        # normalize sigular values in the truncated space
        #normalize_S_TSVD()

        # run QMC to get H_postsvd
        block_time = 300	# in sec
        total_time = 1800	# in sec
        set_vmc_params()
        ezfio.set_properties_hij_fm(False)
        ezfio.set_properties_hij_sm(False)
        ezfio.set_properties_xij_diag(False)
        ezfio.set_properties_h_selected_matrix(True)
        ezfio.set_properties_overlop_selected_matrix(True)
        run_qmc()

        # read QMC=CHEM results
        t_read = time.time()
        print(' getting QMCCHEM results from {}'.format(EZFIO_file) )
        results =  subprocess.check_output(['qmcchem', 'result', EZFIO_file, '>> results.dat'],  encoding='UTF-8')
        print(' getting results after {} minutes \n'.format( (time.time()-t_read)/60. ))

        # < E_loc >
        Eloc, ErrEloc = get_energy()
        print(' Eloc = {} +/- {}'.format(Eloc, ErrEloc))

        # get H and overlop from QMC=CHEM
        h_selected_matrix, h_selected_stater = get_h_selected_matrix()
        o_selected_matrix, o_selected_stater = get_o_selected_matrix()

        # ground state from H S = E S
        E_postsvd, sigma_postsvd = get_Epostsvd()
        print(' post svd energy = {}'.format(E_postsvd+energ_nuc) )
        
        # perform new SVD: U x sigma_postsvd x Vt --> U' x S' x Vt'
        U_FSVD, S_FSVD, Vt_FSVD = SVD_postsvd()
        V_FSVD = Vt_FSVD.T

        # save in EZFIO
        FSVD_save_EZFIO()

        # run QMC to get hij, e and xij_diag from QMC=CHEM
        block_time = 300	# in sec
        total_time = 1800	# in sec
        set_vmc_params()
        ezfio.set_properties_h_selected_matrix(False)
        ezfio.set_properties_overlop_selected_matrix(False)
        ezfio.set_properties_hij_fm(True)
        ezfio.set_properties_hij_sm(True)
        ezfio.set_properties_xij_diag(True)
        run_qmc()

        # read QMC=CHEM results
        t_read = time.time()
        print(' getting QMCCHEM results from {}'.format(EZFIO_file) )
        results =  subprocess.check_output(['qmcchem', 'result', EZFIO_file, '>> results.dat'],  encoding='UTF-8')
        print(' getting results after {} minutes \n'.format( (time.time()-t_read)/60. ))

        # < E_loc >
        Eloc, ErrEloc = get_energy()
        print(' Eloc = {} +/- {}'.format(Eloc, ErrEloc))

        # hij  =  < psi_svd J | H | J l l' > / < psi_svd J | psi_svd J > 
        #      =  < H (J l l')/(psi_svd J) >    ( first  method: fm )
        #      =  < E_loc (l l') / psi_svd >    ( second method: sm )
        hij_fm, hif_fm_stater = get_hij_fm()
        hij_sm, hif_sm_stater = get_hij_sm()

        # get xij_diag
        xij_diag, xij_diag_stater = get_xij_diag()
        
        # first method
        dij_fm        = np.zeros( (n_toselect) )
        dij_fm_stater = np.zeros( (n_toselect) )
        for i in range(n_toselect):
            dij_fm[i] = hij_fm[i] / ( Eloc - xij_diag[i] ) 
            # statistic error 
            a = ( hij_fm_stater[i]*hij_fm_stater[i] ) / ( hij_fm[i]*hij_fm[i] )
            b = ( ErrEloc*ErrEloc + xij_diag_stater[i]*xij_diag_stater[i] ) / ( (Eloc-xij_diag[i])*(Eloc-xij_diag[i]) ) 
            dij_fm_stater[i] = abs(dij_fm[i]) * np.sqrt( a + b )

        cc = np.concatenate((dij_fm,dij_fm_stater),axis=1)
        np.savetxt('dij_fm.txt',cc)

        # second method
        dij_sm        = np.zeros( (n_toselect) )
        dij_sm_stater = np.zeros( (n_toselect) )
        for i in range(n_toselect):
            dij_sm[i] = hij_sm[i] / ( Eloc - xij_diag[i] ) 
            # statistic error 
            a = ( hij_fm_stater[i]*hij_fm_stater[i] ) / ( hij_sm[i]*hij_sm[i] )
            b = ( ErrEloc*ErrEloc + xij_diag_stater[i]*xij_diag_stater[i] ) / ( (Eloc-xij_diag[i])*(Eloc-xij_diag[i]) ) 
            dij_sm_stater[i] = abs(dij_sm[i]) * np.sqrt( a + b )

        cc = np.concatenate((dij_sm,dij_sm_stater),axis=1)
        np.savetxt('dij_sm.txt',cc)

        # TODO
        # choose fm or sm & perform a new SVD 


    #_________________________________________________________________________________________

    print("end after {:.3f} minutes".format((time.time()-t0)/60.) )