LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity
2c8b68fada
quantum_package/scripts/generate_h_apply.py
Anthony Scemama ca973a1e92 Revert "Bugs to fix (#50)" (#51) 
This reverts commit 94f01c0892.
2017-04-20 08:45:56 +02:00

514 lines
16 KiB
Python
Executable File
Raw Blame History

#!/usr/bin/env python
import os
keywords = """
check_double_excitation
copy_buffer
declarations
decls_main
deinit_thread
skip
init_main
filter_integrals
filter2p
filter2h2p_double
filter2h2p_single
filter1h
filter1p
only_2p_single
only_2p_double
only_2h_single
only_2h_double
only_1h_single
only_1h_double
only_1p_single
only_1p_double
only_2h1p_single
only_2h1p_double
filter_only_1h1p_single
filter_only_1h1p_double
filter_only_1h2p_single
filter_only_1h2p_double
filter_only_2h2p_single
filter_only_2h2p_double
filterhole
filter_integrals
filter_only_1h1p_double
filter_only_1h1p_single
filterparticle
filter_vvvv_excitation
finalization
generate_psi_guess
initialization
init_main
init_thread
keys_work
omp_barrier
omp_do
omp_enddo
omp_end_master
omp_end_parallel
omp_master
omp_parallel
only_2p_double
only_2p_single
parameters
params_main
printout_always
printout_now
skip
subroutine
""".split()
class H_apply(object):
def read_template(self):
file = open(os.environ["QP_ROOT"]+'/src/Determinants/H_apply.template.f','r')
self.template = file.read()
file.close()
file = open(os.environ["QP_ROOT"]+'/src/Determinants/H_apply_nozmq.template.f','r')
self.template += file.read()
file.close()
def __init__(self,sub,SingleRef=False,do_mono_exc=True, do_double_exc=True):
self.read_template()
s = {}
for k in keywords:
s[k] = ""
s["subroutine"] = "H_apply_%s"%(sub)
s["params_post"] = ""
self.selection_pt2 = None
self.energy = "CI_electronic_energy"
self.perturbation = None
self.do_double_exc = do_double_exc
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
s["omp_parallel"] = """ PROVIDE elec_num_tab
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, &
!$OMP occ_particle,occ_hole,j_a,k_a,other_spin, &
!$OMP hole_save,ispin,jj,l_a,ib_jb_pairs,array_pairs, &
!$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, &
!$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,&
!$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, &
!$OMP N_elec_in_key_hole_2,ia_ja_pairs,key_union_hole_part) &
!$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, &
!$OMP hole_1, particl_1, hole_2, particl_2, &
!$OMP elec_alpha_num,i_generator) FIRSTPRIVATE(iproc)"""
s["omp_end_parallel"] = "!$OMP END PARALLEL"
s["omp_master"] = "!$OMP MASTER"
s["omp_end_master"] = "!$OMP END MASTER"
s["omp_barrier"] = "!$OMP BARRIER"
s["omp_do"] = "!$OMP DO SCHEDULE (static,1)"
s["omp_enddo"] = "!$OMP ENDDO NOWAIT"
d = { True : '.True.', False : '.False.'}
s["do_mono_excitations"] = d[do_mono_exc]
s["do_double_excitations"] = d[do_double_exc]
s["keys_work"] += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)"
s["filter_integrals"] = "array_pairs = .True."
if SingleRef:
s["filter_integrals"] = """
call get_mo_bielec_integrals_existing_ik(i_a,j_a,mo_tot_num,array_pairs,mo_integrals_map)
"""
s["generate_psi_guess"] = """
! Sort H_jj to find the N_states lowest states
integer :: i
integer, allocatable :: iorder(:)
double precision, allocatable :: H_jj(:)
double precision, external :: diag_h_mat_elem
allocate(H_jj(N_det),iorder(N_det))
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(psi_det,N_int,H_jj,iorder,N_det) &
!$OMP PRIVATE(i)
!$OMP DO
do i = 1, N_det
H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
iorder(i) = i
enddo
!$OMP END DO
!$OMP END PARALLEL
call dsort(H_jj,iorder,N_det)
do k=1,N_states
psi_coef(iorder(k),k) = 1.d0
enddo
deallocate(H_jj,iorder)
"""
s["size_max"] = "8192"
s["copy_buffer"] = """call copy_H_apply_buffer_to_wf
if (s2_eig) then
call make_s2_eigenfunction
endif
SOFT_TOUCH psi_det psi_coef N_det
"""
s["printout_now"] = """write(output_determinants,*) &
100.*float(i_generator)/float(N_det_generators), '% in ', wall_1-wall_0, 's'"""
self.data = s
def __setitem__(self,key,value):
self.data[key] = value
def __getitem__(self,key):
return self.data[key]
def __repr__(self):
buffer = self.template
for key,value in self.data.items():
buffer = buffer.replace('$'+key, value)
return buffer
def unset_double_excitations(self):
self["do_double_excitations"] = ".False."
self["check_double_excitation"] = """
check_double_excitation = .False.
"""
def filter_vvvv_excitation(self):
self["filter_vvvv_excitation"] = """
key_union_hole_part = 0_bit_kind
call set_bit_to_integer(i_a,key_union_hole_part,N_int)
call set_bit_to_integer(j_a,key_union_hole_part,N_int)
call set_bit_to_integer(i_b,key_union_hole_part,N_int)
call set_bit_to_integer(j_b,key_union_hole_part,N_int)
do jtest_vvvv = 1, N_int
if(iand(key_union_hole_part(jtest_vvvv),virt_bitmask(jtest_vvvv,1).ne.key_union_hole_part(jtest_vvvv)))then
b_cycle = .False.
endif
enddo
if(b_cycle) cycle
"""
def set_filter_holes(self):
self["filterhole"] = """
if(iand(ibset(0_bit_kind,j),hole(k,other_spin)).eq.0_bit_kind )cycle
"""
def set_filter_particl(self):
self["filterparticle"] = """
if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind )cycle
"""
def filter_1h(self):
self["filter1h"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h(hole)) cycle
"""
def filter_2p(self):
self["filter2p"] = """
! ! DIR$ FORCEINLINE
if (is_a_2p(hole)) cycle
"""
def filter_1p(self):
self["filter1p"] = """
! ! DIR$ FORCEINLINE
if (is_a_1p(hole)) cycle
"""
def filter_only_2h(self):
self["only_2h_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_2h(hole).eqv. .False.) cycle
"""
self["only_2h_double"] = """
! ! DIR$ FORCEINLINE
if ( is_a_2h(key).eqv. .False. )cycle
"""
def filter_only_1h(self):
self["only_1h_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h(hole) .eqv. .False.) cycle
"""
self["only_1h_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h(key) .eqv. .False.) cycle
"""
def filter_only_1p(self):
self["only_1p_single"] = """
! ! DIR$ FORCEINLINE
if ( is_a_1p(hole) .eqv. .False.) cycle
"""
self["only_1p_double"] = """
! ! DIR$ FORCEINLINE
if ( is_a_1p(key) .eqv. .False.) cycle
"""
def filter_only_2h1p(self):
self["only_2h1p_single"] = """
! ! DIR$ FORCEINLINE
if ( is_a_2h1p(hole) .eqv. .False.) cycle
"""
self["only_2h1p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_2h1p(key) .eqv. .False.) cycle
"""
def filter_only_2p(self):
self["only_2p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_2p(hole).eqv. .False.) cycle
"""
self["only_2p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_2p(key).eqv. .False.) cycle
"""
def filter_only_1h1p(self):
self["filter_only_1h1p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(hole).eqv..False.) cycle
"""
self["filter_only_1h1p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(key).eqv..False.) cycle
"""
def filter_only_2h2p(self):
self["filter_only_2h2p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_two_holes_two_particles(hole).eqv..False.) cycle
"""
self["filter_only_2h2p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_two_holes_two_particles(key).eqv..False.) cycle
"""
def filter_only_1h2p(self):
self["filter_only_1h2p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h2p(hole).eqv..False.) cycle
"""
self["filter_only_1h2p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h2p(key).eqv..False.) cycle
"""
def unset_skip(self):
self["skip"] = """
"""
def set_filter_2h_2p(self):
self["filter2h2p_double"] = """
if (is_a_two_holes_two_particles(key)) cycle
"""
self["filter2h2p_single"] = """
if (is_a_two_holes_two_particles(hole)) cycle
"""
def set_perturbation(self,pert):
if self.perturbation is not None:
raise
self.perturbation = pert
if pert is not None:
self.data["parameters"] = ",sum_e_2_pert_in,sum_norm_pert_in,sum_H_pert_diag_in,N_st,Nint"
self.data["declarations"] = """
integer, intent(in) :: N_st,Nint
double precision, intent(inout) :: sum_e_2_pert_in(N_st)
double precision, intent(inout) :: sum_norm_pert_in(N_st)
double precision, intent(inout) :: sum_H_pert_diag_in(N_st)
double precision :: sum_e_2_pert(N_st)
double precision :: sum_norm_pert(N_st)
double precision :: sum_H_pert_diag(N_st)
double precision, allocatable :: e_2_pert_buffer(:,:)
double precision, allocatable :: coef_pert_buffer(:,:)
ASSERT (Nint == N_int)
"""
self.data["init_thread"] = """
allocate (e_2_pert_buffer(N_st,size_max), coef_pert_buffer(N_st,size_max))
do k=1,N_st
sum_e_2_pert(k) = 0.d0
sum_norm_pert(k) = 0.d0
sum_H_pert_diag(k) = 0.d0
enddo
"""
self.data["deinit_thread"] = """
! OMP CRITICAL
do k=1,N_st
sum_e_2_pert_in(k) = sum_e_2_pert_in(k) + sum_e_2_pert(k)
sum_norm_pert_in(k) = sum_norm_pert_in(k) + sum_norm_pert(k)
sum_H_pert_diag_in(k) = sum_H_pert_diag_in(k) + sum_H_pert_diag(k)
enddo
! OMP END CRITICAL
deallocate (e_2_pert_buffer, coef_pert_buffer)
"""
self.data["size_max"] = "8192"
self.data["initialization"] = """
PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
"""
if self.do_double_exc == True:
self.data["keys_work"] = """
! if(check_double_excitation)then
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s)
"""%(pert,self.energy)
else:
self.data["keys_work"] = """
call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s)
"""%(pert,self.energy)
self.data["finalization"] = """
"""
self.data["copy_buffer"] = ""
self.data["generate_psi_guess"] = ""
self.data["params_main"] = "pt2, norm_pert, H_pert_diag, N_st"
self.data["params_post"] = ","+self.data["params_main"] +", N_int"
self.data["decls_main"] = """ integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag(N_st)
double precision :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st)
PROVIDE N_det_generators
do k=1,N_st
pt2(k) = 0.d0
norm_pert(k) = 0.d0
H_pert_diag(k) = 0.d0
norm_psi(k) = 0.d0
delta_pt2(k) = 0.d0
pt2_old(k) = 0.d0
enddo
write(output_determinants,'(A12, 1X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
'N_generators', 'Norm', 'Delta PT2', 'PT2', 'Est. PT2', 'secs'
write(output_determinants,'(A12, 1X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
'============', '========', '=========', '=========', '=========', &
'========='
"""
self.data["printout_always"] = """
do k=1,N_st
norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k)
delta_pt2(k) = pt2(k) - pt2_old(k)
enddo
"""
self.data["printout_now"] = """
do k=1,N_st
write(output_determinants,'(I10, 4(2X, F9.6), 2X, F8.1)') &
i_generator, norm_psi(k), delta_pt2(k), pt2(k), &
pt2(k)/(norm_psi(k)*norm_psi(k)), &
wall_1-wall_0
pt2_old(k) = pt2(k)
enddo
"""
self.data["omp_parallel"] += """&
!$OMP SHARED(N_st) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
!$OMP PRIVATE(sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
def set_selection_pt2(self,pert):
if self.selection_pt2 is not None:
raise
self.set_perturbation(pert)
self.selection_pt2 = pert
if pert is not None:
self.data["parameters"] += ",select_max_out"
self.data["declarations"] += """
double precision, intent(inout) :: select_max_out"""
self.data["params_post"] += ", select_max(min(i_generator,size(select_max,1)))"
self.data["size_max"] = "8192"
self.data["copy_buffer"] = """
call copy_H_apply_buffer_to_wf
if (s2_eig) then
call make_s2_eigenfunction
endif
SOFT_TOUCH psi_det psi_coef N_det
selection_criterion_min = min(selection_criterion_min, maxval(select_max))*0.1d0
selection_criterion = selection_criterion_min
call write_double(output_determinants,selection_criterion,'Selection criterion')
"""
self.data["keys_work"] = """
e_2_pert_buffer = 0.d0
coef_pert_buffer = 0.d0
""" + self.data["keys_work"]
self.data["keys_work"] += """
call fill_H_apply_buffer_selection(key_idx,keys_out,e_2_pert_buffer, &
coef_pert_buffer,N_st,N_int,iproc,select_max_out)
"""
self.data["omp_parallel"] += """&
!$OMP REDUCTION (max:select_max_out)"""
self.data["skip"] = """
if (i_generator < size_select_max) then
if (select_max(i_generator) < selection_criterion_min*selection_criterion_factor) then
! OMP CRITICAL
do k=1,N_st
norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k)
pt2_old(k) = 0.d0
enddo
! OMP END CRITICAL
cycle
endif
select_max(i_generator) = 0.d0
endif
"""
def unset_openmp(self):
for k in keywords:
if k.startswith("omp_"):
self[k] = ""
class H_apply_zmq(H_apply):
def read_template(self):
file = open(os.environ["QP_ROOT"]+'/src/Determinants/H_apply.template.f','r')
self.template = file.read()
file.close()
file = open(os.environ["QP_ROOT"]+'/src/Determinants/H_apply_zmq.template.f','r')
self.template += file.read()
file.close()
def set_perturbation(self,pert):
H_apply.set_perturbation(self,pert)
self.data["printout_now"] = ""
self.data["printout_always"] = ""
self.data["decls_main"] = """ integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag(N_st)
double precision :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st)
PROVIDE N_det_generators
do k=1,N_st
pt2(k) = 0.d0
norm_pert(k) = 0.d0
H_pert_diag(k) = 0.d0
norm_psi(k) = 0.d0
energy(k) = %s(k)
enddo
""" % (self.energy)
self.data["copy_buffer"] = """
do i=1,N_det_generators
do k=1,N_st
pt2(k) = pt2(k) + pt2_generators(k,i)
norm_pert(k) = norm_pert(k) + norm_pert_generators(k,i)
H_pert_diag(k) = H_pert_diag(k) + H_pert_diag_generators(k,i)
enddo
enddo
"""
def set_selection_pt2(self,pert):
H_apply.set_selection_pt2(self,pert)
self.data["skip"] = """
if (i_generator < size_select_max) then
if (select_max(i_generator) < selection_criterion_min*selection_criterion_factor) then
do k=1,N_st
pt2(k) = select_max(i_generator)
enddo
cycle
endif
select_max(i_generator) = 0.d0
endif
"""
Reference in New Issue View Git Blame Copy Permalink