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.) )