#!/usr/bin/env python from ezfio import ezfio import numpy as np import sys import time from datetime import datetime import matplotlib.pyplot as plt #import seaborn as sns from matplotlib.colors import LogNorm # ____________________________________________________________________________________________ # def read_AO_2eCoulomb(): lines = [] with open(ao_file) as f_in: for line in f_in: lines.append( float(line) ) J = np.zeros((ao_num2,ao_num2)) i_l = 0 for mu in range(ao_num): for nu in range(ao_num): ii = nu + ao_num * mu for sig in range(ao_num): for lam in range(ao_num): jj = lam + ao_num * sig J[ii,jj] = lines[i_l] i_l += 1 return(J) # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # def read_H_det(): H_det = np.zeros((n_det_alpha*n_det_beta,n_det_alpha*n_det_beta)) with open("fort.7000") as f_in: for line in f_in: ll = line.split() i = int( ll[0] ) - 1 j = int( ll[1] ) - 1 ii = j + n_det_beta * i k = int( ll[2] ) - 1 l = int( ll[3] ) - 1 jj = l + n_det_beta * k H_det[ii,jj] = float(ll[4]) return(H_det) # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # def read_T(): lines = [] with open("ij_T_kl.dat") as f_in: for line in f_in: lines.append( float(line) ) T = np.zeros( (n_det_alpha*n_det_beta,n_det_alpha*n_det_beta) ) i_l = 0 for i in range(n_det_alpha): for j in range(n_det_beta): ii = j + n_det_beta * i for k in range(n_det_alpha): for l in range(n_det_beta): jj = l + n_det_beta * k T[ii,jj] = lines[i_l] i_l += 1 return(T) # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # def check_symmetric(a, rtol=1e-05, atol=1e-08): return( np.allclose(a, a.T, rtol=rtol, atol=atol) ) # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # def plot_mat(M, xstep, ystep, namefile): Mp = np.abs(M) + 1e-15 ax = sns.heatmap(Mp, linewidths=0.0, rasterized=True, norm=LogNorm()) # x axis ax.xaxis.tick_top() xgrid = np.arange(1, M.shape[1]+1, xstep-1) plt.xticks( np.arange(0, M.shape[1], xstep-1) + .5, xgrid ) # y axis ygrid = np.arange(1, M.shape[0]+1, ystep-1) plt.yticks( np.arange(0, M.shape[0], ystep-1) + .5, ygrid ) plt.savefig(namefile) plt.clf() # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # def incomplete_Cholesky(M, ma, mb, theta_CD): D = np.zeros((ma,mb)) for lam in range(ma): for sig in range(mb): ii = sig + mb * lam D[sig,lam] = M[ii,ii] # find the maximum val_max = np.amax(D) ind_max = np.where(D == val_max) sig_max = ind_max[0][0] lam_max = ind_max[1][0] ii_max = sig_max + mb * lam_max #print( sig_max, lam_max, D[sig_max,lam_max], val_max ) # define the initial resudual mab = ma * mb R = np.zeros((mab)) R[:] = M[:,ii_max] m = 0 L = np.zeros((mab)) List_vec = [] while(val_max>theta_CD): if(m>0): # update the residual for ii in range(mab): sumii = 0 for i in range(m): sumii += List_vec[i][ii] * List_vec[i][ii_max] R[ii] = M[ii,ii_max] - sumii # compute Cholesky vector L = R / np.sqrt(val_max) List_vec.append(L) # update diagonals for lam in range(ma): for sig in range(mb): ii = sig + mb * lam D[sig,lam] = D[sig,lam] - L[ii]*L[ii] m += 1 # find the new maximum val_max = np.amax(D) ind_max = np.where(D == val_max) sig_max = ind_max[0][0] lam_max = ind_max[1][0] ii_max = sig_max + mb * lam_max return(m, List_vec) # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # def err_p_estim(M, List_vec): mM = M.shape[0] mL = len(List_vec) normM = np.linalg.norm(M) L_CD = np.zeros((mM,mL)) for i in range(mL): L_CD[:,i] = List_vec[i] print(" size of L_CD: {} MB\n".format( sys.getsizeof(L_CD) / 1024. ) ) err_p = 100. * np.linalg.norm( M - np.dot(L_CD, L_CD.T) ) / normM return( L_CD, err_p , normM) # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # def save_L_CD(List_vec, file_List): with open(file_List, 'w') as the_file: for i in range(List_vec.shape[0]): for j in range(List_vec.shape[1]): ff.write(" {:+e} \n".format(List_vec[i,j])) # ____________________________________________________________________________________________ # ____________________________________________________________________________________________ # ____________________________________________________________________________________________________________________________________ # ____________________________________________________________________________________________________________________________________ # if __name__=="__main__": t_beg = time.time() print("") print(" date and time = {} \n".format(datetime.now().strftime("%d/%m/%Y %H:%M:%S"))) EZFIO_name = "/home/aammar/qp2/src/svdwf/f2_work/f2_fci" ezfio.set_file(EZFIO_name) # ------------------------------------------------------------------------------------------------------ # 2-e integrals in AO basis # ------------------------------------------------------------------------------------------------------ ao_num = ezfio.get_ao_basis_ao_num() ao_num2 = ao_num * ao_num print(" nb of AO = {} \n".format(ao_num)) ao_file = "/home/aammar/qp2/src/svdwf/f2_work/ao_2eCoulombIntegrals.txt" J = read_AO_2eCoulomb() print(" size of J: {} MB ".format( sys.getsizeof(J) / 1024. ) ) print(" J is symmetric ? {} ".format(check_symmetric(J)) ) #plot_mat(J, ao_num, ao_num, 'J.pdf') # pivoted incomplete Cholesky decomposition theta_CD_J = 1e-02 m_J, List_vec_J = incomplete_Cholesky(J, ao_num, ao_num, theta_CD_J) print(" nb of L_CD vectors = {} for theta = {}".format(m_J,theta_CD_J)) # check error L_CD_J, err_p_J, normJ = err_p_estim(J, List_vec_J) print(" norm of J = {} \n".format(normJ)) print(" error estimation on J = {} % \n".format(err_p_J)) save_L_CD(L_CD_J, 'L_CD_J.txt') #plot_mat(L_CD_J, 2, ao_num, 'L_CD_J.pdf') # ------------------------------------------------------------------------------------------------------ # ------------------------------------------------------------------------------------------------------ print('\n \n \n') # ------------------------------------------------------------------------------------------------------ # < d_i d_j | T | d_k d_l > # ------------------------------------------------------------------------------------------------------ # n_det_alpha = ezfio.get_spindeterminants_n_det_alpha() # n_det_beta = ezfio.get_spindeterminants_n_det_beta () # n_ab = n_det_alpha * n_det_beta # print(" n_det_alpha = {} ".format(n_det_alpha)) # print(" n_det_beta = {} \n".format(n_det_beta )) # # T = read_T() # print(" T is symmetric ? {} \n".format(check_symmetric(T)) ) # # plot_mat(T, 50, 50, 'T.pdf') # # # # pivoted incomplete Cholesky decomposition # theta_CD_T = 1.5 # m_T, List_vec_T = incomplete_Cholesky(T, n_det_alpha, n_det_beta, theta_CD_T) # print(" nb of L_CD vectors = {} for theta = {}".format(m_T,theta_CD_T)) # # # check error # L_CD_T, err_p_T, normT = err_p_estim(T, List_vec_T) # print(" norm of T = {} ".format(normT)) # print(" error estimation on T = {} % \n".format(err_p_T)) # # plot_mat(L_CD_T, 2, n_det_beta, 'L_CD_T.pdf') # ------------------------------------------------------------------------------------------------------ # ------------------------------------------------------------------------------------------------------ print('\n \n \n') # ------------------------------------------------------------------------------------------------------ # < d_i d_j | H | d_k d_l > # # WARNING: this isn't semipositive definite # ------------------------------------------------------------------------------------------------------ # # # get H_det # n_det_alpha = ezfio.get_spindeterminants_n_det_alpha() # n_det_beta = ezfio.get_spindeterminants_n_det_beta () # n_ab = n_det_alpha * n_det_beta # print(" n_det_alpha = {} ".format(n_det_alpha)) # print(" n_det_beta = {} \n".format(n_det_beta )) # # H_det = read_H_det() # print(" size of H_det: {} MB ".format( sys.getsizeof(H_det) / 1024. ) ) # print(" H_det is symmetric ? {} \n".format(check_symmetric(H_det)) ) # # # plot H_det # H_det = H_det + 1e-15 # plot_mat(H_det, 50, 50, 'H_det.pdf') # # # # pivoted incomplete Cholesky decomposition # theta_CD_H_det = 1e-02 # m_H_det, List_vec_H_det = incomplete_Cholesky(H_det, n_det_alpha, n_det_beta, theta_CD_H_det) # print(" nb of L_CD_H_det vectors = {} for a precision of {}".format(m_H_det,theta_CD_H_det)) # #print(" size of L_CD_H_det: {} MB\n".format( sys.getsizeof(List_vec_H_det) / 1024. ) ) # # # check error # L_CD_H_det, err_p_H_det, normH_det = err_p_estim(H_det, List_vec_H_det) # print(" norm of H_det = {} \n".format(normH_det)) # print(" error estimation on H_det = {} % \n".format(err_p_H_det)) # # plot_mat(L_CD_H_det, 2, n_det_beta, 'L_CD_H_det.pdf') # ------------------------------------------------------------------------------------------------------ # ------------------------------------------------------------------------------------------------------ print('\n\n') t_end = time.time() print(' end after {} min'.format((t_end-t_beg)/60.)) # ____________________________________________________________________________________________________________________________________ # ____________________________________________________________________________________________________________________________________ # ____________________________________________________________________________________________________________________________________