From f1b2c902e1b187d9a21dd4a59173d97c3e7ef641 Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Fri, 15 Jan 2016 16:59:08 +0100 Subject: [PATCH 01/29] Beter anouar script --- scripts/qp_convert_qmcpack_from_ezfio.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/qp_convert_qmcpack_from_ezfio.py b/scripts/qp_convert_qmcpack_from_ezfio.py index 7443be68..9d194314 100755 --- a/scripts/qp_convert_qmcpack_from_ezfio.py +++ b/scripts/qp_convert_qmcpack_from_ezfio.py @@ -48,7 +48,7 @@ print "Atomic coord in Bohr" for i,t in enumerate(zip(l_label,l_charge,l_coord_str)): try : - l = (t[0],t[1]+zcore[i],t[1]) + l = (t[0],t[1]+zcore[i],t[2]) except NameError: l = t print " ".join(map(str,l)) @@ -193,11 +193,11 @@ if do_pseudo: l_str.append(l_dump) str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE" - print str_.format(a,i+1,zcore[i],len(l_str)) + print str_.format(a,i+1,int(zcore[i]),int(len(l_str)-1)) for i, l in enumerate(l_str): str_ = "FOR L= {0} COEFF N ZETA" - print str_.format(len(l_str)-i-1) + print str_.format(int(len(l_str)-i-1)) for ii, ll in enumerate(l): print " ",ii+1, ll From 99a9ba01e12fca7b328c7c315554d8e7432dfbf9 Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Wed, 27 Jan 2016 10:48:10 +0100 Subject: [PATCH 02/29] Same phase for the MOs than gamess --- .../qmcpack/qp_convert_qmcpack_from_ezfio.py | 356 +++++++++++------- 1 file changed, 225 insertions(+), 131 deletions(-) diff --git a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py index 010e277d..dc37d6ae 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py @@ -2,6 +2,11 @@ print "#QP -> QMCPACK" +# ___ +# | ._ o _|_ +# _|_ | | | |_ +# + from ezfio import ezfio import sys @@ -9,76 +14,99 @@ ezfio_path = sys.argv[1] ezfio.set_file(ezfio_path) - do_pseudo = ezfio.get_pseudo_do_pseudo() if do_pseudo: - print "do_pseudo True" - zcore = ezfio.get_pseudo_nucl_charge_remove() + print "do_pseudo True" + zcore = ezfio.get_pseudo_nucl_charge_remove() else: - print "do_pseudo False" + print "do_pseudo False" try: - n_det =ezfio.get_determinants_n_det() + n_det = ezfio.get_determinants_n_det() except IOError: - n_det = 1 + n_det = 1 if n_det == 1: - print "multi_det False" + print "multi_det False" else: - print "multi_det True" + print "multi_det True" +# +# |\/| o _ _ +# | | | _> (_ +# ao_num = ezfio.get_ao_basis_ao_num() print "ao_num", ao_num mo_num = ezfio.get_mo_basis_mo_tot_num() print "mo_num", mo_num - alpha = ezfio.get_electrons_elec_alpha_num() beta = ezfio.get_electrons_elec_beta_num() print "elec_alpha_num", alpha print "elec_beta_num", beta -print "elec_tot_num", alpha + beta -print "spin_multiplicity", 2*(alpha-beta)+1 +print "elec_tot_num", alpha + beta +print "spin_multiplicity", 2 * (alpha - beta) + 1 l_label = ezfio.get_nuclei_nucl_label() l_charge = ezfio.get_nuclei_nucl_charge() l_coord = ezfio.get_nuclei_nucl_coord() -l_coord_str = [" ".join(map(str,i)) for i in l_coord] +l_coord_str = [" ".join(map(str, i)) for i in l_coord] -print "nucl_num",len(l_label) +print "nucl_num", len(l_label) + +# _ +# / _ _ ._ _| +# \_ (_) (_) | (_| +# print "Atomic coord in Bohr" -for i,t in enumerate(zip(l_label,l_charge,l_coord_str)): - try : - l = (t[0],t[1]+zcore[i],t[2]) - except NameError: - l = t - print " ".join(map(str,l)) - +for i, t in enumerate(zip(l_label, l_charge, l_coord_str)): + try: + l = (t[0], t[1] + zcore[i], t[2]) + except NameError: + l = t + print " ".join(map(str, l)) +# +# Call externet process to get the sysmetry +# import subprocess -process = subprocess.Popen(['qp_print_basis', ezfio_path], stdout=subprocess.PIPE) +process = subprocess.Popen( + ['qp_print_basis', ezfio_path], + stdout=subprocess.PIPE) out, err = process.communicate() -basis_raw, sym_raw, mo_raw = out.split("\n\n\n") +basis_raw, sym_raw, _ = out.split("\n\n\n") + +# _ __ +# |_) _. _ o _ (_ _ _|_ +# |_) (_| _> | _> __) (/_ |_ +# + basis_without_header = "\n".join(basis_raw.split("\n")[7:]) -for i,l in enumerate(l_label): - basis_without_header=basis_without_header.replace('Atom {0}'.format(i+1),l) +for i, l in enumerate(l_label): + basis_without_header = basis_without_header.replace('Atom {0}'.format(i + 1), l) print "BEGIN_BASIS_SET" print "" print basis_without_header print "END_BASIS_SET" + # _ # |\/| / \ _ # | | \_/ _> -# +# + +# +# Function +# def same_character(item1): - return item1==item1[0]* len(item1) + return item1 == item1[0] * len(item1) + def compare_gamess_style(item1, item2): if len(item1) < len(item2): @@ -92,123 +120,190 @@ def compare_gamess_style(item1, item2): return 1 elif same_character(item1) and not same_character(item2): return -1 - elif not same_character(item1) and same_character(item2): + elif not same_character(item1) and same_character(item2): return 1 else: - return compare_gamess_style(item1[:-1],item2[:-1]) + return compare_gamess_style(item1[:-1], item2[:-1]) + def expend_and_order_sym(str_): - #Expend - for i,c in enumerate(str_): - try: - n = int(c) - except ValueError: - pass - else: - str_ = str_[:i-1] + str_[i-1]*n + str_[i+1:] + #Expend + for i, c in enumerate(str_): + try: + n = int(c) + except ValueError: + pass + else: + str_ = str_[:i - 1] + str_[i - 1] * n + str_[i + 1:] + + #Order by frequency + return "".join(sorted(str_, key=str_.count, reverse=True)) - #Order by frequency - return "".join(sorted(str_,key=str_.count,reverse=True)) def get_nb_permutation(str_): - l = len(str_)-1 - if l==0: - return 1 - else: - return 2*(2*l + 1) + l = len(str_) - 1 + if l == 0: + return 1 + else: + return 2 * (2 * l + 1) + +#========================== +# We will order the symetry +#========================== -## We will order the symetry l_sym_without_header = sym_raw.split("\n")[3:-2] l_l_sym = [i.split() for i in l_sym_without_header] for l in l_l_sym: - l[2] = expend_and_order_sym(l[2]) + l[2] = expend_and_order_sym(l[2]) l_l_sym_iter = iter(l_l_sym) -for i,l in enumerate(l_l_sym_iter): - n = get_nb_permutation(l[2]) - if n !=1: - l_l_sym[i:i+n] = sorted(l_l_sym[i:i+n],key=lambda x : x[2], cmp=compare_gamess_style) - for next_ in range(n-1): - next(l_l_sym_iter) - -#Is orderd now - -l_block = mo_raw.split("\n\n")[5:-1] +for i, l in enumerate(l_l_sym_iter): + n = get_nb_permutation(l[2]) + if n != 1: + l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n], + key=lambda x: x[2], + cmp=compare_gamess_style) + for next_ in range(n - 1): + next(l_l_sym_iter) -l_block_format=[] +#======== +#MO COEF +#======== +def order_phase(mo_coef): + #Order + mo_coef_phase = [] + import math -print "" -print "BEGIN_MO" -for block in l_block: - print "" - l_ligne = block.split("\n") - print l_ligne.pop(0) + for i in mo_coef: + if abs(max(i)) > abs(min(i)): + sign_max = math.copysign(1, max(i)) + else: + sign_max = math.copysign(1, min(i)) - for l in l_l_sym: - i = int(l[0]) - 1 - i_a = int(l[1]) - 1 - sym = l[2] + if sign_max == -1: + ii = [-1 * l for l in i] + else: + ii = i - print l_label[i_a],sym,l_ligne[i] + mo_coef_phase.append(ii) + return mo_coef_phase -print "END_MO" +def order_by_sim(mo_coef, l_l_sym): + l_sym_oder = [int(l[0]) - 1 for l in l_l_sym] + mo_coef_order = [[x for (y, x) in sorted(zip(l_sym_oder, i))] + for i in mo_coef] + return mo_coef_order + + +def chunked(mo_coef, chunks_size): + mo_coef_block = [] + for i in mo_coef: + chunks = [i[x:x + chunks_size] for x in xrange(0, len(i), chunks_size)] + mo_coef_block.append(chunks) + return mo_coef_block + + +def print_mo_coef(mo_coef_block, l_l_sym): + print "" + print "BEGIN_MO" + print "" + len_block_curent = 0 + nb_block = len(mo_coef_block[0]) + for i_block in range(0, nb_block): + a = [i[i_block] for i in mo_coef_block] + + print " ".join([str( + i + 1) for i in range(len_block_curent, len_block_curent + len(a[ + 0]))]) + + len_block_curent += len(a[0]) + + for l in l_l_sym: + i = int(l[0]) - 1 + i_a = int(l[1]) - 1 + sym = l[2] + + print l_label[i_a], sym, " ".join('{: 3.8f}'.format(i) + for i in a[i]) + + if i_block != nb_block - 1: + print "" + else: + print "END_MO" + + +mo_coef = ezfio.get_mo_basis_mo_coef() +mo_coef_phase = order_phase(mo_coef) +mo_coef_phase_order = order_by_sim(mo_coef_phase, l_l_sym) +mo_coef_transp = zip(*mo_coef_phase_order) +mo_coef_block = chunked(mo_coef_transp, 4) +print_mo_coef(mo_coef_block, l_l_sym) + +# _ +# |_) _ _ _| _ +# | _> (/_ |_| (_| (_) +# if do_pseudo: - print "" - print "BEGIN_PSEUDO" - klocmax = ezfio.get_pseudo_pseudo_klocmax() - kmax = ezfio.get_pseudo_pseudo_kmax() - lmax = ezfio.get_pseudo_pseudo_lmax() - - - n_k = ezfio.get_pseudo_pseudo_n_k() - v_k = ezfio.get_pseudo_pseudo_v_k() - dz_k = ezfio.get_pseudo_pseudo_dz_k() - - n_kl = ezfio.get_pseudo_pseudo_n_kl() - v_kl = ezfio.get_pseudo_pseudo_v_kl() - dz_kl = ezfio.get_pseudo_pseudo_dz_kl() - - def list_to_string(l): - return " ".join(map(str,l)) - - for i,a in enumerate(l_label): - - l_str = [] - - l_dump = [] - for k in range(klocmax): - if v_k[k][i]: - l_ = list_to_string([v_k[k][i], n_k[k][i]+2, dz_k[k][i]]) - l_dump.append(l_) - - l_str.append(l_dump) - for l in range(lmax+1): - l_dump = [] - for k in range(kmax): - if v_kl[l][k][i]: - l_ = list_to_string([v_kl[l][k][i], n_kl[l][k][i]+2, dz_kl[l][k][i]]) - l_dump.append(l_) - if l_dump: - l_str.append(l_dump) - - str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE" - print str_.format(a,i+1,int(zcore[i]),int(len(l_str)-1)) - - for i, l in enumerate(l_str): - str_ = "FOR L= {0} COEFF N ZETA" - print str_.format(int(len(l_str)-i-1)) - for ii, ll in enumerate(l): - print " ",ii+1, ll - - str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS." - print str_.format(sum(zcore)) - print "END_PSEUDO" - + print "" + print "BEGIN_PSEUDO" + klocmax = ezfio.get_pseudo_pseudo_klocmax() + kmax = ezfio.get_pseudo_pseudo_kmax() + lmax = ezfio.get_pseudo_pseudo_lmax() + + n_k = ezfio.get_pseudo_pseudo_n_k() + v_k = ezfio.get_pseudo_pseudo_v_k() + dz_k = ezfio.get_pseudo_pseudo_dz_k() + + n_kl = ezfio.get_pseudo_pseudo_n_kl() + v_kl = ezfio.get_pseudo_pseudo_v_kl() + dz_kl = ezfio.get_pseudo_pseudo_dz_kl() + + def list_to_string(l): + return " ".join(map(str, l)) + + for i, a in enumerate(l_label): + + l_str = [] + + l_dump = [] + for k in range(klocmax): + if v_k[k][i]: + l_ = list_to_string([v_k[k][i], n_k[k][i] + 2, dz_k[k][i]]) + l_dump.append(l_) + + l_str.append(l_dump) + for l in range(lmax + 1): + l_dump = [] + for k in range(kmax): + if v_kl[l][k][i]: + l_ = list_to_string([v_kl[l][k][i], n_kl[l][k][i] + 2, + dz_kl[l][k][i]]) + l_dump.append(l_) + if l_dump: + l_str.append(l_dump) + + str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE" + print str_.format(a, i + 1, int(zcore[i]), int(len(l_str) - 1)) + + for i, l in enumerate(l_str): + str_ = "FOR L= {0} COEFF N ZETA" + print str_.format(int(len(l_str) - i - 1)) + for ii, ll in enumerate(l): + print " ", ii + 1, ll + + str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS." + print str_.format(sum(zcore)) + print "END_PSEUDO" + +# _ +# | \ _ _|_ +# |_/ (/_ |_ +# print "" print "BEGIN_DET" print "" @@ -219,18 +314,17 @@ print "" psi_det = ezfio.get_determinants_psi_det() psi_coef = ezfio.get_determinants_psi_coef()[0] +for c, (l_det_bit_alpha, l_det_bit_beta) in zip(psi_coef, psi_det): + print c + for det in l_det_bit_alpha: + bin_det_raw = "{0:b}".format(det)[::-1] + bin_det = bin_det_raw + "0" * (mo_num - len(bin_det_raw)) + print bin_det -for c, (l_det_bit_alpha, l_det_bit_beta) in zip(psi_coef,psi_det): - print c - for det in l_det_bit_alpha: - bin_det_raw = "{0:b}".format(det)[::-1] - bin_det = bin_det_raw+"0"*(mo_num-len(bin_det_raw)) - print bin_det - - for det in l_det_bit_beta: - bin_det_raw = "{0:b}".format(det)[::-1] - bin_det = bin_det_raw+"0"*(mo_num-len(bin_det_raw)) - print bin_det - print "" + for det in l_det_bit_beta: + bin_det_raw = "{0:b}".format(det)[::-1] + bin_det = bin_det_raw + "0" * (mo_num - len(bin_det_raw)) + print bin_det + print "" print "END_DET" From 9d3900c7ee39df5fbfc16ca955589ea02b479611 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 27 Jan 2016 17:15:57 +0100 Subject: [PATCH 03/29] qp_module is now case insensitive --- configure | 2 +- scripts/module/qp_module.py | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/configure b/configure index a6c5bc2e..2344264a 100755 --- a/configure +++ b/configure @@ -533,7 +533,7 @@ def recommendation(): print " source {0}".format(path) print "" print "Then, install the modules you want to install using :" - print " qp_install_module.py " + print " qp_module.py" print "" print "Finally :" print " ninja" diff --git a/scripts/module/qp_module.py b/scripts/module/qp_module.py index 3eda69f3..63649df5 100755 --- a/scripts/module/qp_module.py +++ b/scripts/module/qp_module.py @@ -175,7 +175,11 @@ def main(arguments): d_child = d_local.copy() d_child.update(d_plugin) - l_name = arguments[""] + normalize_case = {} + for name in d_local.keys() + d_plugin.keys(): + normalize_case [ name.lower() ] = name + + l_name = [ normalize_case[name.lower()] for name in arguments[""] ] for name in l_name: if name in d_local: @@ -206,6 +210,7 @@ def main(arguments): raise print "[ OK ]" + print "" print "You can now compile as usual" print "`cd {0} ; ninja` for exemple".format(QP_ROOT) print " or --in developement mode-- you can cd in a directory and compile here" From ccc061fbeb0099f4abacbf0ee3836af6727cb66f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 28 Jan 2016 14:50:24 +0100 Subject: [PATCH 04/29] Introduced MR-CEPA. Not working yet --- plugins/MRCC_Utils/H_apply.irp.f | 22 ++ plugins/MRCC_Utils/mrcc_utils.irp.f | 4 + plugins/MRCC_Utils/mrcepa_dress.irp.f | 260 ++++++++++++++++++++++++ plugins/MRCC_Utils/mrcepa_general.irp.f | 97 +++++++++ 4 files changed, 383 insertions(+) create mode 100644 plugins/MRCC_Utils/mrcepa_dress.irp.f create mode 100644 plugins/MRCC_Utils/mrcepa_general.irp.f diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index 7314674c..1cafc8de 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -24,5 +24,27 @@ s.data["size_max"] = "3072" print s +s = H_apply("mrcepa") +s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref" +s.data["declarations"] += """ + integer, intent(in) :: Ndet_ref,Ndet_non_ref + double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(in) :: delta_ii_(Ndet_ref,*) +""" +s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" +s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" +s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" +s.data["decls_main"] += """ + integer, intent(in) :: Ndet_ref,Ndet_non_ref + double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(in) :: delta_ii_(Ndet_ref,*) +""" +s.data["finalization"] = "" +s.data["copy_buffer"] = "" +s.data["generate_psi_guess"] = "" +s.data["size_max"] = "3072" +# print s + + END_SHELL diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 8fe6c411..1e2f974d 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -26,7 +26,11 @@ ! TODO --- Test perturbatif ------ do k=1,N_states lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + ! TODO : i_h_psi peut sortir de la boucle? call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current) + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif tmp = psi_non_ref_coef(i,k)/ihpsi_current(k) i_pert = 0 ! Perturbation only if 1st order < 0.5 x second order diff --git a/plugins/MRCC_Utils/mrcepa_dress.irp.f b/plugins/MRCC_Utils/mrcepa_dress.irp.f new file mode 100644 index 00000000..9789b818 --- /dev/null +++ b/plugins/MRCC_Utils/mrcepa_dress.irp.f @@ -0,0 +1,260 @@ +use omp_lib +use bitmasks + +subroutine mrcepa_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet_ref, Ndet_non_ref + double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) + double precision, intent(inout) :: delta_ii_(Ndet_ref,*) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,l + integer :: degree_alpha(psi_det_size) + integer :: idx_alpha(0:psi_det_size) + logical :: good, fullMatch + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq, c_ref ,degree + + double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision, allocatable :: dIa_hia(:,:) + double precision :: haj, phase, phase2 + double precision :: f(N_states), ci_inv(N_states) + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + integer(bit_kind) :: tmp_det(Nint,2) + integer(bit_kind) :: tmp_det_0(Nint,2) + integer :: iint, ipos + integer :: i_state, i_sd, k_sd, l_sd, i_I, i_alpha + + integer(bit_kind),allocatable :: miniList(:,:,:) + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + integer,allocatable :: idx_miniList(:) + integer :: N_miniList, ni, leng + integer(bit_kind) :: isum + + double precision :: hia + integer, allocatable :: index_sorted(:) + + + leng = max(N_det_generators, N_det_non_ref) + allocate(miniList(Nint, 2, leng), idx_miniList(leng), index_sorted(N_det)) + + !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) + call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) + + if(fullMatch) then + return + end if + + + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + + allocate (dIa_hia(N_states,Ndet_non_ref)) + + ! |I> + + ! |alpha> + + if(N_tq > 0) then + call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint) + end if + + + do i_alpha=1,N_tq + ! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) + call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) + + integer, external :: get_index_in_psi_det_sorted_bit + index_sorted = huge(-1) + do j=1,idx_alpha(0) + idx_alpha(j) = idx_miniList(idx_alpha(j)) + index_sorted( get_index_in_psi_det_sorted_bit( psi_non_ref(1,1,idx_alpha(j)), N_int ) ) = idx_alpha(j) + end do + + ! |I> + do i_I=1,N_det_ref + ! Find triples and quadruple grand parents + call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint) + if (degree > 4) then + cycle + endif + + do i_state=1,N_states + dIa(i_state) = 0.d0 + enddo + + !TODO: MR + do i_sd=1,idx_alpha(0) + call get_excitation_degree(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),degree,Nint) + if (degree > 2) then + cycle + endif + call get_excitation(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_det_0 = 0_bit_kind + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos) + if (degree == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos) + endif + + call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(i_sd)),Nint,hia) + + ! |alpha> + do k_sd=1,idx_alpha(0) + call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint) + if (degree > 2) then + cycle + endif + + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_det = 0_bit_kind + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + if (degree == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + endif + + isum = 0_bit_kind + do iint = 1,N_int + isum = isum + iand(tmp_det(iint,1), tmp_det_0(iint,1)) & + + iand(tmp_det(iint,2), tmp_det_0(iint,2)) + enddo + + if (isum /= 0_bit_kind) then + cycle + endif + + ! + ! + call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk) + do i_state=1,N_states + dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + enddo + ! |l> = Exc(k -> alpha) |I> + call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do k=1,N_int + tmp_det(k,1) = psi_ref(k,1,i_I) + tmp_det(k,2) = psi_ref(k,2,i_I) + enddo + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + if (degree == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + endif + + ! + do i_state=1,N_states + dka(i_state) = 0.d0 + enddo + + +! l_sd = index_sorted( get_index_in_psi_det_sorted_bit( tmp_det, N_int ) ) +! call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,Nint) +! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),Nint,hIl) +! do i_state=1,N_states +! dka(i_state) = hIl * lambda_mrcc(i_state,l_sd) * phase * phase2 +! enddo + + do l_sd=1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint) + if (degree == 0) then + call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) + call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) + do i_state=1,N_states + dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + enddo + exit + endif + enddo + do i_state=1,N_states + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) + enddo + enddo + + do i_state=1,N_states + ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state) + enddo + + k_sd = idx_alpha(i_sd) + do i_state=1,N_states + dIa_hia(i_state,k_sd) = dIa(i_state) * hia + enddo + + call omp_set_lock( psi_ref_lock(i_I) ) + do i_state=1,N_states + delta_ij_(i_I,k_sd,i_state) += dIa_hia(i_state,k_sd) + + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + delta_ii_(i_I,i_state) -= dIa_hia(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state) + else + delta_ii_(i_I,i_state) = 0.d0 + endif + enddo + call omp_unset_lock( psi_ref_lock(i_I) ) + enddo + enddo + + enddo + deallocate (dIa_hia,index_sorted) + deallocate(miniList, idx_miniList) +end + + + + + + + + diff --git a/plugins/MRCC_Utils/mrcepa_general.irp.f b/plugins/MRCC_Utils/mrcepa_general.irp.f new file mode 100644 index 00000000..3479548b --- /dev/null +++ b/plugins/MRCC_Utils/mrcepa_general.irp.f @@ -0,0 +1,97 @@ +subroutine run_mrcepa + implicit none + call set_generators_bitmasks_as_holes_and_particles + call mrcepa_iterations +end + +subroutine mrcepa_iterations + implicit none + + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration,i_oscillations + double precision :: E_past(4) + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + j = 1 + i_oscillations = 0 + do while (delta_E > 1.d-7) + iteration += 1 + print *, '===========================' + print *, 'MRCEPA Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCEPA energy") + call diagonalize_ci_dressed + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + + E_past(j) = E_new + j +=1 + if(j>4)then + j=1 + endif + if(iteration > 4) then + if(delta_E > 1.d-10)then + if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then + print*,'OSCILLATIONS !!!' + oscillations = .True. + i_oscillations +=1 + lambda_mrcc_tmp = lambda_mrcc + endif + endif + endif + call save_wavefunction +! if (i_oscillations > 5) then +! exit +! endif + if (iteration > 200) then + exit + endif + print*,'------------' + print*,'VECTOR' + do i = 1, N_det_ref + print*,'' + print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) + print*,'delta_ii(i,1) = ',delta_ii(i,1) + enddo + print*,'------------' + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCEPA energy") + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + call save_wavefunction + +end + +subroutine set_generators_bitmasks_as_holes_and_particles + implicit none + integer :: i,k + do k = 1, N_generators_bitmask + do i = 1, N_int + ! Pure single part + generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha + generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha + generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta + generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta + + ! Double excitation + generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha + generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha + generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta + generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta + + generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha + generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha + generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta + generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta + + enddo + enddo + touch generators_bitmask + + + +end From 7453b47f142c3e6ab9d67220ea37001435467967 Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Tue, 2 Feb 2016 16:20:12 +0100 Subject: [PATCH 05/29] qmcpack No phase reorder and only one basis by atom --- .../qmcpack/qp_convert_qmcpack_from_ezfio.py | 115 +++++++++++------- 1 file changed, 71 insertions(+), 44 deletions(-) diff --git a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py index dc37d6ae..76386e36 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py @@ -5,7 +5,7 @@ print "#QP -> QMCPACK" # ___ # | ._ o _|_ # _|_ | | | |_ -# +# from ezfio import ezfio @@ -34,7 +34,12 @@ else: # # |\/| o _ _ # | | | _> (_ -# +# + +def list_to_string(l): + return " ".join(map(str, l)) + + ao_num = ezfio.get_ao_basis_ao_num() print "ao_num", ao_num @@ -52,14 +57,14 @@ l_label = ezfio.get_nuclei_nucl_label() l_charge = ezfio.get_nuclei_nucl_charge() l_coord = ezfio.get_nuclei_nucl_coord() -l_coord_str = [" ".join(map(str, i)) for i in l_coord] +l_coord_str = [list_to_string(i) for i in zip(*l_coord)] print "nucl_num", len(l_label) # _ # / _ _ ._ _| # \_ (_) (_) | (_| -# +# print "Atomic coord in Bohr" for i, t in enumerate(zip(l_label, l_charge, l_coord_str)): @@ -67,7 +72,7 @@ for i, t in enumerate(zip(l_label, l_charge, l_coord_str)): l = (t[0], t[1] + zcore[i], t[2]) except NameError: l = t - print " ".join(map(str, l)) + print list_to_string(l) # # Call externet process to get the sysmetry @@ -78,29 +83,41 @@ process = subprocess.Popen( stdout=subprocess.PIPE) out, err = process.communicate() -basis_raw, sym_raw, _ = out.split("\n\n\n") +basis_raw, sym_raw, _= out.split("\n\n\n") # _ __ # |_) _. _ o _ (_ _ _|_ # |_) (_| _> | _> __) (/_ |_ -# - +# basis_without_header = "\n".join(basis_raw.split("\n")[7:]) -for i, l in enumerate(l_label): - basis_without_header = basis_without_header.replace('Atom {0}'.format(i + 1), l) -print "BEGIN_BASIS_SET" -print "" -print basis_without_header +import re +l_basis_raw = re.split('\n\s*\n', basis_without_header) + +a_already_print = [] + +l_basis_clean = [] + + +for i, (a,b) in enumerate(zip(l_label,l_basis_raw)): + + if a not in a_already_print: + l_basis_clean.append(b.replace('Atom {0}'.format(i + 1), a)) + a_already_print.append(a) + else: + continue + +print "BEGIN_BASIS_SET\n" +print "\n\n".join(l_basis_clean) print "END_BASIS_SET" - # _ # |\/| / \ _ # | | \_/ _> # + # # Function # @@ -126,8 +143,9 @@ def compare_gamess_style(item1, item2): return compare_gamess_style(item1[:-1], item2[:-1]) -def expend_and_order_sym(str_): - #Expend +def expend_sym_str(str_): + #Expend x2 -> xx + # yx2 -> xxy for i, c in enumerate(str_): try: n = int(c) @@ -140,6 +158,13 @@ def expend_and_order_sym(str_): return "".join(sorted(str_, key=str_.count, reverse=True)) +def expend_sym_l(l_l_sym): + for l in l_l_sym: + l[2] = expend_sym_str(l[2]) + + return l_l_sym + + def get_nb_permutation(str_): l = len(str_) - 1 @@ -148,26 +173,29 @@ def get_nb_permutation(str_): else: return 2 * (2 * l + 1) + +def order_l_l_sym(l_l_sym): + l_l_sym_iter = iter(l_l_sym) + for i, l in enumerate(l_l_sym_iter): + n = get_nb_permutation(l[2]) + if n != 1: + l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n], + key=lambda x: x[2], + cmp=compare_gamess_style) + for next_ in range(n - 1): + next(l_l_sym_iter) + return l_l_sym + #========================== # We will order the symetry #========================== l_sym_without_header = sym_raw.split("\n")[3:-2] -l_l_sym = [i.split() for i in l_sym_without_header] +l_l_sym_raw = [i.split() for i in l_sym_without_header] +l_l_sym_expend_sym = expend_sym_l(l_l_sym_raw) -for l in l_l_sym: - l[2] = expend_and_order_sym(l[2]) - -l_l_sym_iter = iter(l_l_sym) -for i, l in enumerate(l_l_sym_iter): - n = get_nb_permutation(l[2]) - if n != 1: - l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n], - key=lambda x: x[2], - cmp=compare_gamess_style) - for next_ in range(n - 1): - next(l_l_sym_iter) +l_l_sym_ordered = order_l_l_sym(l_l_sym_expend_sym) #======== @@ -200,12 +228,12 @@ def order_by_sim(mo_coef, l_l_sym): return mo_coef_order -def chunked(mo_coef, chunks_size): - mo_coef_block = [] - for i in mo_coef: +def chunked(l, chunks_size): + l_block = [] + for i in l: chunks = [i[x:x + chunks_size] for x in xrange(0, len(i), chunks_size)] - mo_coef_block.append(chunks) - return mo_coef_block + l_block.append(chunks) + return l_block def print_mo_coef(mo_coef_block, l_l_sym): @@ -216,10 +244,8 @@ def print_mo_coef(mo_coef_block, l_l_sym): nb_block = len(mo_coef_block[0]) for i_block in range(0, nb_block): a = [i[i_block] for i in mo_coef_block] - - print " ".join([str( - i + 1) for i in range(len_block_curent, len_block_curent + len(a[ - 0]))]) + r_ = range(len_block_curent, len_block_curent + len(a[0])) + print " ".join([str(i + 1) for i in r_]) len_block_curent += len(a[0]) @@ -238,11 +264,12 @@ def print_mo_coef(mo_coef_block, l_l_sym): mo_coef = ezfio.get_mo_basis_mo_coef() -mo_coef_phase = order_phase(mo_coef) -mo_coef_phase_order = order_by_sim(mo_coef_phase, l_l_sym) +#mo_coef_phase = order_phase(mo_coef) +mo_coef_phase = mo_coef +mo_coef_phase_order = order_by_sim(mo_coef_phase, l_l_sym_ordered) mo_coef_transp = zip(*mo_coef_phase_order) mo_coef_block = chunked(mo_coef_transp, 4) -print_mo_coef(mo_coef_block, l_l_sym) +print_mo_coef(mo_coef_block, l_l_sym_ordered) # _ # |_) _ _ _| _ @@ -263,13 +290,11 @@ if do_pseudo: v_kl = ezfio.get_pseudo_pseudo_v_kl() dz_kl = ezfio.get_pseudo_pseudo_dz_kl() - def list_to_string(l): - return " ".join(map(str, l)) - for i, a in enumerate(l_label): l_str = [] + #Local l_dump = [] for k in range(klocmax): if v_k[k][i]: @@ -277,6 +302,8 @@ if do_pseudo: l_dump.append(l_) l_str.append(l_dump) + + #Non local for l in range(lmax + 1): l_dump = [] for k in range(kmax): From d94138cfedeede81283819a3f65de71ed9bd41ba Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Feb 2016 00:37:03 +0100 Subject: [PATCH 06/29] Corrected bug when writing pseudo to EZFIO --- ocaml/Qputils.ml | 20 +- ocaml/qp_create_ezfio_from_xyz.ml | 553 ++++++++++++++++-------------- 2 files changed, 308 insertions(+), 265 deletions(-) diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 7f840fde..36dc1102 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core.Std (* let rec transpose = function @@ -24,5 +24,21 @@ let input_to_sexp s = print_endline ("("^result^")"); "("^result^")" |> Sexp.of_string -;; + +let rmdir dirname = + let rec remove_one dir = + Sys.chdir dir; + Sys.readdir "." + |> Array.iter ~f:(fun x -> + match (Sys.is_directory x, Sys.is_file x) with + | (`Yes, _) -> remove_one x + | (_, `Yes) -> Sys.remove x + | _ -> failwith ("Unable to remove file "^x^".") + ); + Sys.chdir ".."; + Unix.rmdir dir + in + remove_one dirname + + diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index bb7fd847..a6820d88 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -1,6 +1,6 @@ -open Qputils;; -open Qptypes;; -open Core.Std;; +open Qputils +open Qptypes +open Core.Std let spec = let open Command.Spec in @@ -316,289 +316,316 @@ let run ?o b c d m p cart xyz_file = if Sys.file_exists_exn ezfio_file then failwith (ezfio_file^" already exists"); - (* Create EZFIO *) - Ezfio.set_file ezfio_file; + let write_file () = + (* Create EZFIO *) + Ezfio.set_file ezfio_file; - (* Write Pseudo *) - let pseudo = - List.map nuclei ~f:(fun x -> - match pseudo_channel x.Atom.element with - | Some channel -> Pseudo.read_element channel x.Atom.element - | None -> Pseudo.empty x.Atom.element - ) - in + (* Write Pseudo *) + let pseudo = + List.map nuclei ~f:(fun x -> + match pseudo_channel x.Atom.element with + | Some channel -> Pseudo.read_element channel x.Atom.element + | None -> Pseudo.empty x.Atom.element + ) + in - let molecule = - let n_elec_to_remove = - List.fold pseudo ~init:0 ~f:(fun accu x -> - accu + (Positive_int.to_int x.Pseudo.n_elec)) - in - { Molecule.elec_alpha = - (Elec_alpha_number.to_int molecule.Molecule.elec_alpha) - - n_elec_to_remove/2 - |> Elec_alpha_number.of_int; - Molecule.elec_beta = - (Elec_beta_number.to_int molecule.Molecule.elec_beta) - - (n_elec_to_remove - n_elec_to_remove/2) - |> Elec_beta_number.of_int; - Molecule.nuclei = - let charges = - List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec - |> Float.of_int) - |> Array.of_list + let molecule = + let n_elec_to_remove = + List.fold pseudo ~init:0 ~f:(fun accu x -> + accu + (Positive_int.to_int x.Pseudo.n_elec)) in - List.mapi molecule.Molecule.nuclei ~f:(fun i x -> - { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) - |> Charge.of_float } - ) - } - in - let nuclei = - molecule.Molecule.nuclei @ dummy - in - - - (* Write Electrons *) - Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int - molecule.Molecule.elec_alpha ) ; - Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int - molecule.Molecule.elec_beta ) ; - - (* Write Nuclei *) - let labels = - List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei - and charges = - List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei - and coords = - (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @ - (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @ - (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in - let nucl_num = (List.length labels) in - Ezfio.set_nuclei_nucl_num nucl_num ; - Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| nucl_num |] ~data:labels); - Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| nucl_num |] ~data:charges); - Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); - - - (* Write pseudopotential *) - let () = - match p with - | None -> Ezfio.set_pseudo_do_pseudo false - | _ -> Ezfio.set_pseudo_do_pseudo true - in - - let klocmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.length x.Pseudo.local + { Molecule.elec_alpha = + (Elec_alpha_number.to_int molecule.Molecule.elec_alpha) + - n_elec_to_remove/2 + |> Elec_alpha_number.of_int; + Molecule.elec_beta = + (Elec_beta_number.to_int molecule.Molecule.elec_beta) + - (n_elec_to_remove - n_elec_to_remove/2) + |> Elec_beta_number.of_int; + Molecule.nuclei = + let charges = + List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec + |> Float.of_int) + |> Array.of_list + in + List.mapi molecule.Molecule.nuclei ~f:(fun i x -> + { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) + |> Charge.of_float } + ) + } in - if (x > accu) then x - else accu - ) - and kmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.length x.Pseudo.non_local + let nuclei = + molecule.Molecule.nuclei @ dummy in - if (x > accu) then x - else accu - ) - and lmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) -> + + + (* Write Electrons *) + Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int + molecule.Molecule.elec_alpha ) ; + Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int + molecule.Molecule.elec_beta ) ; + + (* Write Nuclei *) + let labels = + List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei + and charges = + List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei + and coords = + (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @ + (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @ + (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in + let nucl_num = (List.length labels) in + Ezfio.set_nuclei_nucl_num nucl_num ; + Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:labels); + Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:charges); + Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); + + + (* Write pseudopotential *) + let () = + match p with + | None -> Ezfio.set_pseudo_do_pseudo false + | _ -> Ezfio.set_pseudo_do_pseudo true + in + + let klocmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> let x = - Positive_int.to_int x.Pseudo.Primitive_non_local.proj + List.length x.Pseudo.local + in + if (x > accu) then x + else accu + ) + and lmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> + let x = + List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) -> + let x = + Positive_int.to_int x.Pseudo.Primitive_non_local.proj + in + if (x > accu) then x + else accu + ) in if (x > accu) then x else accu ) in - if (x > accu) then x - else accu - ) - in - - let () = - Ezfio.set_pseudo_pseudo_klocmax klocmax; - Ezfio.set_pseudo_pseudo_kmax kmax; - Ezfio.set_pseudo_pseudo_lmax lmax; - let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k = - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0 - in - List.iteri pseudo ~f:(fun j x -> - List.iteri x.Pseudo.local ~f:(fun i (y,c) -> - tmp_array_v_k.(i).(j) <- AO_coef.to_float c; - let y, z = - AO_expo.to_float y.Pseudo.Primitive_local.expo, - R_power.to_int y.Pseudo.Primitive_local.r_power + + let kmax = + Array.init (lmax+1) ~f:(fun i-> + List.map pseudo ~f:(fun x -> + List.filter x.Pseudo.non_local ~f:(fun (y,_) -> + (Positive_int.to_int y.Pseudo.Primitive_non_local.proj) = i) + |> List.length ) + |> List.fold ~init:0 ~f:(fun accu x -> + if accu > x then accu else x) + ) + |> Array.fold ~init:0 ~f:(fun accu i -> + if i > accu then i else accu) + in + + + let () = + Ezfio.set_pseudo_pseudo_klocmax klocmax; + Ezfio.set_pseudo_pseudo_kmax kmax; + Ezfio.set_pseudo_pseudo_lmax lmax; + let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k = + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0 in - tmp_array_dz_k.(i).(j) <- y; - tmp_array_n_k.(i).(j) <- z; - ) - ); - let concat_2d tmp_array = - let data = - Array.map tmp_array ~f:Array.to_list - |> Array.to_list - |> List.concat - in - Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data - in - concat_2d tmp_array_v_k - |> Ezfio.set_pseudo_pseudo_v_k ; - concat_2d tmp_array_dz_k - |> Ezfio.set_pseudo_pseudo_dz_k; - concat_2d tmp_array_n_k - |> Ezfio.set_pseudo_pseudo_n_k; + List.iteri pseudo ~f:(fun j x -> + List.iteri x.Pseudo.local ~f:(fun i (y,c) -> + tmp_array_v_k.(i).(j) <- AO_coef.to_float c; + let y, z = + AO_expo.to_float y.Pseudo.Primitive_local.expo, + R_power.to_int y.Pseudo.Primitive_local.r_power + in + tmp_array_dz_k.(i).(j) <- y; + tmp_array_n_k.(i).(j) <- z; + ) + ); + let concat_2d tmp_array = + let data = + Array.map tmp_array ~f:Array.to_list + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data + in + concat_2d tmp_array_v_k + |> Ezfio.set_pseudo_pseudo_v_k ; + concat_2d tmp_array_dz_k + |> Ezfio.set_pseudo_pseudo_dz_k; + concat_2d tmp_array_n_k + |> Ezfio.set_pseudo_pseudo_n_k; - let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ), - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ), - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 ) - in - List.iteri pseudo ~f:(fun j x -> - List.iteri x.Pseudo.non_local ~f:(fun i (y,c) -> - let k, y, z = - Positive_int.to_int y.Pseudo.Primitive_non_local.proj, - AO_expo.to_float y.Pseudo.Primitive_non_local.expo, - R_power.to_int y.Pseudo.Primitive_non_local.r_power + let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )), + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )), + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 )) in - tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c; - tmp_array_dz_kl.(k).(i).(j) <- y; - tmp_array_n_kl.(k).(i).(j) <- z; - ) - ); - let concat_3d tmp_array = - let data = - Array.map tmp_array ~f:(fun x -> - Array.map x ~f:Array.to_list - |> Array.to_list - |> List.concat) - |> Array.to_list + List.iteri pseudo ~f:(fun j x -> + let last_idx = + Array.create ~len:(lmax+1) 0 + in + List.iter x.Pseudo.non_local ~f:(fun (y,c) -> + let k, y, z = + Positive_int.to_int y.Pseudo.Primitive_non_local.proj, + AO_expo.to_float y.Pseudo.Primitive_non_local.expo, + R_power.to_int y.Pseudo.Primitive_non_local.r_power + in + let i = + last_idx.(k) + in + tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c; + tmp_array_dz_kl.(k).(i).(j) <- y; + tmp_array_n_kl.(k).(i).(j) <- z; + last_idx.(k) <- i+1; + ) + ); + let concat_3d tmp_array = + let data = + Array.map tmp_array ~f:(fun x -> + Array.map x ~f:Array.to_list + |> Array.to_list + |> List.concat) + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data + in + concat_3d tmp_array_v_kl + |> Ezfio.set_pseudo_pseudo_v_kl ; + concat_3d tmp_array_dz_kl + |> Ezfio.set_pseudo_pseudo_dz_kl ; + concat_3d tmp_array_n_kl + |> Ezfio.set_pseudo_pseudo_n_kl ; + in + + + + + (* Write Basis set *) + let basis = + + let nmax = + Nucl_number.get_max () + in + let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function + | [] -> accu + | e::tail -> + let new_accu = + (e,(Nucl_number.of_int ~max:nmax n))::accu + in + do_work new_accu (n+1) tail + in + let result = do_work [] 1 nuclei + |> List.rev + |> List.map ~f:(fun (x,i) -> + try + let e = + match x.Atom.element with + | Element.X -> Element.H + | e -> e + in + Basis.read_element (basis_channel x.Atom.element) i e + with + | End_of_file -> failwith + ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.") + ) |> List.concat + in + (* close all in_channels *) + result in - Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data - in - concat_3d tmp_array_v_kl - |> Ezfio.set_pseudo_pseudo_v_kl ; - concat_3d tmp_array_dz_kl - |> Ezfio.set_pseudo_pseudo_dz_kl ; - concat_3d tmp_array_n_kl - |> Ezfio.set_pseudo_pseudo_n_kl ; - in - - - - - (* Write Basis set *) - let basis = - - let nmax = - Nucl_number.get_max () - in - let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function - | [] -> accu - | e::tail -> - let new_accu = - (e,(Nucl_number.of_int ~max:nmax n))::accu + let long_basis = Long_basis.of_basis basis in + let ao_num = List.length long_basis in + Ezfio.set_ao_basis_ao_num ao_num; + Ezfio.set_ao_basis_ao_basis b; + let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) + and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n) + and ao_power= + let l = List.map long_basis ~f:(fun (x,_,_) -> x) in + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) in - do_work new_accu (n+1) tail - in - let result = do_work [] 1 nuclei - |> List.rev - |> List.map ~f:(fun (x,i) -> - try - let e = - match x.Atom.element with - | Element.X -> Element.H - | e -> e - in - Basis.read_element (basis_channel x.Atom.element) i e - with - | End_of_file -> failwith - ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.") - ) - |> List.concat - in - (* close all in_channels *) - result - in - let long_basis = Long_basis.of_basis basis in - let ao_num = List.length long_basis in - Ezfio.set_ao_basis_ao_num ao_num; - Ezfio.set_ao_basis_ao_basis b; - let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) - and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n) - and ao_power= - let l = List.map long_basis ~f:(fun (x,_,_) -> x) in - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@ - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@ - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) - in - let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x -> - if x > s then x - else s) ao_prim_num - in - let gtos = - List.map long_basis ~f:(fun (_,x,_) -> x) - in - - let create_expo_coef ec = - let coefs = - begin match ec with - | `Coefs -> List.map gtos ~f:(fun x-> - List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) ) - | `Expos -> List.map gtos ~f:(fun x-> - List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float - prim.Primitive.expo) ) - end + let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x -> + if x > s then x + else s) ao_prim_num in - let rec get_n n accu = function - | [] -> List.rev accu - | h::tail -> - let y = - begin match List.nth h n with - | Some x -> x - | None -> 0. + let gtos = + List.map long_basis ~f:(fun (_,x,_) -> x) + in + + let create_expo_coef ec = + let coefs = + begin match ec with + | `Coefs -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) ) + | `Expos -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float + prim.Primitive.expo) ) end - in - get_n n (y::accu) tail + in + let rec get_n n accu = function + | [] -> List.rev accu + | h::tail -> + let y = + begin match List.nth h n with + | Some x -> x + | None -> 0. + end + in + get_n n (y::accu) tail + in + let rec build accu = function + | n when n=ao_prim_num_max -> accu + | n -> build ( accu @ (get_n n [] coefs) ) (n+1) + in + build [] 0 in - let rec build accu = function - | n when n=ao_prim_num_max -> accu - | n -> build ( accu @ (get_n n [] coefs) ) (n+1) - in - build [] 0 - in - let ao_coef = create_expo_coef `Coefs - and ao_expo = create_expo_coef `Expos + let ao_coef = create_expo_coef `Coefs + and ao_expo = create_expo_coef `Expos + in + let () = + Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; + Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; + Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; + Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; + Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; + Ezfio.set_ao_basis_ao_cartesian(cart); + in + match Input.Ao_basis.read () with + | None -> failwith "Error in basis" + | Some x -> Input.Ao_basis.write x in let () = - Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; - Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; - Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; - Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; - Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; - Ezfio.set_ao_basis_ao_cartesian(cart); - in - match Input.Ao_basis.read () with - | None -> failwith "Error in basis" - | Some x -> Input.Ao_basis.write x + try write_file () with + | ex -> + begin + begin + match Sys.is_directory ezfio_file with + | `Yes -> rmdir ezfio_file + | _ -> () + end; + raise ex; + end + in () From 0dbd2ad08b9394a1f55e52502a2b25ca132070ff Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Wed, 3 Feb 2016 17:46:14 +0100 Subject: [PATCH 07/29] Update qp_convert_qmcpack_from_ezfio.py --- .../qmcpack/qp_convert_qmcpack_from_ezfio.py | 403 ++++++++++++------ 1 file changed, 264 insertions(+), 139 deletions(-) diff --git a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py index 7443be68..76386e36 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py @@ -2,6 +2,11 @@ print "#QP -> QMCPACK" +# ___ +# | ._ o _|_ +# _|_ | | | |_ +# + from ezfio import ezfio import sys @@ -9,19 +14,31 @@ ezfio_path = sys.argv[1] ezfio.set_file(ezfio_path) - do_pseudo = ezfio.get_pseudo_do_pseudo() if do_pseudo: - print "do_pseudo True" - zcore = ezfio.get_pseudo_nucl_charge_remove() + print "do_pseudo True" + zcore = ezfio.get_pseudo_nucl_charge_remove() else: - print "do_pseudo False" + print "do_pseudo False" + +try: + n_det = ezfio.get_determinants_n_det() +except IOError: + n_det = 1 -n_det =ezfio.get_determinants_n_det() if n_det == 1: - print "multi_det False" + print "multi_det False" else: - print "multi_det True" + print "multi_det True" + +# +# |\/| o _ _ +# | | | _> (_ +# + +def list_to_string(l): + return " ".join(map(str, l)) + ao_num = ezfio.get_ao_basis_ao_num() print "ao_num", ao_num @@ -29,52 +46,84 @@ print "ao_num", ao_num mo_num = ezfio.get_mo_basis_mo_tot_num() print "mo_num", mo_num - alpha = ezfio.get_electrons_elec_alpha_num() beta = ezfio.get_electrons_elec_beta_num() print "elec_alpha_num", alpha print "elec_beta_num", beta -print "elec_tot_num", alpha + beta -print "spin_multiplicity", 2*(alpha-beta)+1 +print "elec_tot_num", alpha + beta +print "spin_multiplicity", 2 * (alpha - beta) + 1 l_label = ezfio.get_nuclei_nucl_label() l_charge = ezfio.get_nuclei_nucl_charge() l_coord = ezfio.get_nuclei_nucl_coord() -l_coord_str = [" ".join(map(str,i)) for i in l_coord] +l_coord_str = [list_to_string(i) for i in zip(*l_coord)] -print "nucl_num",len(l_label) +print "nucl_num", len(l_label) + +# _ +# / _ _ ._ _| +# \_ (_) (_) | (_| +# print "Atomic coord in Bohr" -for i,t in enumerate(zip(l_label,l_charge,l_coord_str)): - try : - l = (t[0],t[1]+zcore[i],t[1]) - except NameError: - l = t - print " ".join(map(str,l)) - +for i, t in enumerate(zip(l_label, l_charge, l_coord_str)): + try: + l = (t[0], t[1] + zcore[i], t[2]) + except NameError: + l = t + print list_to_string(l) +# +# Call externet process to get the sysmetry +# import subprocess -process = subprocess.Popen(['qp_print_basis', ezfio_path], stdout=subprocess.PIPE) +process = subprocess.Popen( + ['qp_print_basis', ezfio_path], + stdout=subprocess.PIPE) out, err = process.communicate() -basis_raw, sym_raw, mo_raw = out.split("\n\n\n") +basis_raw, sym_raw, _= out.split("\n\n\n") + +# _ __ +# |_) _. _ o _ (_ _ _|_ +# |_) (_| _> | _> __) (/_ |_ +# basis_without_header = "\n".join(basis_raw.split("\n")[7:]) -for i,l in enumerate(l_label): - basis_without_header=basis_without_header.replace('Atom {0}'.format(i+1),l) -print "BEGIN_BASIS_SET" -print "" -print basis_without_header +import re +l_basis_raw = re.split('\n\s*\n', basis_without_header) + +a_already_print = [] + +l_basis_clean = [] + + +for i, (a,b) in enumerate(zip(l_label,l_basis_raw)): + + if a not in a_already_print: + l_basis_clean.append(b.replace('Atom {0}'.format(i + 1), a)) + a_already_print.append(a) + else: + continue + +print "BEGIN_BASIS_SET\n" +print "\n\n".join(l_basis_clean) print "END_BASIS_SET" # _ # |\/| / \ _ # | | \_/ _> -# +# + + +# +# Function +# def same_character(item1): - return item1==item1[0]* len(item1) + return item1 == item1[0] * len(item1) + def compare_gamess_style(item1, item2): if len(item1) < len(item2): @@ -88,123 +137,200 @@ def compare_gamess_style(item1, item2): return 1 elif same_character(item1) and not same_character(item2): return -1 - elif not same_character(item1) and same_character(item2): + elif not same_character(item1) and same_character(item2): return 1 else: - return compare_gamess_style(item1[:-1],item2[:-1]) + return compare_gamess_style(item1[:-1], item2[:-1]) -def expend_and_order_sym(str_): - #Expend - for i,c in enumerate(str_): - try: - n = int(c) - except ValueError: - pass - else: - str_ = str_[:i-1] + str_[i-1]*n + str_[i+1:] - #Order by frequency - return "".join(sorted(str_,key=str_.count,reverse=True)) +def expend_sym_str(str_): + #Expend x2 -> xx + # yx2 -> xxy + for i, c in enumerate(str_): + try: + n = int(c) + except ValueError: + pass + else: + str_ = str_[:i - 1] + str_[i - 1] * n + str_[i + 1:] + + #Order by frequency + return "".join(sorted(str_, key=str_.count, reverse=True)) + + +def expend_sym_l(l_l_sym): + for l in l_l_sym: + l[2] = expend_sym_str(l[2]) + + return l_l_sym + def get_nb_permutation(str_): - l = len(str_)-1 - if l==0: - return 1 - else: - return 2*(2*l + 1) + l = len(str_) - 1 + if l == 0: + return 1 + else: + return 2 * (2 * l + 1) + + +def order_l_l_sym(l_l_sym): + l_l_sym_iter = iter(l_l_sym) + for i, l in enumerate(l_l_sym_iter): + n = get_nb_permutation(l[2]) + if n != 1: + l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n], + key=lambda x: x[2], + cmp=compare_gamess_style) + for next_ in range(n - 1): + next(l_l_sym_iter) + return l_l_sym + +#========================== +# We will order the symetry +#========================== -## We will order the symetry l_sym_without_header = sym_raw.split("\n")[3:-2] -l_l_sym = [i.split() for i in l_sym_without_header] +l_l_sym_raw = [i.split() for i in l_sym_without_header] +l_l_sym_expend_sym = expend_sym_l(l_l_sym_raw) -for l in l_l_sym: - l[2] = expend_and_order_sym(l[2]) - -l_l_sym_iter = iter(l_l_sym) -for i,l in enumerate(l_l_sym_iter): - n = get_nb_permutation(l[2]) - if n !=1: - l_l_sym[i:i+n] = sorted(l_l_sym[i:i+n],key=lambda x : x[2], cmp=compare_gamess_style) - for next_ in range(n-1): - next(l_l_sym_iter) - -#Is orderd now - -l_block = mo_raw.split("\n\n")[5:-1] +l_l_sym_ordered = order_l_l_sym(l_l_sym_expend_sym) -l_block_format=[] +#======== +#MO COEF +#======== +def order_phase(mo_coef): + #Order + mo_coef_phase = [] + import math -print "" -print "BEGIN_MO" -for block in l_block: - print "" - l_ligne = block.split("\n") - print l_ligne.pop(0) + for i in mo_coef: + if abs(max(i)) > abs(min(i)): + sign_max = math.copysign(1, max(i)) + else: + sign_max = math.copysign(1, min(i)) - for l in l_l_sym: - i = int(l[0]) - 1 - i_a = int(l[1]) - 1 - sym = l[2] + if sign_max == -1: + ii = [-1 * l for l in i] + else: + ii = i - print l_label[i_a],sym,l_ligne[i] + mo_coef_phase.append(ii) + return mo_coef_phase -print "END_MO" +def order_by_sim(mo_coef, l_l_sym): + l_sym_oder = [int(l[0]) - 1 for l in l_l_sym] + mo_coef_order = [[x for (y, x) in sorted(zip(l_sym_oder, i))] + for i in mo_coef] + return mo_coef_order + + +def chunked(l, chunks_size): + l_block = [] + for i in l: + chunks = [i[x:x + chunks_size] for x in xrange(0, len(i), chunks_size)] + l_block.append(chunks) + return l_block + + +def print_mo_coef(mo_coef_block, l_l_sym): + print "" + print "BEGIN_MO" + print "" + len_block_curent = 0 + nb_block = len(mo_coef_block[0]) + for i_block in range(0, nb_block): + a = [i[i_block] for i in mo_coef_block] + r_ = range(len_block_curent, len_block_curent + len(a[0])) + print " ".join([str(i + 1) for i in r_]) + + len_block_curent += len(a[0]) + + for l in l_l_sym: + i = int(l[0]) - 1 + i_a = int(l[1]) - 1 + sym = l[2] + + print l_label[i_a], sym, " ".join('{: 3.8f}'.format(i) + for i in a[i]) + + if i_block != nb_block - 1: + print "" + else: + print "END_MO" + + +mo_coef = ezfio.get_mo_basis_mo_coef() +#mo_coef_phase = order_phase(mo_coef) +mo_coef_phase = mo_coef +mo_coef_phase_order = order_by_sim(mo_coef_phase, l_l_sym_ordered) +mo_coef_transp = zip(*mo_coef_phase_order) +mo_coef_block = chunked(mo_coef_transp, 4) +print_mo_coef(mo_coef_block, l_l_sym_ordered) + +# _ +# |_) _ _ _| _ +# | _> (/_ |_| (_| (_) +# if do_pseudo: - print "" - print "BEGIN_PSEUDO" - klocmax = ezfio.get_pseudo_pseudo_klocmax() - kmax = ezfio.get_pseudo_pseudo_kmax() - lmax = ezfio.get_pseudo_pseudo_lmax() - - - n_k = ezfio.get_pseudo_pseudo_n_k() - v_k = ezfio.get_pseudo_pseudo_v_k() - dz_k = ezfio.get_pseudo_pseudo_dz_k() - - n_kl = ezfio.get_pseudo_pseudo_n_kl() - v_kl = ezfio.get_pseudo_pseudo_v_kl() - dz_kl = ezfio.get_pseudo_pseudo_dz_kl() - - def list_to_string(l): - return " ".join(map(str,l)) - - for i,a in enumerate(l_label): - - l_str = [] - - l_dump = [] - for k in range(klocmax): - if v_k[k][i]: - l_ = list_to_string([v_k[k][i], n_k[k][i]+2, dz_k[k][i]]) - l_dump.append(l_) - - l_str.append(l_dump) - for l in range(lmax+1): - l_dump = [] - for k in range(kmax): - if v_kl[l][k][i]: - l_ = list_to_string([v_kl[l][k][i], n_kl[l][k][i]+2, dz_kl[l][k][i]]) - l_dump.append(l_) - if l_dump: - l_str.append(l_dump) - - str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE" - print str_.format(a,i+1,zcore[i],len(l_str)) - - for i, l in enumerate(l_str): - str_ = "FOR L= {0} COEFF N ZETA" - print str_.format(len(l_str)-i-1) - for ii, ll in enumerate(l): - print " ",ii+1, ll - - str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS." - print str_.format(sum(zcore)) - print "END_PSEUDO" - + print "" + print "BEGIN_PSEUDO" + klocmax = ezfio.get_pseudo_pseudo_klocmax() + kmax = ezfio.get_pseudo_pseudo_kmax() + lmax = ezfio.get_pseudo_pseudo_lmax() + + n_k = ezfio.get_pseudo_pseudo_n_k() + v_k = ezfio.get_pseudo_pseudo_v_k() + dz_k = ezfio.get_pseudo_pseudo_dz_k() + + n_kl = ezfio.get_pseudo_pseudo_n_kl() + v_kl = ezfio.get_pseudo_pseudo_v_kl() + dz_kl = ezfio.get_pseudo_pseudo_dz_kl() + + for i, a in enumerate(l_label): + + l_str = [] + + #Local + l_dump = [] + for k in range(klocmax): + if v_k[k][i]: + l_ = list_to_string([v_k[k][i], n_k[k][i] + 2, dz_k[k][i]]) + l_dump.append(l_) + + l_str.append(l_dump) + + #Non local + for l in range(lmax + 1): + l_dump = [] + for k in range(kmax): + if v_kl[l][k][i]: + l_ = list_to_string([v_kl[l][k][i], n_kl[l][k][i] + 2, + dz_kl[l][k][i]]) + l_dump.append(l_) + if l_dump: + l_str.append(l_dump) + + str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE" + print str_.format(a, i + 1, int(zcore[i]), int(len(l_str) - 1)) + + for i, l in enumerate(l_str): + str_ = "FOR L= {0} COEFF N ZETA" + print str_.format(int(len(l_str) - i - 1)) + for ii, ll in enumerate(l): + print " ", ii + 1, ll + + str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS." + print str_.format(sum(zcore)) + print "END_PSEUDO" + +# _ +# | \ _ _|_ +# |_/ (/_ |_ +# print "" print "BEGIN_DET" print "" @@ -215,18 +341,17 @@ print "" psi_det = ezfio.get_determinants_psi_det() psi_coef = ezfio.get_determinants_psi_coef()[0] +for c, (l_det_bit_alpha, l_det_bit_beta) in zip(psi_coef, psi_det): + print c + for det in l_det_bit_alpha: + bin_det_raw = "{0:b}".format(det)[::-1] + bin_det = bin_det_raw + "0" * (mo_num - len(bin_det_raw)) + print bin_det -for c, (l_det_bit_alpha, l_det_bit_beta) in zip(psi_coef,psi_det): - print c - for det in l_det_bit_alpha: - bin_det_raw = "{0:b}".format(det)[::-1] - bin_det = bin_det_raw+"0"*(mo_num-len(bin_det_raw)) - print bin_det + for det in l_det_bit_beta: + bin_det_raw = "{0:b}".format(det)[::-1] + bin_det = bin_det_raw + "0" * (mo_num - len(bin_det_raw)) + print bin_det + print "" - for det in l_det_bit_beta: - bin_det_raw = "{0:b}".format(det)[::-1] - bin_det = bin_det_raw+"0"*(mo_num-len(bin_det_raw)) - print bin_det - print "" - -print "END_DET" \ No newline at end of file +print "END_DET" From e08bc59742c9bc68f0fb264b28c95aefcd7f26f8 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Wed, 3 Feb 2016 17:56:48 +0100 Subject: [PATCH 08/29] Update README.md --- README.md | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index f8b7e308..a21189ba 100644 --- a/README.md +++ b/README.md @@ -56,10 +56,14 @@ This file contains all the environment variables needed by the quantum package b ### Optional) Add some modules - Usage: qp_module.py list (--installed|--avalaible-local|--avalaible-remote) - qp_module.py install ... - qp_module.py create -n [...] +``` +Usage: + qp_module.py create -n [...] qp_module.py download -n [...] + qp_module.py install ... + qp_module.py list (--installed | --available-local) + qp_module.py uninstall +``` For exemple you can type : `qp_module.py install Full_CI` From d93debe2dba57addefb2a420799112edb79e2cba Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Wed, 3 Feb 2016 18:01:28 +0100 Subject: [PATCH 09/29] Update configure --- configure | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure b/configure index 2344264a..b5c46269 100755 --- a/configure +++ b/configure @@ -539,7 +539,7 @@ def recommendation(): print " ninja" print " make -C ocaml" print "" - print "You can install more plugin with the qp_install command" + print "You can install more plugin with the qp_module.py install command" print "PS : For more info on compiling the code, read the README.md" From 6e6a8ac82ac3294ec255dfaca08cfece8bc5e5d5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 Feb 2016 10:56:39 +0100 Subject: [PATCH 10/29] Corrected bug in spherical guess MOs --- src/Determinants/connected_to_ref.irp.f | 5 +- src/MO_Basis/ao_ortho_canonical.irp.f | 150 ++++++++++++++++++------ src/Utils/LinearAlgebra.irp.f | 1 - 3 files changed, 116 insertions(+), 40 deletions(-) diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index dc7698b5..b93b18b6 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -95,12 +95,13 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) exit endif enddo - i += 1 - if (i > N_det) then + if (i >= N_det) then return endif + i += 1 + !DIR$ FORCEINLINE do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref) if ( (key(1,1) /= psi_det_sorted_bit(1,1,i)).or. & diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index 62a72149..1916bfe4 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -1,3 +1,87 @@ + BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num_align,ao_num)] +&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ] + implicit none + BEGIN_DOC +! matrix of the coefficients of the mos generated by the +! orthonormalization by the S^{-1/2} canonical transformation of the aos +! ao_cart_to_sphe_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital + END_DOC + integer :: i + integer, external :: ao_power_index + integer :: ibegin,j,k + ao_cart_to_sphe_coef(:,:) = 0.d0 + ! Assume order provided by ao_power_index + i = 1 + ao_cart_to_sphe_num = 0 + do while (i <= ao_num) + select case ( ao_l(i) ) + case (0) + ao_cart_to_sphe_coef(i,i) = 1.d0 + i += 1 + ao_cart_to_sphe_num += 1 + BEGIN_TEMPLATE + case ($SHELL) + if (ao_power(i,1) == $SHELL) then + do k=1,size(cart_to_sphe_$SHELL,2) + do j=1,size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k) + enddo + enddo + i += size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_num += size(cart_to_sphe_$SHELL,2) + endif + SUBST [ SHELL ] + 1;; + 2;; + 3;; + 4;; + 5;; + 6;; + 7;; + 8;; + 9;; + END_TEMPLATE + case default + stop 'Error in ao_cart_to_sphe' + end select + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num,ao_cart_to_sphe_num) ] + implicit none + BEGIN_DOC + ! AO overlap matrix in the spherical basis set + END_DOC + double precision, allocatable :: S(:,:) + allocate (S(ao_cart_to_sphe_num,ao_num)) + + call dgemm('T','N',ao_cart_to_sphe_num,ao_num,ao_num, 1.d0, & + ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), & + ao_overlap,size(ao_overlap,1), 0.d0, & + S, size(S,1)) + + call dgemm('N','N',ao_cart_to_sphe_num,ao_cart_to_sphe_num,ao_num, 1.d0, & + S, size(S,1), & + ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), 0.d0, & + ao_cart_to_sphe_overlap,size(ao_cart_to_sphe_overlap,1)) + + deallocate(S) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_num) ] + implicit none + BEGIN_DOC + ! AO_cart_to_sphe_coef^(-1) + END_DOC + + call get_pseudo_inverse(ao_cart_to_sphe_coef,ao_num,ao_cart_to_sphe_num, & + ao_cart_to_sphe_inv, size(ao_cart_to_sphe_coef,1)) +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)] &BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ] implicit none @@ -8,47 +92,39 @@ END_DOC integer :: i ao_ortho_canonical_coef(:,:) = 0.d0 + do i=1,ao_num + ao_ortho_canonical_coef(i,i) = 1.d0 + enddo + if (ao_cartesian) then - do i=1,ao_num - ao_ortho_canonical_coef(i,i) = 1.d0 - enddo + + ao_ortho_canonical_num = ao_num + call ortho_canonical(ao_overlap,size(ao_overlap,1), & + ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), & + ao_ortho_canonical_num) + + else - integer, external :: ao_power_index - integer :: ibegin,j,k - ! Assume order provided by ao_power_index - i = 1 - do while (i <= ao_num) - select case ( ao_l(i) ) - case (0) - ao_ortho_canonical_coef(i,i) = 1.d0 - i += 1 - BEGIN_TEMPLATE - case ($SHELL) - if (ao_power(i,1) == $SHELL) then - do k=1,size(cart_to_sphe_$SHELL,2) - do j=1,size(cart_to_sphe_$SHELL,1) - ao_ortho_canonical_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k) - enddo - enddo - i += size(cart_to_sphe_$SHELL,1) - endif - SUBST [ SHELL ] - 1;; - 2;; - 3;; - 4;; - 5;; - 6;; - 7;; - 8;; - 9;; - END_TEMPLATE - case default - stop 'Error in sphe_to_cart' - end select + + double precision, allocatable :: S(:,:) + + allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num)) + S = 0.d0 + do i=1,ao_cart_to_sphe_num + S(i,i) = 1.d0 enddo + + ao_ortho_canonical_num = ao_cart_to_sphe_num + call ortho_canonical (ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & + ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num) + + call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, & + ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & + S, size(S,1), & + 0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1)) + + deallocate(S) endif - call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,ao_ortho_canonical_num) END_PROVIDER BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)] diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index bcfd43b8..dfa5d982 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -88,7 +88,6 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) endif enddo do i=m+1,n - print *, D(i) D(i) = 0.d0 enddo From ceb4ef8d08b027084870e56f1ee45bfb70248798 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 Feb 2016 12:43:24 +0100 Subject: [PATCH 11/29] Updated tests --- tests/bats/qp.bats | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index ab686d71..71604759 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -100,11 +100,11 @@ function run_FCI() { } @test "SCF H2O cc-pVDZ" { - run_HF h2o.ezfio -76.0270218692377 + run_HF h2o.ezfio -76.0230273511649 } @test "FCI H2O cc-pVDZ" { - run_FCI h2o.ezfio 10000 -76.2382562245107 -76.2433933430971 + run_FCI h2o.ezfio 10000 -76.2315960087713 -76.237428352199 } @test "CAS_SD H2O cc-pVDZ" { @@ -113,10 +113,10 @@ function run_FCI() { ezfio set_file $INPUT ezfio set perturbation do_pt2_end False ezfio set determinants n_det_max 1000 - qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-25]" + qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-23]" qp_run cas_sd_selected $INPUT energy="$(ezfio get cas_sd energy)" - eq $energy -76.2219816183745 1.E-6 + eq $energy -76.2085511296300 1.E-6 } @test "MRCC H2O cc-pVDZ" { @@ -128,7 +128,8 @@ function run_FCI() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2313853324571 1.E-3 + eq $energy -76.2165731870755 1.E-3 + } @@ -138,11 +139,11 @@ function run_FCI() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -16.9483703904632 + run_HF h2o_pseudo.ezfio -16.9457263818675 } @test "FCI H2O VDZ pseudo" { - run_FCI h2o_pseudo.ezfio 2000 -17.1550015533975 -17.1645044211228 + run_FCI h2o_pseudo.ezfio 2000 -17.1476897854369 -17.1598005211929 } #=== Convert From 33d70af2e5493e2699b4582a8fb2cd70ce285d4d Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Thu, 4 Feb 2016 18:50:40 +0100 Subject: [PATCH 12/29] Correct the qmcpack script --- ocaml/qp_create_ezfio_from_xyz.ml | 2 +- plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index a6820d88..710523e4 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -85,7 +85,7 @@ let list_basis () = ) in List.sort basis_list ~cmp:String.ascending - |> String.concat ~sep:"\t" + |> String.concat ~sep:"\n" (** Run the program *) diff --git a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py index 76386e36..6555b3ab 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py @@ -90,7 +90,7 @@ basis_raw, sym_raw, _= out.split("\n\n\n") # |_) (_| _> | _> __) (/_ |_ # -basis_without_header = "\n".join(basis_raw.split("\n")[7:]) +basis_without_header = "\n".join(basis_raw.split("\n")[11:]) import re l_basis_raw = re.split('\n\s*\n', basis_without_header) From 647d4ad6de609f0cb90a68b24a975d4f2f1f004d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 Feb 2016 23:33:40 +0100 Subject: [PATCH 13/29] Corrected bug with sperical guess --- src/MO_Basis/ao_ortho_canonical.irp.f | 6 +++--- tests/bats/qp.bats | 14 +++++++------- tests/run_tests.sh | 6 ++++++ 3 files changed, 16 insertions(+), 10 deletions(-) create mode 100755 tests/run_tests.sh diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index 1916bfe4..95a771b0 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -16,15 +16,15 @@ do while (i <= ao_num) select case ( ao_l(i) ) case (0) - ao_cart_to_sphe_coef(i,i) = 1.d0 - i += 1 ao_cart_to_sphe_num += 1 + ao_cart_to_sphe_coef(i,ao_cart_to_sphe_num) = 1.d0 + i += 1 BEGIN_TEMPLATE case ($SHELL) if (ao_power(i,1) == $SHELL) then do k=1,size(cart_to_sphe_$SHELL,2) do j=1,size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k) + ao_cart_to_sphe_coef(i+j-1,ao_cart_to_sphe_num+k) = cart_to_sphe_$SHELL(j,k) enddo enddo i += size(cart_to_sphe_$SHELL,1) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 71604759..ae385fe8 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -100,11 +100,11 @@ function run_FCI() { } @test "SCF H2O cc-pVDZ" { - run_HF h2o.ezfio -76.0230273511649 + run_HF h2o.ezfio -0.760270218692179E+02 } @test "FCI H2O cc-pVDZ" { - run_FCI h2o.ezfio 10000 -76.2315960087713 -76.237428352199 + run_FCI h2o.ezfio 10000 -0.762382562429778E+02 -0.762433933485226E+02 } @test "CAS_SD H2O cc-pVDZ" { @@ -113,10 +113,10 @@ function run_FCI() { ezfio set_file $INPUT ezfio set perturbation do_pt2_end False ezfio set determinants n_det_max 1000 - qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-23]" + qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" qp_run cas_sd_selected $INPUT energy="$(ezfio get cas_sd energy)" - eq $energy -76.2085511296300 1.E-6 + eq $energy -0.762219854008117E+02 1.E-5 } @test "MRCC H2O cc-pVDZ" { @@ -128,7 +128,7 @@ function run_FCI() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2165731870755 1.E-3 + eq $energy -0.762303253805911E+02 1.E-3 } @@ -139,11 +139,11 @@ function run_FCI() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -16.9457263818675 + run_HF h2o_pseudo.ezfio -0.169483703904991E+02 } @test "FCI H2O VDZ pseudo" { - run_FCI h2o_pseudo.ezfio 2000 -17.1476897854369 -17.1598005211929 + run_FCI h2o_pseudo.ezfio 2000 -0.171550015498807E+02 -0.171645044185009E+02 } #=== Convert diff --git a/tests/run_tests.sh b/tests/run_tests.sh new file mode 100755 index 00000000..96f299f9 --- /dev/null +++ b/tests/run_tests.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +rm -rf work +exec bats bats/qp.bats + + From b1944d91d1c11612e331421d099cf828701d3600 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 6 Feb 2016 18:45:11 +0100 Subject: [PATCH 14/29] Added Youtube demo --- README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.md b/README.md index a21189ba..4f9f6719 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,12 @@ Set of quantum chemistry programs and libraries. For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>), or below for the installation instructions. +Demo +==== + +[![Full-CI energy of C2 in 2 minutes](http://img.youtube.com/vi/b89MpnQGDVo/0.jpg)](http://www.youtube.com/watch?v=b89MpnQGDVo "Quantum Package Demo") + + # Installation From 0cbc91667f2c07047378a1eeddfa94aad1d63dfc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 7 Feb 2016 22:00:40 +0100 Subject: [PATCH 15/29] README : Migrated video from YouTube to Vimeo --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 4f9f6719..b378962a 100644 --- a/README.md +++ b/README.md @@ -14,7 +14,7 @@ For more information, you can visit the [wiki of the project](http://github.com/ Demo ==== -[![Full-CI energy of C2 in 2 minutes](http://img.youtube.com/vi/b89MpnQGDVo/0.jpg)](http://www.youtube.com/watch?v=b89MpnQGDVo "Quantum Package Demo") +[![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_130x73.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo") # Installation From 4205cadd35433bf1a2ee9eb36af54b1b4b8bf992 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 7 Feb 2016 22:01:50 +0100 Subject: [PATCH 16/29] README : larger image --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index b378962a..ce90ee61 100644 --- a/README.md +++ b/README.md @@ -14,7 +14,7 @@ For more information, you can visit the [wiki of the project](http://github.com/ Demo ==== -[![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_130x73.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo") +[![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_295x166.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo") # Installation From 5d2434e19c2a5fb15f4fed19790819cc453ad6c7 Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Mon, 8 Feb 2016 09:03:00 +0100 Subject: [PATCH 17/29] qp_convert_qmcpack_from_ezfio -> qp_convert_qmcpack_to_ezfio --- ...ert_qmcpack_from_ezfio.py => qp_convert_qmcpack_to_ezfio.py} | 0 plugins/qmcpack/save_for_qmcpack.irp.f | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) rename plugins/qmcpack/{qp_convert_qmcpack_from_ezfio.py => qp_convert_qmcpack_to_ezfio.py} (100%) diff --git a/plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py similarity index 100% rename from plugins/qmcpack/qp_convert_qmcpack_from_ezfio.py rename to plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py diff --git a/plugins/qmcpack/save_for_qmcpack.irp.f b/plugins/qmcpack/save_for_qmcpack.irp.f index b66fdf45..95e3320c 100644 --- a/plugins/qmcpack/save_for_qmcpack.irp.f +++ b/plugins/qmcpack/save_for_qmcpack.irp.f @@ -15,6 +15,6 @@ program qmcpack enddo call ezfio_set_ao_basis_ao_coef(ao_coef) call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5') - call system('$QP_ROOT/src/qmcpack/qp_convert_qmcpack_from_ezfio.py '//trim(ezfio_filename)) + call system('$QP_ROOT/src/qmcpack/qp_convert_qmcpack_to_ezfio.py '//trim(ezfio_filename)) end From ad21a3b8cea737a867ef91a9dbd8c7d67bde7ff6 Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Mon, 8 Feb 2016 09:49:40 +0100 Subject: [PATCH 18/29] Solve OpenMP segfault lock init --- src/Utils/map_module.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index 24f5a328..5c883c9b 100644 --- a/src/Utils/map_module.f90 +++ b/src/Utils/map_module.f90 @@ -69,7 +69,6 @@ subroutine cache_map_init(map,sze) implicit none type (cache_map_type), intent(inout) :: map integer(cache_map_size_kind) :: sze - call omp_init_lock(map%lock) call omp_set_lock(map%lock) map%n_elements = 0_8 map%map_size = 0_8 @@ -101,6 +100,9 @@ subroutine map_init(map,keymax) stop 5 endif sze = 2 + do i=0_8,map%map_size + call omp_init_lock(map%map(i)%lock) + enddo !$OMP PARALLEL DEFAULT(NONE) SHARED(map,sze) PRIVATE(i) !$OMP DO SCHEDULE(STATIC,512) do i=0_8,map%map_size From 8c5338a081dd8a8599513a4c0c7fb81a1a191176 Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Mon, 8 Feb 2016 11:10:24 +0100 Subject: [PATCH 19/29] Working MO git add plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py! --- .../qmcpack/qp_convert_qmcpack_to_ezfio.py | 35 ++++++++----------- 1 file changed, 14 insertions(+), 21 deletions(-) diff --git a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py index 6555b3ab..a4a15745 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py @@ -175,29 +175,31 @@ def get_nb_permutation(str_): def order_l_l_sym(l_l_sym): - l_l_sym_iter = iter(l_l_sym) - for i, l in enumerate(l_l_sym_iter): - n = get_nb_permutation(l[2]) + n = 1 + for i in range(len(l_l_sym)): if n != 1: - l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n], - key=lambda x: x[2], - cmp=compare_gamess_style) - for next_ in range(n - 1): - next(l_l_sym_iter) + n += -1 + continue + + l = l_l_sym[i] + n = get_nb_permutation(l[2]) + + l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n], + key=lambda x: x[2], + cmp=compare_gamess_style) + return l_l_sym + #========================== # We will order the symetry #========================== l_sym_without_header = sym_raw.split("\n")[3:-2] - l_l_sym_raw = [i.split() for i in l_sym_without_header] l_l_sym_expend_sym = expend_sym_l(l_l_sym_raw) - l_l_sym_ordered = order_l_l_sym(l_l_sym_expend_sym) - #======== #MO COEF #======== @@ -220,14 +222,6 @@ def order_phase(mo_coef): mo_coef_phase.append(ii) return mo_coef_phase - -def order_by_sim(mo_coef, l_l_sym): - l_sym_oder = [int(l[0]) - 1 for l in l_l_sym] - mo_coef_order = [[x for (y, x) in sorted(zip(l_sym_oder, i))] - for i in mo_coef] - return mo_coef_order - - def chunked(l, chunks_size): l_block = [] for i in l: @@ -264,9 +258,8 @@ def print_mo_coef(mo_coef_block, l_l_sym): mo_coef = ezfio.get_mo_basis_mo_coef() -#mo_coef_phase = order_phase(mo_coef) mo_coef_phase = mo_coef -mo_coef_phase_order = order_by_sim(mo_coef_phase, l_l_sym_ordered) +mo_coef_phase_order = mo_coef_phase mo_coef_transp = zip(*mo_coef_phase_order) mo_coef_block = chunked(mo_coef_transp, 4) print_mo_coef(mo_coef_block, l_l_sym_ordered) From 7c91cfbe52356d5f84731f3b84c7ab7ae3328c1c Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Mon, 8 Feb 2016 14:39:23 +0100 Subject: [PATCH 20/29] qp -> qmcpack support now more than 64 MO for the det representation --- .../qmcpack/qp_convert_qmcpack_to_ezfio.py | 27 ++++++++++++------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py index a4a15745..9349b9e2 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py @@ -258,9 +258,7 @@ def print_mo_coef(mo_coef_block, l_l_sym): mo_coef = ezfio.get_mo_basis_mo_coef() -mo_coef_phase = mo_coef -mo_coef_phase_order = mo_coef_phase -mo_coef_transp = zip(*mo_coef_phase_order) +mo_coef_transp = zip(*mo_coef) mo_coef_block = chunked(mo_coef_transp, 4) print_mo_coef(mo_coef_block, l_l_sym_ordered) @@ -336,14 +334,25 @@ psi_coef = ezfio.get_determinants_psi_coef()[0] for c, (l_det_bit_alpha, l_det_bit_beta) in zip(psi_coef, psi_det): print c - for det in l_det_bit_alpha: - bin_det_raw = "{0:b}".format(det)[::-1] - bin_det = bin_det_raw + "0" * (mo_num - len(bin_det_raw)) + + bin_det = "" + for i,int_det in enumerate(l_det_bit_alpha): + bin_det_raw = "{0:b}".format(int_det)[::-1] + if mo_num - 64*(i+1) > 0: + bin_det += bin_det_raw + "0" * (64*(i+1) - len(bin_det_raw)) + else: + bin_det += bin_det_raw + "0" * (mo_num-64*i - len(bin_det_raw)) + print bin_det - for det in l_det_bit_beta: - bin_det_raw = "{0:b}".format(det)[::-1] - bin_det = bin_det_raw + "0" * (mo_num - len(bin_det_raw)) + bin_det = "" + for i,int_det in enumerate(l_det_bit_beta): + bin_det_raw = "{0:b}".format(int_det)[::-1] + if mo_num - 64*(i+1) > 0: + bin_det += bin_det_raw + "0" * (64*(i+1) - len(bin_det_raw)) + else: + bin_det += bin_det_raw + "0" * (mo_num-64*i - len(bin_det_raw)) + print bin_det print "" From 6d30dabc8b6aa7cbd983a6897f743f7e1e5352b2 Mon Sep 17 00:00:00 2001 From: Lorenzo Tenti Date: Fri, 12 Feb 2016 12:27:13 +0100 Subject: [PATCH 21/29] Added the Orbital Entanglement plugin, added the list of active orbitals in bitmasks.irp.f. --- ocaml/qp_edit.ml | 36 +- .../NEEDED_CHILDREN_MODULES | 1 + .../Orbital_Entanglement.irp.f | 345 ++++++++++++++++++ plugins/Orbital_Entanglement/README.rst | 65 ++++ .../print_entanglement.irp.f | 46 +++ plugins/Orbital_Entanglement/tree_dependency | 21 ++ .../Orbital_Entanglement/tree_dependency.png | 0 scripts/entanglement.py | 140 +++++++ src/Bitmask/bitmasks.irp.f | 28 ++ 9 files changed, 664 insertions(+), 18 deletions(-) create mode 100644 plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES create mode 100644 plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f create mode 100644 plugins/Orbital_Entanglement/README.rst create mode 100644 plugins/Orbital_Entanglement/print_entanglement.irp.f create mode 100644 plugins/Orbital_Entanglement/tree_dependency create mode 100644 plugins/Orbital_Entanglement/tree_dependency.png create mode 100755 scripts/entanglement.py diff --git a/ocaml/qp_edit.ml b/ocaml/qp_edit.ml index a693aa2f..53e3ea59 100644 --- a/ocaml/qp_edit.ml +++ b/ocaml/qp_edit.ml @@ -17,12 +17,12 @@ type keyword = | Electrons | Mo_basis | Nuclei -| Determinants -| Integrals_bielec +| Hartree_fock | Pseudo +| Integrals_bielec | Perturbation | Properties -| Hartree_fock +| Determinants ;; @@ -32,12 +32,12 @@ let keyword_to_string = function | Electrons -> "Electrons" | Mo_basis -> "MO basis" | Nuclei -> "Molecule" -| Determinants -> "Determinants" -| Integrals_bielec -> "Integrals_bielec" +| Hartree_fock -> "Hartree_fock" | Pseudo -> "Pseudo" +| Integrals_bielec -> "Integrals_bielec" | Perturbation -> "Perturbation" | Properties -> "Properties" -| Hartree_fock -> "Hartree_fock" +| Determinants -> "Determinants" ;; @@ -86,18 +86,18 @@ let get s = f Ao_basis.(read, to_rst) | Determinants_by_hand -> f Determinants_by_hand.(read, to_rst) - | Determinants -> - f Determinants.(read, to_rst) - | Integrals_bielec -> - f Integrals_bielec.(read, to_rst) + | Hartree_fock -> + f Hartree_fock.(read, to_rst) | Pseudo -> f Pseudo.(read, to_rst) + | Integrals_bielec -> + f Integrals_bielec.(read, to_rst) | Perturbation -> f Perturbation.(read, to_rst) | Properties -> f Properties.(read, to_rst) - | Hartree_fock -> - f Hartree_fock.(read, to_rst) + | Determinants -> + f Determinants.(read, to_rst) end with | Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "") @@ -135,12 +135,12 @@ let set str s = in let open Input in match s with - | Determinants -> write Determinants.(of_rst, write) s - | Integrals_bielec -> write Integrals_bielec.(of_rst, write) s + | Hartree_fock -> write Hartree_fock.(of_rst, write) s | Pseudo -> write Pseudo.(of_rst, write) s + | Integrals_bielec -> write Integrals_bielec.(of_rst, write) s | Perturbation -> write Perturbation.(of_rst, write) s | Properties -> write Properties.(of_rst, write) s - | Hartree_fock -> write Hartree_fock.(of_rst, write) s + | Determinants -> write Determinants.(of_rst, write) s | Electrons -> write Electrons.(of_rst, write) s | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s | Nuclei -> write Nuclei.(of_rst, write) s @@ -188,12 +188,12 @@ let run check_only ezfio_filename = Nuclei ; Ao_basis; Electrons ; - Determinants ; - Integrals_bielec ; + Hartree_fock ; Pseudo ; + Integrals_bielec ; Perturbation ; Properties ; - Hartree_fock ; + Determinants ; Mo_basis; Determinants_by_hand ; ] diff --git a/plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES b/plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..aae89501 --- /dev/null +++ b/plugins/Orbital_Entanglement/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants diff --git a/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f b/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f new file mode 100644 index 00000000..74d7d8ae --- /dev/null +++ b/plugins/Orbital_Entanglement/Orbital_Entanglement.irp.f @@ -0,0 +1,345 @@ +use bitmasks + +BEGIN_PROVIDER [integer, mo_inp_num] + implicit none + BEGIN_DOC +! This is the number of orbitals involved in the entanglement calculation. +! It is taken equal to the number of active orbitals n_act_orb. + END_DOC + mo_inp_num = n_act_orb + +END_PROVIDER + + BEGIN_PROVIDER [integer, mo_inp_list, (N_int*bit_kind_size)] +&BEGIN_PROVIDER [integer, mo_inp_list_rev, (mo_tot_num)] +&BEGIN_PROVIDER [integer(bit_kind), mo_inp_bit_list, (N_int)] + implicit none + BEGIN_DOC +! mo_inp_list is the list of the orbitals involved in the entanglement calculation. +! It is taken equal to the list of active orbitals list_act. +! mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i. + END_DOC + integer :: i + + do i = 1, mo_inp_num + mo_inp_list(i)=list_act(i) + enddo + + do i = 1, mo_inp_num + mo_inp_list_rev(mo_inp_list(i))=i + enddo + call list_to_bitstring( mo_inp_bit_list, mo_inp_list, mo_inp_num, N_int) +END_PROVIDER + + + BEGIN_PROVIDER [double precision, entropy_one_orb, (mo_inp_num)] +&BEGIN_PROVIDER [double precision, entropy_two_orb, (mo_inp_num,mo_inp_num)] + implicit none + BEGIN_DOC +! entropy_one_orb is the one-orbital von Neumann entropy S(1)_i +! entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij. + END_DOC + + double precision, allocatable :: ro1(:,:),ro2(:,:,:) + integer :: i,j,k,l,ii,jj,iii,istate,kl,info + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha, n_occ_beta + logical, allocatable :: zocc(:,:) + logical :: zalpha, zbeta, zalpha2, zbeta2 + integer :: exc(0:2,2,2),degree,h1,p1,h2,p2,spin1,spin2 + double precision :: phase + integer(bit_kind) :: key_tmp(N_int), key_tmp2(N_int) + integer :: ip + double precision, parameter :: eps=10.d0**(-14) + double precision :: w(16), work(3*16) + + + allocate(ro1(4,mo_inp_num),ro2(16,16,(mo_inp_num*(mo_inp_num-1)/2))) + + entropy_one_orb = 0.d0 + entropy_two_orb = 0.d0 + ro1 = 0.d0 + ro2 = 0.d0 + + allocate (occ(N_int*bit_kind_size,2)) + allocate (zocc(mo_tot_num,2)) + + istate = 1 !Only GS, to be generalized... + do ii=1,N_det + ! We get the occupation of the alpha electrons in occ(:,1) + call bitstring_to_list(psi_det(1,1,ii), occ(1,1), n_occ_alpha, N_int) + ! We get the occupation of the beta electrons in occ(:,2) + call bitstring_to_list(psi_det(1,2,ii), occ(1,2), n_occ_beta, N_int) + zocc = .false. + do i=1,n_occ_alpha + zocc(occ(i,1),1)=.true. + enddo + do i=1,n_occ_beta + zocc(occ(i,2),2)=.true. + enddo + + do k=1,mo_inp_num + zalpha = zocc(mo_inp_list(k),1) + zbeta = zocc(mo_inp_list(k),2) +! mono start + if (zbeta.and.zalpha) then + ro1(4,k) = ro1(4,k) + psi_coef(ii,istate)**2 ! double occupied + elseif (zalpha) then + ro1(2,k) = ro1(2,k) + psi_coef(ii,istate)**2 ! single alpha + elseif (zbeta) then + ro1(3,k) = ro1(3,k) + psi_coef(ii,istate)**2 ! single beta + else + ro1(1,k) = ro1(1,k) + psi_coef(ii,istate)**2 ! empty + endif +! mono stop +! double start + if (k.eq.mo_inp_num) cycle + do l=k+1,mo_inp_num + kl=(l-1)*(l-2)/2+k + zalpha2 = zocc(mo_inp_list(l),1) + zbeta2 = zocc(mo_inp_list(l),2) + + if (zbeta.and.zalpha.and.zbeta2.and.zalpha2) then + ro2(16,16,kl) = ro2(16,16,kl) + psi_coef(ii,istate)**2 ! both double occupied + else if (zbeta.and.zalpha.and.zbeta2) then + ro2(15,15,kl) = ro2(15,15,kl) + psi_coef(ii,istate)**2 ! one double, one beta + else if (zbeta.and.zalpha.and.zalpha2) then + ro2(13,13,kl) = ro2(13,13,kl) + psi_coef(ii,istate)**2 ! one double, one alpha + else if (zbeta.and.zbeta2.and.zalpha2) then + ro2(14,14,kl) = ro2(14,14,kl) + psi_coef(ii,istate)**2 ! one beta, one double + else if (zalpha.and.zbeta2.and.zalpha2) then + ro2(12,12,kl) = ro2(12,12,kl) + psi_coef(ii,istate)**2 ! one alpha, one double + else if (zalpha.and.zbeta) then + ro2(11,11,kl) = ro2(11,11,kl) + psi_coef(ii,istate)**2 ! one double, one empty + else if (zbeta2.and.zalpha2) then + ro2(8,8,kl) = ro2(8,8,kl) + psi_coef(ii,istate)**2 ! one empty, one double + else if (zbeta.and.zalpha2) then + ro2(10,10,kl) = ro2(10,10,kl) + psi_coef(ii,istate)**2 ! one beta, one alpha + else if (zalpha.and.zbeta2) then + ro2(9,9,kl) = ro2(9,9,kl) + psi_coef(ii,istate)**2 ! one alpha, one beta + else if (zbeta.and.zbeta2) then + ro2(7,7,kl) = ro2(7,7,kl) + psi_coef(ii,istate)**2 ! one beta, one beta + else if (zalpha.and.zalpha2) then + ro2(6,6,kl) = ro2(6,6,kl) + psi_coef(ii,istate)**2 ! one alpha, one alpha + else if (zbeta) then + ro2(5,5,kl) = ro2(5,5,kl) + psi_coef(ii,istate)**2 ! one beta, one empty + else if (zbeta2) then + ro2(4,4,kl) = ro2(4,4,kl) + psi_coef(ii,istate)**2 ! one empty, one beta + else if (zalpha) then + ro2(3,3,kl) = ro2(3,3,kl) + psi_coef(ii,istate)**2 ! one alpha, one empty + else if (zalpha2) then + ro2(2,2,kl) = ro2(2,2,kl) + psi_coef(ii,istate)**2 ! one empty, one alpha + else + ro2(1,1,kl) = ro2(1,1,kl) + psi_coef(ii,istate)**2 ! both empty + end if + enddo + enddo +! stop double + + if (ii.eq.N_det) cycle +!Off Diagonal Elements + do jj=ii+1,N_det + + call get_excitation_degree(psi_det(1,1,ii),psi_det(1,1,jj),degree,N_int) + if (degree.gt.2) cycle + ip=0 + do iii =1,N_int + key_tmp(iii) = ior(xor(psi_det(iii,1,ii),psi_det(iii,1,jj)),xor(psi_det(iii,2,ii),psi_det(iii,2,jj))) + ip += popcnt(key_tmp(iii)) + enddo + if (ip.ne.2) cycle !They involve more than 2 orbitals. + ip=0 + do iii=1,N_int + ip += popcnt(iand(key_tmp(iii),mo_inp_bit_list(iii))) + enddo + if (ip.ne.2) cycle !They do not involve orbitals of the list. + + if (degree.eq.2) then + call get_double_excitation(psi_det(1,1,ii),psi_det(1,1,jj),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,spin1,spin2) + k=mo_inp_list_rev(h1) + l=mo_inp_list_rev(p1) + if (k.gt.l) then + kl=(k-1)*(k-2)/2+l + else + kl=(l-1)*(l-2)/2+k + endif + + if ((.not.zocc(mo_inp_list(l),1)).and.(.not.zocc(mo_inp_list(l),2))& + .and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then + ro2(8,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(11,8,kl) = ro2(8,11,kl) + endif + + if ((zocc(mo_inp_list(l),1)).and.(.not.zocc(mo_inp_list(l),2))& + .and.(.not.zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then + ro2(9,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative + ro2(10,9,kl) = ro2(9,10,kl) + endif + if ((zocc(mo_inp_list(k),1)).and.(.not.zocc(mo_inp_list(k),2))& + .and.(.not.zocc(mo_inp_list(l),1)).and.(zocc(mo_inp_list(l),2))) then + ro2(9,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative + ro2(10,9,kl) = ro2(9,10,kl) + endif + endif + + if (degree.eq.1) then + call get_mono_excitation(psi_det(1,1,ii),psi_det(1,1,jj),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,spin1,spin2) + k=mo_inp_list_rev(h1) + l=mo_inp_list_rev(p1) + if (k.gt.l) then + kl=(k-1)*(k-2)/2+l + else + kl=(l-1)*(l-2)/2+k + endif + + if ((.not.(zocc(mo_inp_list(l),2))).and.& + (.not.(zocc(mo_inp_list(l),1))).and.(zocc(mo_inp_list(k),1))& + .and.(.not.(zocc(mo_inp_list(k),2)))) then + ro2(2,3,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(3,2,kl)=ro2(2,3,kl) + endif + + if ((.not.(zocc(mo_inp_list(l),2))).and.& + (.not.(zocc(mo_inp_list(l),1))).and.(zocc(mo_inp_list(k),2))& + .and.(.not.(zocc(mo_inp_list(k),1)))) then + ro2(4,5,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(5,4,kl)=ro2(4,5,kl) + endif + + + if ((.not.(zocc(mo_inp_list(l),1))).and.& !k doubly occupied, l empty + (.not.(zocc(mo_inp_list(l),2))).and.& + (zocc(mo_inp_list(k),1)).and.& + (zocc(mo_inp_list(k),2))) then + if (k.gt.l) then + if (spin1.eq.1) then !spin alpha + ro2(8,9,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(9,8,kl)=ro2(8,9,kl) + else !spin beta + ro2(8,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative + ro2(10,8,kl)=ro2(8,10,kl) + endif + else ! k.lt.l + if (spin1.eq.1) then !spin alpha + ro2(10,11,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative + ro2(11,10,kl)=ro2(10,11,kl) + else !spin beta + ro2(9,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(11,9,kl)=ro2(9,11,kl) + endif + endif + endif + + if ((.not.(zocc(mo_inp_list(l),1))).and.& !k alpha, l beta + (.not.(zocc(mo_inp_list(k),2))).and.& + (zocc(mo_inp_list(k),1)).and.& + (zocc(mo_inp_list(l),2))) then + if (k.gt.l) then + if (spin1.eq.1) then !spin alpha + ro2(10,11,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative + ro2(11,10,kl)=ro2(10,11,kl) + else !spin beta + print*, "problem in k alpha l beta k.gt.l spin beta" + endif + else ! k.lt.l + if (spin1.eq.1) then !spin alpha + ro2(8,9,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(9,8,kl)=ro2(8,9,kl) + else !spin beta + print*, "problem in k alpha l beta k.lt.l spin beta" + endif + endif + endif + + if ((.not.(zocc(mo_inp_list(k),1))).and.& !k beta, l alpha + (.not.(zocc(mo_inp_list(l),2))).and.& + (zocc(mo_inp_list(l),1)).and.& + (zocc(mo_inp_list(k),2))) then + if (k.gt.l) then + if (spin1.eq.2) then !spin beta + ro2(9,11,kl) += phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(11,9,kl)=ro2(9,11,kl) + else !spin alpha + print*, "problem in k beta l alpha k.gt.l spin alpha" + endif + else ! k.lt.l + if (spin1.eq.2) then !spin beta + ro2(8,10,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) !negative + ro2(10,8,kl)=ro2(8,10,kl) + else !spin alpha + print*, "problem in k beta l alpha k.lt.l spin alpha" + endif + endif + endif + + if (zocc(mo_inp_list(l),1).and.(.not.zocc(mo_inp_list(l),2))& + .and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then + ro2(12,13,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(13,12,kl) = ro2(12,13,kl) + endif + + if (zocc(mo_inp_list(l),2).and.(.not.zocc(mo_inp_list(l),1))& + .and.(zocc(mo_inp_list(k),1)).and.(zocc(mo_inp_list(k),2))) then + ro2(14,15,kl) -= phase*psi_coef(ii,istate)*psi_coef(jj,istate) + ro2(15,14,kl) = ro2(14,15,kl) + endif + endif + enddo + enddo + + + + entropy_one_orb=0.d0 + do k=1,mo_inp_num + do i=1,4 + if (ro1(i,k).ge.eps) then + entropy_one_orb(k) = entropy_one_orb(k)-ro1(i,k)*log(ro1(i,k)) + endif + enddo + enddo + + entropy_two_orb=0.d0 + do k=1,mo_inp_num + do l=1,k + if (k.eq.l) cycle + kl=(k-1)*(k-2)/2+l + call dsyev('N','U',16,ro2(1,1,kl),16,w,work,3*16,info) + if (info.ne.0) then + write(*,*) "Errore in dsyev" + endif + do j=1,16 + if (w(j).ge.eps) then + entropy_two_orb(k,l) = entropy_two_orb(k,l)-w(j)*log(w(j)) + entropy_two_orb(l,k) = entropy_two_orb(k,l) + elseif ((w(j)).lt.(-eps)) then + write(6,*) "Negative Eigenvalue. You have a big problem..." + write(6,*) w(j) + stop + endif + enddo + enddo + enddo + + + deallocate (occ,zocc,ro1,ro2) + +END_PROVIDER + +BEGIN_PROVIDER [double precision, mutinf, (mo_inp_num,mo_inp_num)] +implicit none +BEGIN_DOC +!mutinf is the mutual information (entanglement), calculated as I_ij=0.5*[S(1)_i+S(1)_j-S(2)_ij] +!see the refence: 10.1016/j.chemphys.2005.10.018 +END_DOC + integer :: i,j +! mutal information: + mutinf = 0.d0 + do i=1,mo_inp_num + do j=1,mo_inp_num + if (j.eq.i) cycle + mutinf(i,j)=-0.5d0*(entropy_two_orb(i,j)-entropy_one_orb(i)-entropy_one_orb(j)) + enddo + enddo +END_PROVIDER diff --git a/plugins/Orbital_Entanglement/README.rst b/plugins/Orbital_Entanglement/README.rst new file mode 100644 index 00000000..9b5b8119 --- /dev/null +++ b/plugins/Orbital_Entanglement/README.rst @@ -0,0 +1,65 @@ +==================== +Orbital_Entanglement +==================== + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +.. image:: tree_dependency.png + +* `Determinants `_ + +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +`entropy_one_orb `_ + entropy_one_orb is the one-orbital von Neumann entropy S(1)_i + entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij. + + +`entropy_two_orb `_ + entropy_one_orb is the one-orbital von Neumann entropy S(1)_i + entropy_two_orb is the two-orbital von Neumann entropy S(2)_ij. + + +`mo_inp_bit_list `_ + mo_inp_list is the list of the orbitals involved in the entanglement calculation. + It is taken equal to the list of active orbitals list_act. + mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i. + + +`mo_inp_list `_ + mo_inp_list is the list of the orbitals involved in the entanglement calculation. + It is taken equal to the list of active orbitals list_act. + mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i. + + +`mo_inp_list_rev `_ + mo_inp_list is the list of the orbitals involved in the entanglement calculation. + It is taken equal to the list of active orbitals list_act. + mo_inp_list_rev is a list such that mo_inp_list_rev(mo_inp_list(i))=i. + + +`mo_inp_num `_ + This is the number of orbitals involved in the entanglement calculation. + It is taken equal to the number of active orbitals n_act_orb. + + +`mutinf `_ + mutinf is the mutual information (entanglement), calculated as I_ij=0.5*[S(1)_i+S(1)_j-S(2)_ij] + see the refence: 10.1016/j.chemphys.2005.10.018 + + +`pouet `_ + Undocumented + + +`routine `_ + Undocumented + diff --git a/plugins/Orbital_Entanglement/print_entanglement.irp.f b/plugins/Orbital_Entanglement/print_entanglement.irp.f new file mode 100644 index 00000000..e7e732cc --- /dev/null +++ b/plugins/Orbital_Entanglement/print_entanglement.irp.f @@ -0,0 +1,46 @@ +program pouet + + implicit none + read_wf = .true. + touch read_wf + call routine +end + +subroutine routine + implicit none + integer:: i,j + + write(6,*) 'Total orbitals: ', mo_tot_num + write(6,*) 'Orbitals for entanglement calculation: ', mo_inp_num + write(6,*) 'Index: ',(mo_inp_list(i),i=1,mo_inp_num) + write(6,*) + write(6,*) "s1: One-Orbital von Neumann entropy" + write(6,'(1000f8.5)') entropy_one_orb + write(6,*) + write(6,*) "s2: Two-Orbital von Neumann entropy" + do i=1,mo_inp_num + write(6,'(1000f8.5)') (entropy_two_orb(i,j),j=1,mo_inp_num) + enddo + write(6,*) + +! mutal information: + write(6,*) "Mutual Information (Entanglement)" + do i=1,mo_inp_num + write(6,'(1000f8.5)') (mutinf(i,j),j=1,mo_inp_num) + enddo + + + open(17,file=(trim(ezfio_filename)//".entanglement"),status='unknown',form='formatted') + + write(17,'(1000f8.5)') entropy_one_orb + do i=1,mo_inp_num + write(17,'(1000f8.5)') (mutinf(i,j),j=1,mo_inp_num) + enddo + + + close(17) + write(6,*) + write(6,*) "The .entanglement file is the input for the entanglement.py script." + write(6,*) "You can find the script in the directory Scripts of QP." + write(6,*) +end diff --git a/plugins/Orbital_Entanglement/tree_dependency b/plugins/Orbital_Entanglement/tree_dependency new file mode 100644 index 00000000..a601745c --- /dev/null +++ b/plugins/Orbital_Entanglement/tree_dependency @@ -0,0 +1,21 @@ +// ['Orbital_Entanglement'] +digraph { + Orbital_Entanglement [fontcolor=red] + Orbital_Entanglement -> Determinants + Determinants -> Integrals_Monoelec + Integrals_Monoelec -> MO_Basis + MO_Basis -> AO_Basis + AO_Basis -> Nuclei + Nuclei -> Ezfio_files + Nuclei -> Utils + MO_Basis -> Electrons + Electrons -> Ezfio_files + Integrals_Monoelec -> Pseudo + Pseudo -> Nuclei + Determinants -> Integrals_Bielec + Integrals_Bielec -> Pseudo + Integrals_Bielec -> Bitmask + Bitmask -> MO_Basis + Integrals_Bielec -> ZMQ + ZMQ -> Utils +} \ No newline at end of file diff --git a/plugins/Orbital_Entanglement/tree_dependency.png b/plugins/Orbital_Entanglement/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/scripts/entanglement.py b/scripts/entanglement.py new file mode 100755 index 00000000..c55cc98e --- /dev/null +++ b/scripts/entanglement.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +import sys +import matplotlib.pyplot as plt +from matplotlib import lines +from pylab import * + +if __name__ == '__main__': + inputfile = sys.argv[1] +with open(str(inputfile),'r') as f: + vec_s1 = [float(x) for x in f.readline().split()] + mat_I = np.array([[float(x) for x in ln.split()] for ln in f]) +#mat_I = np.loadtxt(str(inputfile),skiprows=1) + +print vec_s1 +print mat_I +#enable the following lines if you dont have a working X display +#import matplotlib +#matplotlib.use('Agg') + +#MUTUALI INFORMATION PLOT: + +# edit the following line to select the orbitals... +zselect=False +if zselect: + select =[] + select.extend(range(34,132)) + select = [x-1 for x in select] + nselect = len(select) + vec_s1_s = np.zeros(nselect) + mat_I_s = np.zeros((nselect,nselect)) + print mat_I_s, vec_s1_s + + for i,x in enumerate(select): + vec_s1_s[i]=vec_s1[x] + for i,x in enumerate(select): + for j,y in enumerate(select): + mat_I_s[i,j] = mat_I[x,y] + + + print vec_s1_s + print mat_I_s + vec_s1 = vec_s1_s + mat_I = mat_I_s + + + + +#end selection + + +N= len(mat_I) +print N +theta = np.zeros(N) +r=np.zeros(N) +labels=np.zeros(N) +area=np.zeros(N) + +for i in range(N): + theta[i]=-2*pi/N*i+pi/2 + r[i]=1.0 + if zselect: + labels[i]=select[i]+1 + else: + labels[i]=i+1 + area[i]=vec_s1[i]*500 + +if (len(sys.argv) > 2): + for j in range(2,len(sys.argv)): + if (sys.argv[j] == '-c'): + coordfile = sys.argv[j+1] + with open(str(coordfile),'r') as f: + for i in range(N): + line = f.readline().split() + theta[i] = float(line[0]) +pi/2 + r[i] = float(line[1]) +print 'theta=',theta +print 'r=',r + + + + +ax = plt.subplot(111,polar=True) +ax.set_xticklabels([]) +ax.set_yticklabels([]) +ax.grid(b=False) +c = plt.scatter(theta,r,c="Red",s=area) + + +if (len(sys.argv) > 2): + for j in range(2,len(sys.argv)): + print 'sys.argv',j,sys.argv[j] + if (sys.argv[j] == '-t'): + plt.title(sys.argv[j+1]) + break + else: + plt.title('N = '+str(N/2)+', Results file: '+sys.argv[1]) + +#this is dummy: +#c1 = plt.scatter(theta,(r+0.05),c="red",s=0) + +for i in range(N): +# plt.annotate(int(labels[i]),xy=(theta[i],(r[i]+0.2)),size='xx-large',) + plt.text(theta[i],(r[i]+0.25),int(labels[i]),size='x-small',ha='center',va='center') + for j in range(i,N): + x =[theta[i],theta[j]] + y =[r[i],r[j]] + #y =[1,1] + if mat_I[i,j] >= 0.1: + line1 = lines.Line2D(x, y, linewidth=3, color='black',linestyle='-', alpha=0.8,label='0.1') + ax.add_line(line1) + elif mat_I[i,j] >=0.01: + line2 = lines.Line2D(x, y, linewidth=3, color='green',linestyle='-', alpha=0.8,label='0.01') + ax.add_line(line2) + elif mat_I[i,j] >=0.001: + line3 = lines.Line2D(x, y, linewidth=3, color='gray',linestyle='-', alpha=0.6,label='0.001') + ax.add_line(line3) +params = {'legend.fontsize': 16, + 'legend.linewidth': 2} + +#need to check if the three types of lines really exist. +if 'line1' in locals() and 'line2' in locals() and 'line3' in locals(): + ax.legend([line1, line2, line3],[line1.get_label(),line2.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True) +elif 'line1' in locals() and 'line2' in locals(): + ax.legend([line1, line2],[line1.get_label(),line2.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True) +elif 'line1' in locals() and 'line3' in locals(): + ax.legend([line1, line3],[line1.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True) +elif 'line2' in locals() and 'line3' in locals(): + ax.legend([line2, line3],[line2.get_label(),line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True) +elif 'line1' in locals(): + ax.legend([line1],[line1.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True) +elif 'line2' in locals(): + ax.legend([line3],[line2.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True) +elif 'line3' in locals(): + ax.legend([line2],[line3.get_label()],bbox_to_anchor=(0.00,1.0),fancybox=True,shadow=True) +#probably not the best way to code it. + +#for saving, enable matplotlib.use('Agg') +#plt.savefig(str(sys.argv[1])+'.eps') +plt.show() diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index f54532ca..70a84a42 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -407,3 +407,31 @@ END_PROVIDER unpaired_alpha_electrons(i) = xor(HF_bitmask(i,1),HF_bitmask(i,2)) enddo END_PROVIDER + +BEGIN_PROVIDER [ integer, n_act_orb] + BEGIN_DOC + ! number of active orbitals + END_DOC + implicit none + integer :: i,j + n_act_orb = 0 + do i = 1, N_int + n_act_orb += popcnt(cas_bitmask(i,1,1)) + enddo +END_PROVIDER + +BEGIN_PROVIDER [integer, list_act, (n_act_orb)] + BEGIN_DOC + ! list of active orbitals + END_DOC + implicit none + integer :: occ_act(N_int*bit_kind_size) + integer :: itest,i + occ_act = 0 + call bitstring_to_list(cas_bitmask(1,1,1), occ_act(1), itest, N_int) + ASSERT(itest==n_act_orb) + do i = 1, n_act_orb + list_act(i) = occ_act(i) + enddo + +END_PROVIDER From acb8c1a22596d26af0fe2da89ffd36a49561f242 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 12 Feb 2016 13:48:43 +0100 Subject: [PATCH 22/29] Update README.md --- README.md | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/README.md b/README.md index ce90ee61..e313f444 100644 --- a/README.md +++ b/README.md @@ -148,3 +148,20 @@ You have two or more ezfio configuration files for the same variable. Check file - rm $QP_ROOT/install/EZFIO/config/* - ninja + + +### Error: Seg Fault (139) + +``` +Segmentation fault (core dumped) +Program exited with code 139. +``` + +#### Why ? + +It's caused when we call the DGEM routine of LAPACK. + +##### Fix + +Set `ulimit -s unlimited`, before runing `qp_run`. It seem to fix the problem. + From 5bad973ba625d4176e82d9843475da3969b84e44 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 15 Feb 2016 08:15:43 +0100 Subject: [PATCH 23/29] Don't remove python executable when make clean --- scripts/module/module_handler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/module/module_handler.py b/scripts/module/module_handler.py index 136dc8cf..0667c376 100755 --- a/scripts/module/module_handler.py +++ b/scripts/module/module_handler.py @@ -40,7 +40,7 @@ def is_plugin(path_module_rel): def is_exe(fpath): - return os.path.isfile(fpath) and os.access(fpath, os.X_OK) + return os.path.isfile(fpath) and os.access(fpath, os.X_OK) and not fpath.endswith(".py") def get_dict_child(l_root_abs=None): From 92dd252be00138d1811fbc22493c4bf3c89d3d46 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 16 Feb 2016 11:14:19 +0100 Subject: [PATCH 24/29] added the All_singles module --- plugins/All_singles/NEEDED_CHILDREN_MODULES | 1 + plugins/All_singles/README.rst | 12 ++++ plugins/All_singles/all_1h_1p.irp.f | 75 ++++++++++++++++++++ plugins/All_singles/all_singles.irp.f | 76 +++++++++++++++++++++ scripts/generate_h_apply.py | 14 ++++ src/Bitmask/bitmask_cas_routines.irp.f | 44 ++++++++++++ src/Determinants/H_apply.template.f | 5 ++ 7 files changed, 227 insertions(+) create mode 100644 plugins/All_singles/NEEDED_CHILDREN_MODULES create mode 100644 plugins/All_singles/README.rst create mode 100644 plugins/All_singles/all_1h_1p.irp.f create mode 100644 plugins/All_singles/all_singles.irp.f diff --git a/plugins/All_singles/NEEDED_CHILDREN_MODULES b/plugins/All_singles/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..bb97ddb9 --- /dev/null +++ b/plugins/All_singles/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Generators_restart Perturbation Properties Selectors_no_sorted Utils diff --git a/plugins/All_singles/README.rst b/plugins/All_singles/README.rst new file mode 100644 index 00000000..b4b3f517 --- /dev/null +++ b/plugins/All_singles/README.rst @@ -0,0 +1,12 @@ +=========== +All_singles +=========== + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/All_singles/all_1h_1p.irp.f b/plugins/All_singles/all_1h_1p.irp.f new file mode 100644 index 00000000..e6d48749 --- /dev/null +++ b/plugins/All_singles/all_1h_1p.irp.f @@ -0,0 +1,75 @@ +program restart_more_singles + BEGIN_DOC + ! Generates and select single and double excitations of type 1h-1p + ! on the top of a given restart wave function of type CAS + END_DOC + read_wf = .true. + touch read_wf + print*,'ref_bitmask_energy = ',ref_bitmask_energy + call routine + +end +subroutine routine + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:) + integer :: N_st, degree + integer :: n_det_before + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + i = 0 + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + pt2=-1.d0 + E_before = ref_bitmask_energy + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_1h_1p(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + E_before = CI_energy + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + call save_wavefunction + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + + enddo + + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1)) + enddo + endif + + call save_wavefunction + deallocate(pt2,norm_pert) +end diff --git a/plugins/All_singles/all_singles.irp.f b/plugins/All_singles/all_singles.irp.f new file mode 100644 index 00000000..3b5c5cce --- /dev/null +++ b/plugins/All_singles/all_singles.irp.f @@ -0,0 +1,76 @@ +program restart_more_singles + BEGIN_DOC + ! Generates and select single excitations + ! on the top of a given restart wave function + END_DOC + read_wf = .true. + touch read_wf + call routine + +end +subroutine routine + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + double precision,allocatable :: E_before(:) + integer :: n_det_before + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + i = 0 + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + pt2=-1.d0 + E_before = ref_bitmask_energy + pt2_max = 1.d-10 + n_det_max = 200000 + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_mono(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy + print*,'pt2 = ',pt2 + print*,'E+PT2 = ',E_before + pt2 + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + E_before = CI_energy + call save_wavefunction + + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + enddo + + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1)) + enddo + endif + call save_wavefunction + deallocate(pt2,norm_pert,E_before) +end diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 02524c3d..b52e4445 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -24,6 +24,8 @@ skip init_main filter_integrals filter2h2p +filter_only_1h1p_single +filter_only_1h1p_double filterhole filterparticle do_double_excitations @@ -150,6 +152,18 @@ class H_apply(object): self["filterparticle"] = """ if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind )cycle """ + + def filter_only_1h1p(self): + self["filter_only_1h1p_single"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h1p(hole).eq..False.) cycle + """ + self["filter_only_1h1p_double"] = """ +! ! DIR$ FORCEINLINE + if (is_a_1h1p(key).eq..False.) cycle + """ + + def unset_skip(self): self["skip"] = """ """ diff --git a/src/Bitmask/bitmask_cas_routines.irp.f b/src/Bitmask/bitmask_cas_routines.irp.f index 776d4546..4441fb22 100644 --- a/src/Bitmask/bitmask_cas_routines.irp.f +++ b/src/Bitmask/bitmask_cas_routines.irp.f @@ -445,3 +445,47 @@ integer function number_of_particles_verbose(key_in) + popcnt( iand( iand( xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1))), virt_bitmask(1,2) ), virt_bitmask(1,2)) ) end +logical function is_a_1h1p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1h1p = .False. + if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.1)then + is_a_1h1p = .True. + endif + +end + +logical function is_a_1h(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1h = .False. + if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.0)then + is_a_1h = .True. + endif + +end + +logical function is_a_1p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_1p = .False. + if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.1)then + is_a_1p = .True. + endif + +end + +logical function is_a_2p(key_in) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: number_of_particles, number_of_holes + is_a_2p = .False. + if(number_of_holes(key_in).eq.0 .and. number_of_particles(key_in).eq.2)then + is_a_2p = .True. + endif + +end + diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 48e5d335..3f99f705 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -165,6 +165,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl endif logical :: check_double_excitation + logical :: is_a_1h1p logical :: b_cycle check_double_excitation = .True. iproc = iproc_in @@ -298,6 +299,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,other_spin) = ibset(key(k,other_spin),l) $filter2h2p + $filter_only_1h1p_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -348,6 +350,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,ispin) = ibset(key(k,ispin),l) $filter2h2p + $filter_only_1h1p_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -418,6 +421,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato integer(bit_kind) :: key_mask(N_int, 2) logical :: check_double_excitation + logical :: is_a_1h1p key_mask(:,:) = 0_bit_kind @@ -490,6 +494,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato $filterparticle hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) $filter2h2p + $filter_only_1h1p_single key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = hole(k,1) From 843a9b3b889971a9954dcbcb3e7387e74fefbbc8 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 16 Feb 2016 11:26:24 +0100 Subject: [PATCH 25/29] Added H_apply.irp.f in All_singles --- plugins/All_singles/H_apply.irp.f | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 plugins/All_singles/H_apply.irp.f diff --git a/plugins/All_singles/H_apply.irp.f b/plugins/All_singles/H_apply.irp.f new file mode 100644 index 00000000..d0a41f90 --- /dev/null +++ b/plugins/All_singles/H_apply.irp.f @@ -0,0 +1,18 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("just_1h_1p") +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +s.filter_only_1h1p() +print s + +s = H_apply("just_mono") +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +s.unset_double_excitations() +print s + +END_SHELL + From 21275034057d8a964595010ffddd7ad3de3453a2 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 16 Feb 2016 12:34:04 +0100 Subject: [PATCH 26/29] added tests for all_1h_1p --- data/basis/chipman-dzp | 168 ++++++++++++++++++++++++++++ plugins/All_singles/EZFIO.cfg | 5 + plugins/All_singles/all_1h_1p.irp.f | 1 + tests/bats/qp.bats | 28 +++++ tests/input/dhno.xyz | 7 ++ 5 files changed, 209 insertions(+) create mode 100644 data/basis/chipman-dzp create mode 100644 plugins/All_singles/EZFIO.cfg create mode 100644 tests/input/dhno.xyz diff --git a/data/basis/chipman-dzp b/data/basis/chipman-dzp new file mode 100644 index 00000000..b0178ef4 --- /dev/null +++ b/data/basis/chipman-dzp @@ -0,0 +1,168 @@ +HYDROGEN +S 3 + 1 127.9500000 0.0107360 + 2 19.2406000 0.1195440 + 3 2.8992000 0.9264160 +S 1 + 1 0.6534000 1.0000000 +S 1 + 1 0.1776000 1.0000000 +S 1 + 1 0.0483000 1.0000000 +P 1 + 1 1.0000000 1.0000000 + +BORON +S 5 + 1 2788.4100000 0.0021220 + 2 419.0390000 0.0161710 + 3 96.4683000 0.0783560 + 4 28.0694000 0.2632500 + 5 9.3760000 0.5967290 +S 1 + 1 1.3057000 1.0000000 +S 1 + 1 3.4062000 1.0000000 +S 1 + 1 0.3245000 1.0000000 +S 1 + 1 0.1022000 1.0000000 +S 1 + 1 0.0330000 1.0000000 +P 4 + 1 11.3413000 0.0179870 + 2 2.4360000 0.1103390 + 3 0.6836000 0.3831110 + 4 0.2134000 0.6478600 +P 1 + 1 0.0701000 1.0000000 +P 1 + 1 0.0226000 1.0000000 +D 1 + 1 0.1600000 1.0000000 +D 1 + 1 0.6400000 1.0000000 + +CARBON +S 5 + 1 5909.4400000 0.0020040 + 2 887.4510000 0.0153100 + 3 204.7490000 0.0742930 + 4 59.8376000 0.2533640 + 5 19.9981000 0.6005760 +S 1 + 1 2.6860000 1.0000000 +S 1 + 1 7.1927000 1.0000000 +S 1 + 1 0.7000000 1.0000000 +S 1 + 1 0.2133000 1.0000000 +S 1 + 1 0.0667000 1.0000000 +P 4 + 1 26.7860000 0.0182570 + 2 5.9564000 0.1164070 + 3 1.7074000 0.3901110 + 4 0.5314000 0.6372210 +P 1 + 1 0.1654000 1.0000000 +P 1 + 1 0.0517000 1.0000000 +D 1 + 1 0.3700000 1.0000000 +D 1 + 1 1.4800000 1.0000000 + +NITROGEN +S 5 + 1 5909.4400000 0.0020040 + 2 887.4510000 0.0153100 + 3 204.7490000 0.0742930 + 4 59.8376000 0.2533640 + 5 19.9981000 0.6005760 +S 1 + 1 2.6860000 1.0000000 +S 1 + 1 7.1927000 1.0000000 +S 1 + 1 0.7000000 1.0000000 +S 1 + 1 0.2133000 1.0000000 +S 1 + 1 0.0667000 1.0000000 +P 4 + 1 26.7860000 0.0182570 + 2 5.9564000 0.1164070 + 3 1.7074000 0.3901110 + 4 0.5314000 0.6372210 +P 1 + 1 0.1654000 1.0000000 +P 1 + 1 0.0517000 1.0000000 +D 1 + 1 0.3700000 1.0000000 +D 1 + 1 1.4800000 1.0000000 + +OXYGEN + S 5 + 1 7816.5400000 0.0020310 + 2 1175.8200000 0.0154360 + 3 273.1880000 0.0737710 + 4 81.1696000 0.2476060 + 5 27.1836000 0.6118320 +S 1 + 1 3.4136000 1.0000000 +S 1 + 1 9.5322000 1.0000000 +S 1 + 1 0.9398000 1.0000000 +S 1 + 1 0.2846000 1.0000000 +S 1 + 1 0.0862000 1.0000000 +P 4 + 1 35.1832000 0.0195800 + 2 7.9040000 0.1241890 + 3 2.3051000 0.3947270 + 4 0.7171000 0.6273750 +P 1 + 1 0.2137000 1.0000000 +P 1 + 1 0.0648000 1.0000000 +D 1 + 1 0.5500000 1.0000000 +D 1 + 1 2.2000000 1.0000000 + +FLUORINE +S 5 + 1 9994.7900000 0.0020170 + 2 1506.0300000 0.0152950 + 3 350.2690000 0.0731100 + 4 104.0530000 0.2464200 + 5 34.8432000 0.6125930 +S 1 + 1 4.3688000 1.0000000 +S 1 + 1 12.2164000 1.0000000 +S 1 + 1 1.2078000 1.0000000 +S 1 + 1 0.3634000 1.0000000 +S 1 + 1 0.1101000 1.0000000 +P 4 + 1 44.3555000 0.0208680 + 2 10.0820000 0.1300920 + 3 2.9959000 0.3962190 + 4 0.9383000 0.6203680 +P 1 + 1 0.2733000 1.0000000 +P 1 + 1 0.0828000 1.0000000 +D 1 + 1 0.6650000 1.0000000 +D 1 + 1 2.6600000 1.0000000 diff --git a/plugins/All_singles/EZFIO.cfg b/plugins/All_singles/EZFIO.cfg new file mode 100644 index 00000000..b2498c43 --- /dev/null +++ b/plugins/All_singles/EZFIO.cfg @@ -0,0 +1,5 @@ +[energy] +type: double precision +doc: Calculated Selected all_singles or all_1h_1p energy +interface: ezfio + diff --git a/plugins/All_singles/all_1h_1p.irp.f b/plugins/All_singles/all_1h_1p.irp.f index e6d48749..a2786248 100644 --- a/plugins/All_singles/all_1h_1p.irp.f +++ b/plugins/All_singles/all_1h_1p.irp.f @@ -69,6 +69,7 @@ subroutine routine print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1)) enddo endif + call ezfio_set_all_singles_energy(CI_energy) call save_wavefunction deallocate(pt2,norm_pert) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index ae385fe8..892d1a1b 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -79,11 +79,39 @@ function run_FCI() { eq $energy_pt2 $4 $thresh } +function run_all_1h_1p() { + thresh=1.e-6 + test_exe all_1h_1p || skip + ezfio set_file $1 + ezfio set determinants n_det_max $2 + ezfio set perturbation pt2_max $3 + ezfio set determinants threshold_davidson 1.e-10 + + qp_run all_1h_1p $1 | tee $1.F1h1p.out + energy="$(ezfio get all_singles energy)" + eq $energy $4 $thresh +} + # ___ # | _ _ _|_ # | (/_ _> |_ # + +#=== DHNO +@test "init DHNO chipman-dzp" { + run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio +} + +@test "SCF DHNO chipman-dzp" { + run_HF dhno.ezfio -130.4278777822 +} + +@test "all_1h_1p DHNO chipman-dzp" { + qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio + run_all_1h_1p dhno.ezfio 10000 0.0000000001 -130.4466283766202 +} + #=== HBO @test "init HBO STO-3G" { run_init HBO.xyz "-b STO-3G" hbo.ezfio diff --git a/tests/input/dhno.xyz b/tests/input/dhno.xyz new file mode 100644 index 00000000..367a2fa7 --- /dev/null +++ b/tests/input/dhno.xyz @@ -0,0 +1,7 @@ +4 +XYZ file: coordinates in Angstrom +H -0.877367 -1.047049 0.000000 +N 0.000000 -0.544985 0.000000 +O 0.000000 0.738624 0.000000 +H 0.877367 -1.047049 0.000000 + From 6bbc24c1ecab453b79698dd1c63cf33cf7bebb22 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 16 Feb 2016 17:43:26 +0100 Subject: [PATCH 27/29] Added the diagonalize_s2 option --- .travis.yml | 2 +- config/ifort.cfg | 2 +- plugins/Full_CI/full_ci.irp.f | 16 +++ src/Determinants/EZFIO.cfg | 6 + src/Determinants/diagonalize_CI.irp.f | 191 +++++++++++++++++++++++--- src/Determinants/s2.irp.f | 171 +++++++++++++++++++++++ src/Determinants/slater_rules.irp.f | 27 ++++ 7 files changed, 394 insertions(+), 21 deletions(-) diff --git a/.travis.yml b/.travis.yml index f451f1d6..b5cf3053 100644 --- a/.travis.yml +++ b/.travis.yml @@ -24,7 +24,7 @@ python: script: - ./configure --production ./config/gfortran.cfg - - source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD + - source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles - source ./quantum_package.rc ; ninja - source ./quantum_package.rc ; cd ocaml ; make ; cd - - source ./quantum_package.rc ; cd tests ; bats bats/qp.bats diff --git a/config/ifort.cfg b/config/ifort.cfg index cc848cba..f8c36e4e 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [OPTION] -MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below CACHE : 1 ; Enable cache_compile.py OPENMP : 1 ; Append OpenMP flags diff --git a/plugins/Full_CI/full_ci.irp.f b/plugins/Full_CI/full_ci.irp.f index 074ec7e1..56bf96d9 100644 --- a/plugins/Full_CI/full_ci.irp.f +++ b/plugins/Full_CI/full_ci.irp.f @@ -62,11 +62,27 @@ program full_ci endif print *, 'N_det = ', N_det print *, 'N_states = ', N_states + do k = 1, N_states + print*,'State ',k print *, 'PT2 = ', pt2 print *, 'E = ', CI_energy print *, 'E(before)+PT2 = ', E_CI_before+pt2 + enddo print *, '-----' E_CI_before = CI_energy + if(N_states.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_states + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_states + print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1)) + enddo + endif + E_CI_before = CI_energy call ezfio_set_full_ci_energy(CI_energy) if (abort_all) then exit diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index 4ab84b7a..b1c459ba 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2 interface: ezfio,provider,ocaml default: False +[diagonalize_s2] +type: logical +doc: Diagonalize the S^2 operator within the n_states_diag states required. Notice : the vectors are sorted by increasing S^2 values. +interface: ezfio,provider,ocaml +default: True + [threshold_davidson] type: Threshold doc: Thresholds of Davidson's algorithm diff --git a/src/Determinants/diagonalize_CI.irp.f b/src/Determinants/diagonalize_CI.irp.f index 9af8d413..b533bed2 100644 --- a/src/Determinants/diagonalize_CI.irp.f +++ b/src/Determinants/diagonalize_CI.irp.f @@ -36,17 +36,36 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ] - implicit none BEGIN_DOC ! Eigenvectors/values of the CI matrix END_DOC - integer :: i,j + implicit none + double precision :: ovrlp,u_dot_v + integer :: i_good_state + integer, allocatable :: index_good_state_array(:) + logical, allocatable :: good_state_array(:) + double precision, allocatable :: s2_values_tmp(:) + integer :: i_other_state + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + integer :: i_state + double precision :: s2,e_0 + integer :: i,j,k + double precision, allocatable :: s2_eigvalues(:) + double precision, allocatable :: e_array(:) + integer, allocatable :: iorder(:) - do j=1,N_states_diag + ! Guess values for the "N_states_diag" states of the CI_eigenvectors + do j=1,min(N_states_diag,N_det) do i=1,N_det CI_eigenvectors(i,j) = psi_coef(i,j) enddo enddo + + do j=N_det+1,N_states_diag + do i=1,N_det + CI_eigenvectors(i,j) = 0.d0 + enddo + enddo if (diag_algorithm == "Davidson") then @@ -59,36 +78,77 @@ END_PROVIDER else if (diag_algorithm == "Lapack") then - double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) allocate (eigenvalues(N_det)) call lapack_diag(eigenvalues,eigenvectors, & H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) CI_electronic_energy(:) = 0.d0 - do i=1,N_det - CI_eigenvectors(i,1) = eigenvectors(i,1) - enddo - integer :: i_state - double precision :: s2 if (s2_eig) then i_state = 0 + allocate (s2_eigvalues(N_det)) + allocate(index_good_state_array(N_det),good_state_array(N_det)) + good_state_array = .False. do j=1,N_det call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2) - print*,'s2 = ',s2 + s2_eigvalues(j) = s2 + ! Select at least n_states states with S^2 values closed to "expected_s2" if(dabs(s2-expected_s2).le.0.3d0)then - i_state += 1 - do i=1,N_det - CI_eigenvectors(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy(i_state) = eigenvalues(j) - CI_eigenvectors_s2(i_state) = s2 + i_state +=1 + index_good_state_array(i_state) = j + good_state_array(j) = .True. endif - if (i_state.ge.N_states_diag) then - exit + if(i_state.eq.N_states) then + exit endif enddo + if(i_state .ne.0)then + ! Fill the first "i_state" states that have a correct S^2 value + do j = 1, i_state + do i=1,N_det + CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) + enddo + CI_electronic_energy(j) = eigenvalues(index_good_state_array(j)) + CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) + enddo + i_other_state = 0 + do j = 1, N_det + if(good_state_array(j))cycle + i_other_state +=1 + if(i_state+i_other_state.gt.n_states_diag)then + exit + endif + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2) + do i=1,N_det + CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) + enddo + CI_electronic_energy(i_state+i_other_state) = eigenvalues(j) + CI_eigenvectors_s2(i_state+i_other_state) = s2 + enddo + + deallocate(index_good_state_array,good_state_array) + + else + print*,'' + print*,'!!!!!!!! WARNING !!!!!!!!!' + print*,' Within the ',N_det,'determinants selected' + print*,' and the ',N_states_diag,'states requested' + print*,' We did not find any state with S^2 values close to ',expected_s2 + print*,' We will then set the first N_states eigenvectors of the H matrix' + print*,' as the CI_eigenvectors' + print*,' You should consider more states and maybe ask for diagonalize_s2 to be .True. or just enlarge the CI space' + print*,'' + do j=1,min(N_states_diag,N_det) + do i=1,N_det + CI_eigenvectors(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy(j) = eigenvalues(j) + CI_eigenvectors_s2(j) = s2_eigvalues(j) + enddo + endif + deallocate(s2_eigvalues) else - do j=1,N_states_diag + ! Select the "N_states_diag" states of lowest energy + do j=1,min(N_det,N_states_diag) call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) do i=1,N_det CI_eigenvectors(i,j) = eigenvectors(i,j) @@ -99,7 +159,100 @@ END_PROVIDER endif deallocate(eigenvectors,eigenvalues) endif + + if(diagonalize_s2.and.n_states_diag > 1.and. n_det >= n_states_diag)then + ! Diagonalizing S^2 within the "n_states_diag" states found + allocate(s2_eigvalues(N_states_diag)) + call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors,n_det,size(psi_det,3),size(CI_eigenvectors,1),min(n_states_diag,n_det),s2_eigvalues) + + do j = 1, N_states_diag + do i = 1, N_det + psi_coef(i,j) = CI_eigenvectors(i,j) + enddo + enddo + + if(s2_eig)then + + ! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value + ! closer to the "expected_s2" set as input + + allocate(index_good_state_array(N_det),good_state_array(N_det)) + good_state_array = .False. + i_state = 0 + do j = 1, N_states_diag + if(dabs(s2_eigvalues(j)-expected_s2).le.0.3d0)then + good_state_array(j) = .True. + i_state +=1 + index_good_state_array(i_state) = j + endif + enddo + ! Sorting the i_state good states by energy + allocate(e_array(i_state),iorder(i_state)) + do j = 1, i_state + do i = 1, N_det + CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(j)) + enddo + CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) + call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int) + CI_electronic_energy(j) = e_0 + e_array(j) = e_0 + iorder(j) = j + enddo + call dsort(e_array,iorder,i_state) + do j = 1, i_state + CI_electronic_energy(j) = e_array(j) + CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(iorder(j))) + do i = 1, N_det + CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(iorder(j))) + enddo +! call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int) +! print*,'e = ',CI_electronic_energy(j) +! print*,' = ',e_0 +! call get_s2_u0(psi_det,CI_eigenvectors(1,j),N_det,size(CI_eigenvectors,1),s2) +! print*,'s^2 = ',CI_eigenvectors_s2(j) +! print*,'= ',s2 + enddo + deallocate(e_array,iorder) + + ! Then setting the other states without any specific energy order + i_other_state = 0 + do j = 1, N_states_diag + if(good_state_array(j))cycle + i_other_state +=1 + do i = 1, N_det + CI_eigenvectors(i,i_state + i_other_state) = psi_coef(i,j) + enddo + CI_eigenvectors_s2(i_state + i_other_state) = s2_eigvalues(j) + call u0_H_u_0(e_0,CI_eigenvectors(1,i_state + i_other_state),n_det,psi_det,N_int) + CI_electronic_energy(i_state + i_other_state) = e_0 + enddo + deallocate(index_good_state_array,good_state_array) + + + else + + ! Sorting the N_states_diag by energy, whatever the S^2 value is + + allocate(e_array(n_states_diag),iorder(n_states_diag)) + do j = 1, N_states_diag + call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int) + e_array(j) = e_0 + iorder(j) = j + enddo + call dsort(e_array,iorder,n_states_diag) + do j = 1, N_states_diag + CI_electronic_energy(j) = e_array(j) + do i = 1, N_det + CI_eigenvectors(i,j) = psi_coef(i,iorder(j)) + enddo + CI_eigenvectors_s2(j) = s2_eigvalues(iorder(j)) + enddo + deallocate(e_array,iorder) + endif + deallocate(s2_eigvalues) + endif + END_PROVIDER subroutine diagonalize_CI diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index 6da7b8ec..d84e4578 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -214,4 +214,175 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) deallocate (shortcut, sort_idx, sorted, version) end +subroutine get_uJ_s2_uI(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nstates) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax_keys) + integer, intent(in) :: n,nmax_coefs,nmax_keys,nstates + double precision, intent(in) :: psi_coefs_tmp(nmax_coefs,nstates) + double precision, intent(out) :: s2(nstates,nstates) + double precision :: s2_tmp,accu + integer :: i,j,l,jj,ll,kk + integer, allocatable :: idx(:) + double precision, allocatable :: tmp(:,:) + BEGIN_DOC + ! returns the matrix elements of S^2 "s2(i,j)" between the "nstates" states + ! psi_coefs_tmp(:,i) and psi_coefs_tmp(:,j) + END_DOC + s2 = 0.d0 + do ll = 1, nstates + do jj = 1, nstates + accu = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i,j,kk,idx,tmp,s2_tmp) & + !$OMP SHARED (ll,jj,psi_keys_tmp,psi_coefs_tmp,N_int,n,nstates) & + !$OMP REDUCTION(+:accu) + allocate(idx(0:n)) + !$OMP DO SCHEDULE(dynamic) + do i = 1, n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) + accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(i,jj) + call filter_connected(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx) + do kk=1,idx(0) + j = idx(kk) + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int) + accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(j,jj) + psi_coefs_tmp(i,jj) * s2_tmp * psi_coefs_tmp(j,ll) + enddo + enddo + !$OMP END DO NOWAIT + deallocate(idx) + !$OMP BARRIER + !$OMP END PARALLEL + s2(ll,jj) += accu + enddo + enddo + do i = 1, nstates + do j =i+1,nstates + accu = 0.5d0 * (s2(i,j) + s2(j,i)) + s2(i,j) = accu + s2(j,i) = accu + enddo + enddo +end + +subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nmax_coefs,nstates,s2_eigvalues) + BEGIN_DOC +! You enter with nstates vectors in psi_coefs_inout that may be coupled by S^2 +! The subroutine diagonalize the S^2 operator in the basis of these states. +! The vectors that you obtain in output are no more coupled by S^2, +! which does not necessary mean that they are eigenfunction of S^2. +! n,nmax,nstates = number of determinants, physical dimension of the arrays and number of states +! keys_tmp = array of integer(bit_kind) that represents the determinants +! psi_coefs(i,j) = coeff of the ith determinant in the jth state +! VECTORS ARE SUPPOSED TO BE ORTHONORMAL IN INPUT + END_DOC + implicit none + use bitmasks + integer, intent(in) :: n,nmax_keys,nmax_coefs,nstates + integer(bit_kind), intent(in) :: keys_tmp(N_int,2,nmax_keys) + double precision, intent(inout) :: psi_coefs_inout(nmax_coefs,nstates) + +!integer, intent(in) :: ndets_real,ndets_keys,ndets_coefs,nstates +!integer(bit_kind), intent(in) :: keys_tmp(N_int,2,ndets_keys) +!double precision, intent(inout) :: psi_coefs_inout(ndets_coefs,nstates) + double precision, intent(out) :: s2_eigvalues(nstates) + + + double precision,allocatable :: s2(:,:),overlap(:,:) + double precision, allocatable :: eigvalues(:),eigvectors(:,:) + integer :: i,j,k + double precision, allocatable :: psi_coefs_tmp(:,:) + double precision :: accu,coef_contract + double precision :: u_dot_u,u_dot_v + + print*,'' + print*,'*********************************************************************' + print*,'Cleaning the various vectors by diagonalization of the S^2 matrix ...' + print*,'' + print*,'nstates = ',nstates + allocate(s2(nstates,nstates),overlap(nstates,nstates)) + do i = 1, nstates + overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n) + do j = i+1, nstates + overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n) + overlap(j,i) = overlap(i,j) + enddo + enddo + print*,'Overlap matrix in the basis of the states considered' + do i = 1, nstates + write(*,'(10(F16.10,X))')overlap(i,:) + enddo + call ortho_lowdin(overlap,size(overlap,1),nstates,psi_coefs_inout,size(psi_coefs_inout,1),n) + print*,'passed ortho' + + do i = 1, nstates + overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n) + do j = i+1, nstates + overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n) + overlap(j,i) = overlap(i,j) + enddo + enddo + print*,'Overlap matrix in the basis of the Lowdin orthonormalized states ' + do i = 1, nstates + write(*,'(10(F16.10,X))')overlap(i,:) + enddo + + call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates) + print*,'S^2 matrix in the basis of the states considered' + double precision :: accu_precision_diag,accu_precision_of_diag + accu_precision_diag = 0.d0 + accu_precision_of_diag = 0.d0 + do i = 1, nstates + do j = i+1, nstates + if( ( dabs(s2(i,i) - s2(j,j)) .le.1.d-10 ) .and. (dabs(s2(i,j) + dabs(s2(i,j)))) .le.1.d-10) then + s2(i,j) = 0.d0 + s2(j,i) = 0.d0 + endif + enddo + enddo + do i = 1, nstates + write(*,'(10(F10.6,X))')s2(i,:) + enddo + + print*,'Diagonalizing the S^2 matrix' + + allocate(eigvalues(nstates),eigvectors(nstates,nstates)) + call lapack_diagd(eigvalues,eigvectors,s2,nstates,nstates) + print*,'Eigenvalues of s^2' + do i = 1, nstates + print*,'s2 = ',eigvalues(i) + s2_eigvalues(i) = eigvalues(i) + enddo + + print*,'Building the eigenvectors of the S^2 matrix' + allocate(psi_coefs_tmp(nmax_coefs,nstates)) + psi_coefs_tmp = 0.d0 + do j = 1, nstates + do k = 1, nstates + coef_contract = eigvectors(k,j) ! + do i = 1, n_det + psi_coefs_tmp(i,j) += psi_coefs_inout(i,k) * coef_contract + enddo + enddo + enddo + do j = 1, nstates + accu = 0.d0 + do i = 1, n_det + accu += psi_coefs_tmp(i,j) * psi_coefs_tmp(i,j) + enddo + print*,'Norm of vector = ',accu + accu = 1.d0/dsqrt(accu) + do i = 1, n_det + psi_coefs_inout(i,j) = psi_coefs_tmp(i,j) * accu + enddo + enddo +!call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates) +!print*,'S^2 matrix in the basis of the NEW states considered' +!do i = 1, nstates +! write(*,'(10(F16.10,X))')s2(i,:) +!enddo + + deallocate(s2,eigvalues,eigvectors,psi_coefs_tmp,overlap) + +end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 202d310d..e983ec34 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1507,6 +1507,33 @@ subroutine get_occ_from_key(key,occ,Nint) end +subroutine u0_H_u_0(e_0,u_0,n,keys_tmp,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes e_0 = / + ! + ! n : number of determinants + ! + END_DOC + integer, intent(in) :: n,Nint + double precision, intent(out) :: e_0 + double precision, intent(in) :: u_0(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + + double precision :: H_jj(n) + double precision :: v_0(n) + double precision :: u_dot_u,u_dot_v,diag_H_mat_elem + integer :: i,j + do i = 1, n + H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint) + enddo + + call H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) + e_0 = u_dot_v(v_0,u_0,n)/u_dot_u(u_0,n) +end + + subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) use bitmasks implicit none From 19e276fc0d064cd62f973b9a3675aae399074f42 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 16 Feb 2016 18:32:53 +0100 Subject: [PATCH 28/29] added the mulliken and hyperfine coupling constants analysis --- plugins/Hartree_Fock/.gitignore | 2 - plugins/Properties/hyperfine_constants.irp.f | 135 +++++++++++++++++++ plugins/Properties/mulliken.irp.f | 107 +++++++++++++++ plugins/Properties/print_mulliken.irp.f | 35 +++++ scripts/generate_h_apply.py | 4 +- src/AO_Basis/aos.irp.f | 27 ++++ src/AO_Basis/aos_value.irp.f | 48 +++++++ src/Determinants/density_matrix.irp.f | 51 +++++++ src/Integrals_Bielec/.gitignore | 3 +- src/Integrals_Monoelec/.gitignore | 2 - src/MOGuess/.gitignore | 5 +- 11 files changed, 407 insertions(+), 12 deletions(-) create mode 100644 plugins/Properties/hyperfine_constants.irp.f create mode 100644 plugins/Properties/mulliken.irp.f create mode 100644 plugins/Properties/print_mulliken.irp.f create mode 100644 src/AO_Basis/aos_value.irp.f diff --git a/plugins/Hartree_Fock/.gitignore b/plugins/Hartree_Fock/.gitignore index 9f1c0929..f1a4ff4f 100644 --- a/plugins/Hartree_Fock/.gitignore +++ b/plugins/Hartree_Fock/.gitignore @@ -5,7 +5,6 @@ AO_Basis Bitmask Electrons Ezfio_files -Huckel_guess IRPF90_man IRPF90_temp Integrals_Bielec @@ -16,7 +15,6 @@ Makefile Makefile.depend Nuclei Pseudo -SCF Utils ZMQ ezfio_interface.irp.f diff --git a/plugins/Properties/hyperfine_constants.irp.f b/plugins/Properties/hyperfine_constants.irp.f new file mode 100644 index 00000000..c1d88d2c --- /dev/null +++ b/plugins/Properties/hyperfine_constants.irp.f @@ -0,0 +1,135 @@ +BEGIN_PROVIDER [double precision, spin_density_at_nucleous, (nucl_num)] + implicit none + BEGIN_DOC +! value of the spin density at each nucleus + END_DOC + integer :: i,j,k + do i = 1, nucl_num + double precision :: r(3),accu,aos_array(ao_num) + accu = 0.d0 + r(1:3) = nucl_coord(i,1:3) + call give_all_aos_at_r(r,aos_array) + do j = 1, ao_num + do k = 1, ao_num + accu += one_body_spin_density_ao(k,j) * aos_array(k) * aos_array(j) + enddo + enddo + spin_density_at_nucleous(i) = accu + enddo +END_PROVIDER + + BEGIN_PROVIDER [double precision, spin_density_at_nucleous_from_mo, (nucl_num)] +&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_per_mo, (nucl_num,mo_tot_num)] + implicit none + BEGIN_DOC +! value of the spin density at each nucleus + END_DOC + integer :: i,j,k,l,m + do i = 1, nucl_num + double precision :: r(3),accu,aos_array(ao_num) + double precision :: contrib + double precision :: mo_values(mo_tot_num) + accu = 0.d0 + r(1:3) = nucl_coord(i,1:3) + call give_all_aos_at_r(r,aos_array) + spin_density_at_nucleous_from_mo(i) = 0.d0 + do k = 1, mo_tot_num + mo_values(k) = 0.d0 + do j = 1, ao_num + mo_values(k) += mo_coef(j,k) * aos_array(j) + enddo + enddo + do k = 1, mo_tot_num + spin_density_at_nucleous_contrib_per_mo(i,k) = 0.d0 + do m = 1, mo_tot_num + contrib = one_body_spin_density_mo(k,m) * mo_values(k) * mo_values(m) + spin_density_at_nucleous_from_mo(i) += contrib + spin_density_at_nucleous_contrib_per_mo(i,k) += contrib + enddo + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo, (nucl_num,mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo_test, (nucl_num)] + implicit none + BEGIN_DOC +! value of the spin density at each nucleus + END_DOC + integer :: i,j,k,l,m + spin_density_at_nucleous_contrib_mo_test = 0.d0 + do i = 1, nucl_num + double precision :: r(3),accu,aos_array(ao_num) + double precision :: c_i1,c_j1 + r(1:3) = nucl_coord(i,1:3) + call give_all_aos_at_r(r,aos_array) + do k = 1, mo_tot_num + do m = 1, mo_tot_num + accu = 0.d0 + do j = 1, ao_num + c_i1 = mo_coef(j,k) + do l = 1, ao_num + c_j1 = c_i1*mo_coef(l,m) + accu += one_body_spin_density_mo(k,m) * aos_array(l) * aos_array(j) * c_j1 + enddo + enddo + spin_density_at_nucleous_contrib_mo(i,k,m) = accu + spin_density_at_nucleous_contrib_mo_test(i) += accu + enddo + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [double precision, conversion_factor_mhz_hcc, (100)] +&BEGIN_PROVIDER [double precision, conversion_factor_gauss_hcc, (100)] +&BEGIN_PROVIDER [double precision, conversion_factor_cm_1_hcc, (100)] + BEGIN_DOC +! Conversion factor for the calculation of the hcc, according to the nuclear charge + END_DOC + + conversion_factor_mhz_hcc =0.d0 + conversion_factor_mhz_hcc =0.d0 + conversion_factor_mhz_hcc =0.d0 + + + ! hydrogen + conversion_factor_mhz_hcc(1) = 4469.84692227102460d0 + conversion_factor_gauss_hcc(1) = 1594.95296390862904d0 + conversion_factor_cm_1_hcc(1) = 1490.98044430157870d0 + + ! Li + conversion_factor_mhz_hcc(3) = 1737.2746512855997d0 + conversion_factor_gauss_hcc(3) = 619.9027742370165d0 + conversion_factor_cm_1_hcc(3) = 579.4924475562677d0 + + ! carbon + conversion_factor_mhz_hcc(6) = 1124.18303629792945d0 + conversion_factor_gauss_hcc(6) = 401.136570647523058d0 + conversion_factor_cm_1_hcc(6) = 374.987097339830086d0 + + ! nitrogen + conversion_factor_mhz_hcc(7) = 323.102093833793390d0 + conversion_factor_gauss_hcc(7) = 115.290892768082614d0 + conversion_factor_cm_1_hcc(7) = 107.775257586297698d0 + + ! Oxygen + conversion_factor_mhz_hcc(8) = -606.1958551736545d0 + conversion_factor_gauss_hcc(8) = -216.30574771560407d0 + conversion_factor_cm_1_hcc(8) = -202.20517197179822d0 + +END_PROVIDER + + BEGIN_PROVIDER [double precision, iso_hcc_mhz, (nucl_num)] +&BEGIN_PROVIDER [double precision, iso_hcc_gauss, (nucl_num)] +&BEGIN_PROVIDER [double precision, iso_hcc_cm_1, (nucl_num)] + BEGIN_DOC +! isotropic hyperfine coupling constants among the various atoms + END_DOC + integer :: i + do i = 1, nucl_num + iso_hcc_mhz(i) = conversion_factor_mhz_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !* 0.5d0 + iso_hcc_gauss(i) = conversion_factor_gauss_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i)!* 0.5d0 + iso_hcc_cm_1(i) = conversion_factor_cm_1_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !*0.5d0 + enddo + +END_PROVIDER diff --git a/plugins/Properties/mulliken.irp.f b/plugins/Properties/mulliken.irp.f new file mode 100644 index 00000000..d56c9a44 --- /dev/null +++ b/plugins/Properties/mulliken.irp.f @@ -0,0 +1,107 @@ + +BEGIN_PROVIDER [double precision, spin_population, (ao_num_align,ao_num)] + implicit none + integer :: i,j + BEGIN_DOC +! spin population on the ao basis : +! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * + END_DOC + spin_population = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + spin_population(j,i) = one_body_spin_density_ao(i,j) * ao_overlap(i,j) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)] + implicit none + integer :: i + double precision :: accu + spin_population_angular_momentum = 0.d0 + do i = 1, ao_num + spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, spin_gross_orbital_product, (ao_num)] + implicit none + spin_gross_orbital_product = 0.d0 + integer :: i,j + BEGIN_DOC +! gross orbital product for the spin population + END_DOC + do i = 1, ao_num + do j = 1, ao_num + spin_gross_orbital_product(i) += spin_population(j,i) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [double precision, mulliken_spin_densities, (nucl_num)] + implicit none + integer :: i,j + BEGIN_DOC +!ATOMIC SPIN POPULATION (ALPHA MINUS BETA) + END_DOC + mulliken_spin_densities = 0.d0 + do i = 1, ao_num + mulliken_spin_densities(ao_nucl(i)) += spin_gross_orbital_product(i) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, electronic_population_alpha, (ao_num_align,ao_num)] +&BEGIN_PROVIDER [double precision, electronic_population_beta, (ao_num_align,ao_num)] + implicit none + integer :: i,j + BEGIN_DOC +! spin population on the ao basis : +! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * + END_DOC + electronic_population_alpha = 0.d0 + electronic_population_beta = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + electronic_population_alpha(j,i) = one_body_dm_ao_alpha(i,j) * ao_overlap(i,j) + electronic_population_beta(j,i) = one_body_dm_ao_beta(i,j) * ao_overlap(i,j) + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, gross_orbital_product_alpha, (ao_num)] +&BEGIN_PROVIDER [double precision, gross_orbital_product_beta, (ao_num)] + implicit none + spin_gross_orbital_product = 0.d0 + integer :: i,j + BEGIN_DOC +! gross orbital product + END_DOC + do i = 1, ao_num + do j = 1, ao_num + gross_orbital_product_alpha(i) += electronic_population_alpha(j,i) + gross_orbital_product_beta(i) += electronic_population_beta(j,i) + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, mulliken_densities_alpha, (nucl_num)] +&BEGIN_PROVIDER [double precision, mulliken_densities_beta, (nucl_num)] + implicit none + integer :: i,j + BEGIN_DOC +! + END_DOC + mulliken_densities_alpha = 0.d0 + mulliken_densities_beta = 0.d0 + do i = 1, ao_num + mulliken_densities_alpha(ao_nucl(i)) += gross_orbital_product_alpha(i) + mulliken_densities_beta(ao_nucl(i)) += gross_orbital_product_beta(i) + enddo + +END_PROVIDER diff --git a/plugins/Properties/print_mulliken.irp.f b/plugins/Properties/print_mulliken.irp.f new file mode 100644 index 00000000..100c8556 --- /dev/null +++ b/plugins/Properties/print_mulliken.irp.f @@ -0,0 +1,35 @@ +program print_mulliken + implicit none + read_wf = .True. + touch read_wf + print*,'Mulliken spin densities' + + call test +end +subroutine test + double precision :: accu + integer :: i + integer :: j + accu= 0.d0 + do i = 1, nucl_num + print*,i,nucl_charge(i),mulliken_spin_densities(i) + accu += mulliken_spin_densities(i) + enddo + print*,'Sum of Mulliken SD = ',accu + print*,'AO SPIN POPULATIONS' + accu = 0.d0 + do i = 1, ao_num + accu += spin_gross_orbital_product(i) + write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i) + enddo + print*,'sum = ',accu + accu = 0.d0 + print*,'Angular momentum analysis' + do i = 0, ao_l_max + accu += spin_population_angular_momentum(i) + print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) + print*,'sum = ',accu + enddo + +end + diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index b52e4445..f83d8b41 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -156,11 +156,11 @@ class H_apply(object): def filter_only_1h1p(self): self["filter_only_1h1p_single"] = """ ! ! DIR$ FORCEINLINE - if (is_a_1h1p(hole).eq..False.) cycle + if (is_a_1h1p(hole).eqv..False.) cycle """ self["filter_only_1h1p_double"] = """ ! ! DIR$ FORCEINLINE - if (is_a_1h1p(key).eq..False.) cycle + if (is_a_1h1p(key).eqv..False.) cycle """ diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index 341d1453..8c2db90e 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -146,3 +146,30 @@ integer function ao_power_index(nx,ny,nz) ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1 end + BEGIN_PROVIDER [ integer, ao_l, (ao_num) ] +&BEGIN_PROVIDER [ integer, ao_l_max ] +&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ] + implicit none + BEGIN_DOC +! ao_l = l value of the AO: a+b+c in x^a y^b z^c + END_DOC + integer :: i + do i=1,ao_num + ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) + ao_l_char(i) = l_to_charater(ao_l(i)) + enddo + ao_l_max = maxval(ao_l) +END_PROVIDER + +BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)] + BEGIN_DOC + ! character corresponding to the "L" value of an AO orbital + END_DOC + implicit none + l_to_charater(0)='S' + l_to_charater(1)='P' + l_to_charater(2)='D' + l_to_charater(3)='F' + l_to_charater(4)='G' +END_PROVIDER + diff --git a/src/AO_Basis/aos_value.irp.f b/src/AO_Basis/aos_value.irp.f new file mode 100644 index 00000000..a531ce50 --- /dev/null +++ b/src/AO_Basis/aos_value.irp.f @@ -0,0 +1,48 @@ +double precision function ao_value(i,r) + implicit none + BEGIN_DOC +! return the value of the ith ao at point r + END_DOC + double precision, intent(in) :: r(3) + integer, intent(in) :: i + + integer :: m,num_ao + double precision :: center_ao(3) + double precision :: beta + integer :: power_ao(3) + num_ao = ao_nucl(i) + power_ao(1:3)= ao_power(i,1:3) + center_ao(1:3) = nucl_coord(num_ao,1:3) + double precision :: accu,dx,dy,dz,r2 + dx = (r(1) - center_ao(1)) + dy = (r(2) - center_ao(2)) + dz = (r(3) - center_ao(3)) + r2 = dx*dx + dy*dy + dz*dz + dx = dx**power_ao(1) + dy = dy**power_ao(2) + dz = dz**power_ao(3) + + accu = 0.d0 + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) + enddo + ao_value = accu * dx * dy * dz + +end + +subroutine give_all_aos_at_r(r,aos_array) + implicit none + BEGIN_dOC +! gives the values of aos at a given point r + END_DOC + double precision, intent(in) :: r(3) + double precision, intent(out) :: aos_array(ao_num) + integer :: i + double precision :: ao_value + do i = 1, ao_num + aos_array(i) = ao_value(i,r) + enddo + + +end diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index e5d243f4..62d09381 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -206,3 +206,54 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] state_average_weight = 1.d0/dble(N_states) END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_body_spin_density_ao, (ao_num_align,ao_num) ] + BEGIN_DOC +! one body spin density matrix on the AO basis : rho_AO(alpha) - rho_AO(beta) + END_DOC + implicit none + integer :: i,j,k,l + double precision :: dm_mo + + one_body_spin_density_ao = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, mo_tot_num + do j = 1, mo_tot_num + dm_mo = one_body_spin_density_mo(j,i) +! if(dabs(dm_mo).le.1.d-10)cycle + one_body_spin_density_ao(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo + + enddo + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha, (ao_num_align,ao_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta, (ao_num_align,ao_num) ] + BEGIN_DOC +! one body density matrix on the AO basis : rho_AO(alpha) , rho_AO(beta) + END_DOC + implicit none + integer :: i,j,k,l + double precision :: dm_mo + + one_body_spin_density_ao = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, mo_tot_num + do j = 1, mo_tot_num + dm_mo = one_body_dm_mo_alpha(j,i) +! if(dabs(dm_mo).le.1.d-10)cycle + one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo + one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo + + enddo + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/Integrals_Bielec/.gitignore b/src/Integrals_Bielec/.gitignore index ad6b465d..1d52a821 100644 --- a/src/Integrals_Bielec/.gitignore +++ b/src/Integrals_Bielec/.gitignore @@ -17,5 +17,4 @@ ZMQ ezfio_interface.irp.f irpf90.make irpf90_entities -tags -test_integrals \ No newline at end of file +tags \ No newline at end of file diff --git a/src/Integrals_Monoelec/.gitignore b/src/Integrals_Monoelec/.gitignore index 577068de..e8bd9b05 100644 --- a/src/Integrals_Monoelec/.gitignore +++ b/src/Integrals_Monoelec/.gitignore @@ -12,9 +12,7 @@ Makefile.depend Nuclei Pseudo Utils -check_orthonormality ezfio_interface.irp.f irpf90.make irpf90_entities -save_ortho_mos tags \ No newline at end of file diff --git a/src/MOGuess/.gitignore b/src/MOGuess/.gitignore index e9ba5cf5..797574f4 100644 --- a/src/MOGuess/.gitignore +++ b/src/MOGuess/.gitignore @@ -4,7 +4,6 @@ AO_Basis Electrons Ezfio_files -H_CORE_guess IRPF90_man IRPF90_temp Integrals_Monoelec @@ -15,8 +14,6 @@ Nuclei Pseudo Utils ezfio_interface.irp.f -guess_overlap irpf90.make irpf90_entities -tags -truncate_mos \ No newline at end of file +tags \ No newline at end of file From 27a121bc1ba960bf1e778af2388e8c46b756e9e4 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 17 Feb 2016 10:52:57 +0100 Subject: [PATCH 29/29] added print_hcc and OVB plugin --- config/ifort.cfg | 2 +- plugins/Hartree_Fock/EZFIO.cfg | 7 + plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES | 2 +- plugins/Hartree_Fock/diagonalize_fock.irp.f | 74 +-- plugins/OVB/NEEDED_CHILDREN_MODULES | 1 + plugins/OVB/README.rst | 20 + plugins/OVB/ovb_components.irp.f | 510 +++++++++++++++++++ plugins/OVB/print_ovb.irp.f | 27 + plugins/Properties/print_hcc.irp.f | 17 + plugins/Psiref_Utils/psi_ref_utils.irp.f | 21 +- src/Bitmask/bitmasks.irp.f | 61 ++- src/Determinants/usefull_for_ovb.irp.f | 283 ++++++++++ 12 files changed, 967 insertions(+), 58 deletions(-) create mode 100644 plugins/OVB/NEEDED_CHILDREN_MODULES create mode 100644 plugins/OVB/README.rst create mode 100644 plugins/OVB/ovb_components.irp.f create mode 100644 plugins/OVB/print_ovb.irp.f create mode 100644 plugins/Properties/print_hcc.irp.f create mode 100644 src/Determinants/usefull_for_ovb.irp.f diff --git a/config/ifort.cfg b/config/ifort.cfg index f8c36e4e..cc848cba 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [OPTION] -MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below CACHE : 1 ; Enable cache_compile.py OPENMP : 1 ; Append OpenMP flags diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index d8207cc4..2fa29cf0 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -26,3 +26,10 @@ default: Huckel type: double precision doc: Calculated HF energy interface: ezfio + +[no_oa_or_av_opt] +type: logical +doc: If true, skip the (inactive+core) --> (active) and the (active) --> (virtual) orbital rotations within the SCF procedure +interface: ezfio,provider,ocaml +default: False + diff --git a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES index 85bdd3ad..6fb87e35 100644 --- a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES +++ b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Integrals_Bielec MOGuess +Integrals_Bielec MOGuess Bitmask diff --git a/plugins/Hartree_Fock/diagonalize_fock.irp.f b/plugins/Hartree_Fock/diagonalize_fock.irp.f index 7c424aeb..c80077b3 100644 --- a/plugins/Hartree_Fock/diagonalize_fock.irp.f +++ b/plugins/Hartree_Fock/diagonalize_fock.irp.f @@ -11,63 +11,35 @@ double precision, allocatable :: work(:), F(:,:), S(:,:) -! if (mo_tot_num == ao_num) then -! ! Solve H.C = E.S.C in AO basis set -! -! allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) ) -! do j=1,ao_num -! do i=1,ao_num -! S(i,j) = ao_overlap(i,j) -! F(i,j) = Fock_matrix_ao(i,j) -! enddo -! enddo -! -! n = ao_num -! lwork = 1+6*n + 2*n*n -! liwork = 3 + 5*n -! -! allocate(work(lwork), iwork(liwork) ) -! -! lwork = -1 -! liwork = -1 -! -! call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& -! diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) -! -! if (info /= 0) then -! print *, irp_here//' failed : ', info -! stop 1 -! endif -! lwork = int(work(1)) -! liwork = iwork(1) -! deallocate(work,iwork) -! allocate(work(lwork), iwork(liwork) ) -! -! call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),& -! diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) -! -! if (info /= 0) then -! print *, irp_here//' failed : ', info -! stop 1 -! endif -! do j=1,mo_tot_num -! do i=1,ao_num -! eigenvectors_Fock_matrix_mo(i,j) = F(i,j) -! enddo -! enddo -! -! deallocate(work, iwork, F, S) -! -! else -! - ! Solve H.C = E.C in MO basis set - allocate( F(mo_tot_num_align,mo_tot_num) ) do j=1,mo_tot_num do i=1,mo_tot_num F(i,j) = Fock_matrix_mo(i,j) enddo enddo + if(no_oa_or_av_opt)then + integer :: iorb,jorb + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + enddo + endif + + ! Insert level shift here diff --git a/plugins/OVB/NEEDED_CHILDREN_MODULES b/plugins/OVB/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..05706cac --- /dev/null +++ b/plugins/OVB/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants Psiref_CAS diff --git a/plugins/OVB/README.rst b/plugins/OVB/README.rst new file mode 100644 index 00000000..13488c7e --- /dev/null +++ b/plugins/OVB/README.rst @@ -0,0 +1,20 @@ +======================= +OVB +======================= +The present module proposes an orthogonal Valence Bond analysis +of the wave function, that are the printing of the various Hamiltonian +matrix elements on the basis of the level of ionicity of the components +of the wave function. + +Assumptions : it supposes that you have some orthogonal local orbitals within +the active space and that you performed a CI within the active orbitals. +Such CI might be complete or not, no matter. + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/OVB/ovb_components.irp.f b/plugins/OVB/ovb_components.irp.f new file mode 100644 index 00000000..82f0f931 --- /dev/null +++ b/plugins/OVB/ovb_components.irp.f @@ -0,0 +1,510 @@ + + use bitmasks + BEGIN_PROVIDER [integer, max_number_ionic] +&BEGIN_PROVIDER [integer, min_number_ionic] + BEGIN_DOC + ! Maximum and minimum number of ionization in psi_ref + END_DOC + implicit none + integer :: i,j + integer :: n_closed_shell_cas + max_number_ionic = 0 + min_number_ionic = 100000 + do i = 1, N_det_ref + j = n_closed_shell_cas(psi_ref(1,1,i),n_int) + if(j> max_number_ionic)then + max_number_ionic = j + endif + if(j< min_number_ionic)then + min_number_ionic = j + endif + + enddo + print*,'max_number_ionic = ',max_number_ionic + print*,'min_number_ionic = ',min_number_ionic +END_PROVIDER + + BEGIN_PROVIDER [integer, ionic_index, (min_number_ionic:max_number_ionic,0:N_det_ref)] +&BEGIN_PROVIDER [double precision, normalization_factor_ionic, (min_number_ionic:max_number_ionic, N_states)] + BEGIN_DOC + ! Index of the various determinants in psi_ref according to their level of ionicity + ! ionic_index(i,0) = number of determinants in psi_ref having the degree of ionicity "i" + ! ionic_index(i,j) = index of the determinants having the degree of ionicity "i" + END_DOC + implicit none + integer :: i,j,k + integer :: n_closed_shell_cas + double precision :: accu + ionic_index = 0 + do i = 1, N_det_ref + j = n_closed_shell_cas(psi_ref(1,1,i),n_int) + ionic_index(j,0) +=1 + ionic_index(j,ionic_index(j,0)) = i + enddo + do i = min_number_ionic,max_number_ionic + accu = 0.d0 + do j = 1, N_states + do k = 1, ionic_index(i,0) + accu += psi_ref_coef_diagonalized(ionic_index(i,k),j) * psi_ref_coef_diagonalized(ionic_index(i,k),j) + enddo + normalization_factor_ionic(i,j) = 1.d0/dsqrt(accu) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [double precision, H_OVB_naked, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)] + BEGIN_DOC + ! Hamiltonian matrix expressed in the basis of contracted forms in terms of ionic structures + END_DOC + implicit none + integer :: i,j,istate,k,l + double precision :: accu,hij + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + do istate = 1, N_states + accu = 0.d0 + do k = 1, ionic_index(i,0) + do l = 1, ionic_index(j,0) + hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l)) + accu += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_ionic(i,istate) * & + psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + enddo + enddo + H_OVB_naked(i,j,istate) = accu + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [integer, n_couples_act_orb] + implicit none + n_couples_act_orb = 3 + END_PROVIDER + + BEGIN_PROVIDER [integer, couples_act_orb, (n_couples_act_orb,2) ] + implicit none + + couples_act_orb(1,1) = 20 + couples_act_orb(1,2) = 21 + couples_act_orb(2,1) = 22 + couples_act_orb(2,2) = 23 + couples_act_orb(3,1) = 24 + couples_act_orb(3,2) = 25 + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom , (n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC +! Hamiltonian matrix elements between the various contracted functions +! that have a negative charge on a given active orbital + END_DOC + integer :: i,j,k,l,jj,ii + integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:) + double precision :: accu,hij + double precision :: norm + allocate (key_1(N_int,2),key_2(N_int,2)) + do i = 1, n_act_orb + j = i ! Diagonal part + norm = 0.d0 + accu = 0.d0 + do k = 1, n_det_ionic_on_given_atom(i) + norm += psi_coef_mono_ionic_on_given_atom(k,i) **2 + do ii = 1, N_int + key_1(ii,1) = psi_det_mono_ionic_on_given_atom(ii,1,k,i) + key_1(ii,2) = psi_det_mono_ionic_on_given_atom(ii,2,k,i) + enddo + do l = 1, n_det_ionic_on_given_atom(j) + do jj = 1, N_int + key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j) + key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij + enddo + enddo + H_matrix_between_ionic_on_given_atom(i,j) = accu + + + do j = i+1, n_act_orb ! Extra diagonal part + accu = 0.d0 + do k = 1, n_det_ionic_on_given_atom(i) + do jj = 1, N_int + key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i) + key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i) + enddo + do l = 1, n_det_ionic_on_given_atom(j) + do jj = 1, N_int + key_2(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,l,j) + key_2(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_mono_ionic_on_given_atom(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij + enddo + enddo + H_matrix_between_ionic_on_given_atom(i,j) = accu + H_matrix_between_ionic_on_given_atom(j,i) = accu + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, H_matrix_between_ionic_on_given_atom_and_others , (n_act_orb,min_number_ionic:max_number_ionic)] + implicit none + use bitmasks + BEGIN_DOC +! Hamiltonian matrix elements between the various contracted functions +! that have a negative charge on a given active orbital +! and all the other fully contracted OVB structures + END_DOC + integer :: i,j,k,l,jj,ii + integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:) + double precision :: accu,hij + double precision :: norm + allocate (key_1(N_int,2),key_2(N_int,2)) + + do i = 1, n_act_orb + do j = min_number_ionic,max_number_ionic + if(j==1)then + H_matrix_between_ionic_on_given_atom_and_others(i,j) = 0.d0 + endif + accu = 0.d0 + do k = 1, n_det_ionic_on_given_atom(i) + do jj = 1, N_int + key_1(jj,1) = psi_det_mono_ionic_on_given_atom(jj,1,k,i) + key_1(jj,2) = psi_det_mono_ionic_on_given_atom(jj,2,k,i) + enddo + do l = 1, ionic_index(j,0) + do ii = 1, N_int + key_2(ii,1) = psi_det_ovb(ii,1,l,j) + key_2(ii,2) = psi_det_ovb(ii,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_ovb(l,j) * psi_coef_mono_ionic_on_given_atom(k,i) * hij + enddo + enddo + H_matrix_between_ionic_on_given_atom_and_others(i,j) = accu + enddo + enddo + + print*,'H_matrix_between_ionic_on_given_atom_and_others' + print*,'' + do i = 1, n_act_orb + write(*,'(I3,X,100(F16.7))'),H_matrix_between_ionic_on_given_atom_and_others(i,:) + enddo + + + +END_PROVIDER + + BEGIN_PROVIDER [integer, n_det_ionic_on_given_atom, (n_act_orb)] +&BEGIN_PROVIDER [double precision, normalization_factor_ionic_on_given_atom, (n_act_orb) ] +&BEGIN_PROVIDER [double precision, psi_coef_mono_ionic_on_given_atom, (N_det_ref,n_act_orb) ] +&BEGIN_PROVIDER [integer(bit_kind), psi_det_mono_ionic_on_given_atom, (N_int,2,N_det_ref,n_act_orb)] + implicit none + use bitmasks + BEGIN_DOC +! number of determinants that are mono ionic with the negative charge +! on a given atom, normalization_factor, array of determinants,and coefficients + END_DOC + integer :: i,j,k,l + ionicity_level = 1 + integer :: ionicity_level + logical :: doubly_occupied_array(n_act_orb) + n_det_ionic_on_given_atom = 0 + normalization_factor_ionic_on_given_atom = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + call give_index_of_doubly_occ_in_active_space(psi_det(1,1,ionic_index(ionicity_level,i)),doubly_occupied_array) + do j = 1, n_act_orb + if(doubly_occupied_array(j))then + n_det_ionic_on_given_atom(j) += 1 + normalization_factor_ionic_on_given_atom(j) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2 + do k = 1, N_int + psi_det_mono_ionic_on_given_atom(k,1,n_det_ionic_on_given_atom(j),j) = psi_det(k,1,ionic_index(ionicity_level,i)) + psi_det_mono_ionic_on_given_atom(k,2,n_det_ionic_on_given_atom(j),j) = psi_det(k,2,ionic_index(ionicity_level,i)) + enddo + psi_coef_mono_ionic_on_given_atom(n_det_ionic_on_given_atom(j),j) = psi_ref_coef_diagonalized(ionic_index(1,i),1) + endif + enddo + enddo + integer :: i_count + i_count = 0 + do j = 1, n_act_orb + i_count += n_det_ionic_on_given_atom(j) + normalization_factor_ionic_on_given_atom(j) = 1.d0/dsqrt(normalization_factor_ionic_on_given_atom(j)) + enddo + if(i_count.ne.ionic_index(ionicity_level,0))then + print*,'PB with n_det_ionic_on_given_atom' + print*,'i_count = ',i_count + print*,'ionic_index(ionicity_level,0)',ionic_index(ionicity_level,0) + stop + endif + do j = 1, n_act_orb + do i = 1, n_det_ionic_on_given_atom(j) + psi_coef_mono_ionic_on_given_atom(i,j) = psi_coef_mono_ionic_on_given_atom(i,j) * normalization_factor_ionic_on_given_atom(j) + enddo + enddo + + + END_PROVIDER + + BEGIN_PROVIDER [integer(bit_kind), psi_det_ovb, (N_int,2,N_det_ref,min_number_ionic:max_number_ionic)] +&BEGIN_PROVIDER [double precision, psi_coef_ovb, (N_det_ref,min_number_ionic:max_number_ionic) ] + implicit none + BEGIN_DOC +! Array of the determinants belonging to each ovb structures (neutral, mono ionic, bi ionic etc ...) +! together with the arrays of coefficients + END_DOC + integer :: i,j,k,l + use bitmasks + integer :: ionicity_level,i_count + double precision :: accu + + do ionicity_level = min_number_ionic,max_number_ionic + accu = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + do j = 1, N_int + psi_det_ovb(j,1,i,ionicity_level) = psi_det(j,1,ionic_index(ionicity_level,i)) + psi_det_ovb(j,2,i,ionicity_level) = psi_det(j,2,ionic_index(ionicity_level,i)) + enddo + psi_coef_ovb(i,ionicity_level) = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) * normalization_factor_ionic(ionicity_level,1) + accu += psi_coef_ovb(i,ionicity_level)**2 + enddo + accu = 1.d0/dsqrt(accu) + do i = 1, ionic_index(ionicity_level,0) + psi_coef_ovb(i,ionicity_level) = psi_coef_ovb(i,ionicity_level) * accu + enddo + accu = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + accu += psi_coef_ovb(i,ionicity_level) **2 + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, H_matrix_psi_det_ovb, (min_number_ionic:max_number_ionic,min_number_ionic:max_number_ionic)] + implicit none + BEGIN_DOC +! H matrix between the fully contracted OVB forms + END_DOC + integer :: i,j,k,l,jj,ii + integer(bit_kind), allocatable :: key_1(:,:),key_2(:,:) + use bitmasks + double precision :: accu,hij + double precision :: norm + allocate (key_1(N_int,2),key_2(N_int,2)) + do i = min_number_ionic,max_number_ionic + do j = min_number_ionic,max_number_ionic + accu = 0.d0 + do k = 1, ionic_index(i,0) + do ii = 1, N_int + key_1(ii,1) = psi_det_ovb(ii,1,k,i) + key_1(ii,2) = psi_det_ovb(ii,2,k,i) + enddo + do l = 1, ionic_index(j,0) + do ii = 1, N_int + key_2(ii,1) = psi_det_ovb(ii,1,l,j) + key_2(ii,2) = psi_det_ovb(ii,2,l,j) + enddo + call i_H_j(key_1,key_2,N_int,hij) + accu += psi_coef_ovb(l,j) * psi_coef_ovb(k,i) * hij + enddo + enddo + H_matrix_psi_det_ovb(i,j) = accu + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [integer, number_first_ionic_couples] +&BEGIN_PROVIDER [logical , is_a_first_ionic_couple, (N_det_ref)] +&BEGIN_PROVIDER [double precision, normalization_factor_special_first_ionic, (2)] + implicit none + BEGIN_DOC + ! Number of determinants belonging to the class of first ionic + ! AND that have a couple of positive/negative charge belonging + ! to a couple of orbital couples_act_orb + ! If is_a_first_ionic_couple(i) = .True. then this determinant is a first ionic + ! and have a couple of positive/negative charge belonging + ! to a couple of orbital couples_act_orb + ! normalization factor (1) = 1/(sum c_i^2 .with. is_a_first_ionic_couple(i) = .True.) + ! normalization factor (2) = 1/(sum c_i^2 .with. is_a_first_ionic_couple(i) = .False.) + END_DOC + integer :: i,j + use bitmasks + number_first_ionic_couples = 0 + integer :: ionicity_level + logical :: couples_out(0:n_couples_act_orb) + integer(bit_kind) :: key_tmp(N_int,2) + ionicity_level = 1 + normalization_factor_special_first_ionic = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + do j = 1, N_int + key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i)) + key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i)) + enddo + call doubly_occ_empty_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out) + if(couples_out(0))then + number_first_ionic_couples +=1 + is_a_first_ionic_couple(i) = .True. + normalization_factor_special_first_ionic(1) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2 + else + is_a_first_ionic_couple(i) = .False. + normalization_factor_special_first_ionic(2) += psi_ref_coef_diagonalized(ionic_index(1,i),1) **2 + endif + enddo + normalization_factor_special_first_ionic(1) = 1.d0/dsqrt(normalization_factor_special_first_ionic(1)) + normalization_factor_special_first_ionic(2) = 1.d0/dsqrt(normalization_factor_special_first_ionic(2)) + print*,'number_first_ionic_couples = ',number_first_ionic_couples + END_PROVIDER + + + BEGIN_PROVIDER [integer, number_neutral_no_hund_couples] +&BEGIN_PROVIDER [logical , is_a_neutral_no_hund_couple, (N_det_ref)] +&BEGIN_PROVIDER [double precision, normalization_factor_neutra_no_hund_couple, (2)] +&BEGIN_PROVIDER [double precision, ratio_hund_no_hund ] + implicit none + BEGIN_DOC + ! Number of determinants belonging to the class of neutral determinants + ! AND that have a couple of alpha beta electrons in couple of orbital couples_act_orb + ! If is_a_neutral_no_hund_couple(i) = .True. then this determinant is a neutral determinants + ! and have a a couple of alpha beta electrons in couple of orbital couples_act_orb + ! normalization factor (1) = 1/sqrt(sum c_i^2 .with. is_a_neutral_no_hund_couple(i) = .True.) + ! normalization factor (2) = 1/sqrt(sum c_i^2 .with. is_a_neutral_no_hund_couple(i) = .False.) + END_DOC + integer :: i,j + use bitmasks + number_neutral_no_hund_couples = 0 + integer :: ionicity_level + logical :: couples_out(0:n_couples_act_orb) + integer(bit_kind) :: key_tmp(N_int,2) + integer :: ifirst_hund,ifirst_no_hund + double precision :: coef_ref_hund,coef_ref_no_hund + ifirst_hund = 0 + ifirst_no_hund = 0 + ionicity_level = 0 + normalization_factor_neutra_no_hund_couple = 0.d0 + do i = 1, ionic_index(ionicity_level,0) + do j = 1, N_int + key_tmp(j,1) = psi_det(j,1,ionic_index(ionicity_level,i)) + key_tmp(j,2) = psi_det(j,2,ionic_index(ionicity_level,i)) + enddo + call neutral_no_hund_in_couple(key_tmp,n_couples_act_orb,couples_act_orb,couples_out) + if(couples_out(0))then + if(ifirst_no_hund == 0)then + coef_ref_no_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) + ifirst_no_hund = 1 + endif + number_neutral_no_hund_couples +=1 + is_a_neutral_no_hund_couple(i) = .True. + normalization_factor_neutra_no_hund_couple(1) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2 + else + if(ifirst_hund == 0)then + coef_ref_hund = psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) + ifirst_hund = 1 + endif + is_a_neutral_no_hund_couple(i) = .False. + normalization_factor_neutra_no_hund_couple(2) += psi_ref_coef_diagonalized(ionic_index(ionicity_level,i),1) **2 + endif + enddo + ratio_hund_no_hund = coef_ref_no_hund/coef_ref_hund + + normalization_factor_neutra_no_hund_couple(1) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(1)) + normalization_factor_neutra_no_hund_couple(2) = 1.d0/dsqrt(normalization_factor_neutra_no_hund_couple(2)) + print*,'number_neutral_no_hund_couples = ',number_neutral_no_hund_couples + END_PROVIDER + + BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic, (2,min_number_ionic:max_number_ionic,n_states)] +&BEGIN_PROVIDER [double precision, H_OVB_naked_first_ionic_between_ionic, (2,2,n_states)] + BEGIN_DOC + ! H_OVB_naked_first_ionic(1,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True + ! and the contracted ith ionic form + ! if i == 1 not defined + ! H_OVB_naked_first_ionic(2,i) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False + ! and the contracted ith ionic form + ! if i == 1 not defined + ! H_OVB_naked_first_ionic_between_ionic(1,1) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True + ! and the first ionic determinants belonging to is_a_first_ionic_couple = True + ! H_OVB_naked_first_ionic_between_ionic(1,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = True + ! and the first ionic determinants belonging to is_a_first_ionic_couple = False + ! H_OVB_naked_first_ionic_between_ionic(2,2) = H_matrix element between the first ionic determinants belonging to is_a_first_ionic_couple = False + ! and the first ionic determinants belonging to is_a_first_ionic_couple = False + END_DOC + implicit none + integer :: i,j,istate,k,l + double precision :: accu_1,accu_2,hij + H_OVB_naked_first_ionic = 0.d0 + H_OVB_naked_first_ionic_between_ionic = 0.d0 + i = 1 + do j = min_number_ionic,max_number_ionic + if(j==1)cycle + do istate = 1, N_states + accu_1 = 0.d0 + accu_2 = 0.d0 + do k = 1, ionic_index(i,0) + if(is_a_first_ionic_couple(k))then + do l = 1, ionic_index(j,0) + hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l)) + accu_1 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(1) * & + psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + enddo + else + do l = 1, ionic_index(j,0) + hij = ref_hamiltonian_matrix(ionic_index(i,k),ionic_index(j,l)) + accu_2 += psi_ref_coef_diagonalized(ionic_index(i,k),istate) * normalization_factor_special_first_ionic(2) * & + psi_ref_coef_diagonalized(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij + enddo + endif + enddo + H_OVB_naked_first_ionic(1,j,istate) = accu_1 + H_OVB_naked_first_ionic(2,j,istate) = accu_2 + enddo + enddo + + + do istate = 1, N_states + accu_1 = 0.d0 + accu_2 = 0.d0 + integer :: i_count + i_count = 0 + do k = 1, ionic_index(1,0) + do l = 1, ionic_index(1,0) + hij = ref_hamiltonian_matrix(ionic_index(1,k),ionic_index(1,l)) + accu_1 = hij * psi_ref_coef_diagonalized(ionic_index(1,k),istate) * psi_ref_coef_diagonalized(ionic_index(1,l),istate) + if(is_a_first_ionic_couple(k).and. is_a_first_ionic_couple(l))then + H_OVB_naked_first_ionic_between_ionic(1,1,istate) += accu_1 * normalization_factor_special_first_ionic(1) **2 + elseif(is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then + i_count += 1 + H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 * & + normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2) +! elseif(is_a_first_ionic_couple(l).and. .not.is_a_first_ionic_couple(k))then +! i_count += 1 +! H_OVB_naked_first_ionic_between_ionic(1,2,istate) += accu_1 * & +! normalization_factor_special_first_ionic(1) *normalization_factor_special_first_ionic(2) + elseif(.not.is_a_first_ionic_couple(k).and. .not.is_a_first_ionic_couple(l))then + H_OVB_naked_first_ionic_between_ionic(2,2,istate) += accu_1 * normalization_factor_special_first_ionic(2) **2 + endif + enddo + enddo + enddo + print*,'i_count = ',i_count + print*,'number_first_ionic_couples**2 = ',ionic_index(1,0) * number_first_ionic_couples + + double precision :: convert_hartree_ev + convert_hartree_ev = 27.211399d0 + print*,'Special H matrix' + do i = 1,2 + write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic(i,:,1) + enddo + + print*,'Special H matrix bis' + do i = 1,2 + write(*,'(I4,X,10(F16.8 ,4X))')i, H_OVB_naked_first_ionic_between_ionic(i,:,1) + enddo + + + END_PROVIDER + diff --git a/plugins/OVB/print_ovb.irp.f b/plugins/OVB/print_ovb.irp.f new file mode 100644 index 00000000..9e333c15 --- /dev/null +++ b/plugins/OVB/print_ovb.irp.f @@ -0,0 +1,27 @@ +program print_OVB + implicit none + read_wf = .True. + call provide_all + +end + +subroutine provide_all + implicit none + integer :: i,j,k,l,istate + do istate= 1, N_states + print*,'-------------------' + print*,'ISTATE = ',istate + print*,'-------------------' + print*,'CAS MATRIX ' + print*,'' + do i = min_number_ionic,max_number_ionic + write(*,'(I4,X,10(F8.5 ,4X))')i, H_OVB_naked(i,:,istate) + enddo + print*,'' + print*,'-------------------' + print*,'-------------------' + enddo + + +end + diff --git a/plugins/Properties/print_hcc.irp.f b/plugins/Properties/print_hcc.irp.f new file mode 100644 index 00000000..f0091e1e --- /dev/null +++ b/plugins/Properties/print_hcc.irp.f @@ -0,0 +1,17 @@ +program print_hcc + implicit none + read_wf = .True. + touch read_wf + call test +end +subroutine test + implicit none + double precision :: accu + integer :: i,j + print*,'Z AU GAUSS MHZ cm^-1' + do i = 1, nucl_num + write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i) + enddo + +end + diff --git a/plugins/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f index f9cf1303..fb45b13d 100644 --- a/plugins/Psiref_Utils/psi_ref_utils.irp.f +++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f @@ -125,7 +125,7 @@ BEGIN_PROVIDER [double precision, H_matrix_ref, (N_det_ref,N_det_ref)] enddo END_PROVIDER - BEGIN_PROVIDER [double precision, psi_coef_ref_diagonalized, (N_det_ref,N_states)] + BEGIN_PROVIDER [double precision, psi_ref_coef_diagonalized, (N_det_ref,N_states)] &BEGIN_PROVIDER [double precision, psi_ref_energy_diagonalized, (N_states)] implicit none integer :: i,j @@ -137,9 +137,11 @@ END_PROVIDER do i = 1, N_states psi_ref_energy_diagonalized(i) = eigenvalues(i) do j = 1, N_det_ref - psi_coef_ref_diagonalized(j,i) = eigenvectors(j,i) + psi_ref_coef_diagonalized(j,i) = eigenvectors(j,i) enddo enddo + deallocate (eigenvectors) + deallocate (eigenvalues) END_PROVIDER @@ -264,3 +266,18 @@ integer function get_index_in_psi_ref_sorted_bit(key,Nint) end +BEGIN_PROVIDER [double precision, ref_hamiltonian_matrix, (n_det_ref,n_det_ref)] + BEGIN_DOC + ! H matrix in the Reference space + END_DOC + implicit none + integer :: i,j + double precision :: hij + do i = 1, N_det_ref + do j = 1, N_det_ref + call i_H_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij) + ref_hamiltonian_matrix(i,j) = hij + enddo + enddo +END_PROVIDER + diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 70a84a42..50a95312 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -289,7 +289,12 @@ END_PROVIDER &BEGIN_PROVIDER [ integer, n_virt_orb ] implicit none BEGIN_DOC - ! Bitmasks for the inactive orbitals that are excited in post CAS method + ! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited + ! in post CAS methods + ! n_inact_orb : Number of inactive orbitals + ! virt_bitmask : Bitmaks of vritual orbitals which are supposed to be recieve electrons + ! in post CAS methods + ! n_virt_orb : Number of virtual orbitals END_DOC logical :: exists integer :: j,i @@ -327,8 +332,14 @@ END_PROVIDER - BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)] + BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)] &BEGIN_PROVIDER [ integer, list_virt, (n_virt_orb)] + BEGIN_DOC + ! list_inact : List of the inactive orbitals which are supposed to be doubly excited + ! in post CAS methods + ! list_virt : List of vritual orbitals which are supposed to be recieve electrons + ! in post CAS methods + END_DOC implicit none integer :: occ_inact(N_int*bit_kind_size) integer :: itest,i @@ -348,6 +359,21 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the inactive, active and virtual bitmasks + END_DOC + integer :: i,j + do i = 1, N_int + reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1)) + reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2)) + enddo + END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] implicit none BEGIN_DOC @@ -376,7 +402,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)] implicit none BEGIN_DOC - ! Reunion of the inactive, active and virtual bitmasks + ! Bitmask of the core orbitals that are never excited in post CAS method END_DOC integer :: i,j do i = 1, N_int @@ -385,6 +411,35 @@ END_PROVIDER enddo END_PROVIDER + BEGIN_PROVIDER [integer, list_core, (n_core_orb)] + BEGIN_DOC + ! List of the core orbitals that are never excited in post CAS method + END_DOC + implicit none + integer :: occ_core(N_int*bit_kind_size) + integer :: itest,i + occ_core = 0 + call bitstring_to_list(core_bitmask(1,1), occ_core(1), itest, N_int) + ASSERT(itest==n_core_orb) + do i = 1, n_core_orb + list_core(i) = occ_core(i) + enddo + END_PROVIDER + + BEGIN_PROVIDER [ integer, n_core_orb ] + implicit none + BEGIN_DOC + ! Number of core orbitals that are never excited in post CAS method + END_DOC + logical :: exists + integer :: j,i + integer :: i_hole,i_part,i_gen + + n_core_orb = 0 + do j = 1, N_int + n_core_orb += popcnt(core_bitmask(j,1)) + enddo + END_PROVIDER BEGIN_PROVIDER [ integer, i_bitmask_gen ] diff --git a/src/Determinants/usefull_for_ovb.irp.f b/src/Determinants/usefull_for_ovb.irp.f new file mode 100644 index 00000000..7b89897b --- /dev/null +++ b/src/Determinants/usefull_for_ovb.irp.f @@ -0,0 +1,283 @@ + +integer function n_open_shell(det_in,nint) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(nint,2),nint + integer :: i + n_open_shell = 0 + do i=1,Nint + n_open_shell += popcnt(iand(xor(det_in(i,1),det_in(i,2)),det_in(i,1))) + enddo +end + +integer function n_closed_shell(det_in,nint) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(nint,2),nint + integer :: i + n_closed_shell = 0 + do i=1,Nint + n_closed_shell += popcnt(iand(det_in(i,1),det_in(i,2))) + enddo +end + +integer function n_closed_shell_cas(det_in,nint) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(nint,2),nint + integer(bit_kind) :: det_tmp(nint,2) + integer :: i + n_closed_shell_cas = 0 + do i=1,Nint + det_tmp(i,1) = xor(det_in(i,1),reunion_of_core_inact_bitmask(i,1)) + det_tmp(i,2) = xor(det_in(i,2),reunion_of_core_inact_bitmask(i,2)) + enddo +!call debug_det(det_tmp,nint) + do i=1,Nint + n_closed_shell_cas += popcnt(iand(det_tmp(i,1),det_tmp(i,2))) + enddo +end + +subroutine doubly_occ_empty_in_couple(det_in,n_couples,couples,couples_out) + implicit none + use bitmasks + integer, intent(in) :: n_couples,couples(n_couples,2) + integer(bit_kind),intent(in) :: det_in(N_int,2) + logical, intent(out) :: couples_out(0:n_couples) + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + ! n_couples is the number of couples of orbitals to be checked + ! couples(i,1) = first orbital of the ith couple + ! couples(i,2) = second orbital of the ith couple + ! returns the array couples_out + ! couples_out(i) = .True. if det_in contains + ! an orbital empty in the ith couple AND + ! an orbital doubly occupied in the ith couple + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + couples_out(0) = .False. + do i = 1, n_couples + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb + ! det_tmp is zero except for the two orbitals of the couple + integer :: i_count + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 1)then + couples_out(i) = .False. + cycle + endif + + ! test if orbital there are two electrons or not + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j))) + enddo + if(i_count.ne.1)then + couples_out(i) = .False. + else + couples_out(i) = .True. + couples_out(0) = .True. + endif + enddo +end + +subroutine give_index_of_doubly_occ_in_active_space(det_in,doubly_occupied_array) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: det_in(N_int,2) + logical, intent(out) :: doubly_occupied_array(n_act_orb) + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + do i = 1, n_act_orb + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + i_bite = list_act(i) + call set_bit_to_integer(i_bite,det_tmp_bis,N_int) ! act orb + ! det_tmp is zero except for the orbital "ith" active orbital + integer :: i_count,i_bite + + ! test if orbital there are two electrons or not + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j))) + enddo + if(i_count.ne.1)then + doubly_occupied_array(i) = .False. + else + doubly_occupied_array(i) = .True. + endif + enddo +end + +subroutine doubly_occ_empty_in_couple_and_no_hund_elsewhere(det_in,n_couple_no_hund,couple_ion,couple_no_hund,is_ok) + implicit none + use bitmasks + integer, intent(in) :: n_couple_no_hund,couple_ion(2),couple_no_hund(n_couple_no_hund,2) + integer(bit_kind),intent(in) :: det_in(N_int,2) + logical, intent(out) :: is_ok + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + ! n_couples is the number of couples of orbitals to be checked + ! couples(i,1) = first orbital of the ith couple + ! couples(i,2) = second orbital of the ith couple + ! returns the array couples_out + ! couples_out(i) = .True. if det_in contains + ! an orbital empty in the ith couple AND + ! an orbital doubly occupied in the ith couple + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + is_ok = .False. + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couple_ion(1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couple_ion(2),det_tmp_bis,N_int) ! second orb + ! det_tmp is zero except for the two orbitals of the couple + integer :: i_count + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 1)then + is_ok = .False. + return + endif + + ! test if orbital there are two electrons or not + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j))) + enddo + if(i_count.ne.1)then + is_ok = .False. + return + else + do i = 1, n_couple_no_hund + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couple_no_hund (i,1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couple_no_hund (i,2),det_tmp_bis,N_int) ! second orb + ! det_tmp_bis is zero except for the two orbitals of the couple + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 2)then + is_ok = .False. + return + endif + ! test if orbital there are one alpha and one beta + integer :: i_count_alpha,i_count_beta + i_count_alpha = 0 + i_count_beta = 0 + do j = 1, N_int + i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j))) + i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j))) + enddo + if(i_count_alpha==1.and.i_count_beta==1)then + is_ok = .True. + else + is_ok = .False. + return + endif + enddo + is_ok = .True. + endif +end + + +subroutine neutral_no_hund_in_couple(det_in,n_couples,couples,couples_out) + implicit none + use bitmasks + integer, intent(in) :: n_couples,couples(n_couples,2) + integer(bit_kind),intent(in) :: det_in(N_int,2) + logical, intent(out) :: couples_out(0:n_couples) + integer(bit_kind) :: det_tmp(N_int) + integer(bit_kind) :: det_tmp_bis(N_int) + BEGIN_DOC + ! n_couples is the number of couples of orbitals to be checked + ! couples(i,1) = first orbital of the ith couple + ! couples(i,2) = second orbital of the ith couple + ! returns the array couples_out + ! couples_out(i) = .True. if det_in contains + ! an orbital empty in the ith couple AND + ! an orbital doubly occupied in the ith couple + END_DOC + integer :: i,j,k,l + + ! det_tmp tells you if the orbitals are occupied or not + do j = 1, N_int + det_tmp(j) = ior(det_in(j,1),det_in(j,2)) + enddo + + couples_out(0) = .True. + do i = 1, n_couples + do j = 1, N_int + det_tmp_bis(j) = 0_bit_kind + enddo + call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb + call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb + ! det_tmp_bis is zero except for the two orbitals of the couple + integer :: i_count + i_count = 0 + do j = 1, N_int + i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied + enddo + if(i_count .ne. 2)then + couples_out(i) = .False. + cycle + endif + + ! test if orbital there are one alpha and one beta + integer :: i_count_alpha,i_count_beta + i_count_alpha = 0 + i_count_beta = 0 + do j = 1, N_int + i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j))) + i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j))) + enddo + if(i_count_alpha==1.and.i_count_beta==1)then + couples_out(i) = .True. + else + couples_out(i) = .False. + endif + enddo + do i = 1, n_couples + if(.not.couples_out(i))then + couples_out(0) = .False. + endif + enddo +end + +