mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-24 21:21:47 +02:00
Merge dev
This commit is contained in:
parent
5576bb6ccd
commit
bfc090ff16
3
.gitignore
vendored
3
.gitignore
vendored
@ -5,7 +5,10 @@ build.ninja
|
||||
.ninja_deps
|
||||
bin/
|
||||
lib/
|
||||
lib64/
|
||||
libexec/
|
||||
config/qp_create_ninja.pickle
|
||||
src/*/.gitignore
|
||||
ezfio_interface.irp.f
|
||||
share
|
||||
*.swp
|
||||
|
2
external/.gitignore
vendored
2
external/.gitignore
vendored
@ -1,2 +1,2 @@
|
||||
#*
|
||||
*
|
||||
|
||||
|
9
include/.gitignore
vendored
9
include/.gitignore
vendored
@ -1,8 +1 @@
|
||||
zmq.h
|
||||
gmp.h
|
||||
zconf.h
|
||||
zconf.h
|
||||
zlib.h
|
||||
zmq_utils.h
|
||||
f77_zmq_free.h
|
||||
f77_zmq.h
|
||||
*
|
||||
|
57
ocaml/.gitignore
vendored
57
ocaml/.gitignore
vendored
@ -1,57 +0,0 @@
|
||||
_build
|
||||
element_create_db
|
||||
element_create_db.byte
|
||||
ezfio.ml
|
||||
.gitignore
|
||||
Git.ml
|
||||
Input_ao_one_e_ints.ml
|
||||
Input_ao_two_e_erf_ints.ml
|
||||
Input_ao_two_e_ints.ml
|
||||
Input_auto_generated.ml
|
||||
Input_becke_numerical_grid.ml
|
||||
Input_cassd.ml
|
||||
Input_champ.ml
|
||||
Input_cipsi.ml
|
||||
Input_cipsi_tc_bi_ortho.ml
|
||||
Input_cosgtos_ao_int.ml
|
||||
Input_davidson_keywords.ml
|
||||
Input_davidson.ml
|
||||
Input_density_for_dft.ml
|
||||
Input_determinants.ml
|
||||
Input_dft_keywords.ml
|
||||
Input_dressing.ml
|
||||
Input_fci_tc_bi.ml
|
||||
Input_general_mrci.ml
|
||||
Input_mo_localization.ml
|
||||
Input_mo_one_e_ints.ml
|
||||
Input_mo_two_e_erf_ints.ml
|
||||
Input_mo_two_e_ints.ml
|
||||
Input_mpn.ml
|
||||
Input_mu_of_r.ml
|
||||
Input_nuclei.ml
|
||||
Input_perturbation.ml
|
||||
Input_pseudo.ml
|
||||
Input_qmcchem.ml
|
||||
Input_scf_utils.ml
|
||||
Input_some_mu_of_r.ml
|
||||
Input_tc_keywords.ml
|
||||
Input_three_body_ints.ml
|
||||
Input_trust_region.ml
|
||||
Input_two_body_rdm.ml
|
||||
Input_utils.ml
|
||||
Input_utils_trust_region.ml
|
||||
qp_create_ezfio
|
||||
qp_create_ezfio.native
|
||||
qp_edit
|
||||
qp_edit.ml
|
||||
qp_edit.native
|
||||
qp_print_basis
|
||||
qp_print_basis.native
|
||||
qp_run
|
||||
qp_run.native
|
||||
qp_set_mo_class
|
||||
qp_set_mo_class.native
|
||||
qp_tunnel
|
||||
qp_tunnel.native
|
||||
qptypes_generator.byte
|
||||
Qptypes.ml
|
@ -29,7 +29,7 @@ tests: $(ALL_TESTS)
|
||||
.gitignore: $(MLFILES) $(MLIFILES)
|
||||
@for i in .gitignore ezfio.ml element_create_db Qptypes.ml Git.ml *.byte *.native _build $(ALL_EXE) $(ALL_TESTS) \
|
||||
$(patsubst %.ml,%,$(wildcard test_*.ml)) $(patsubst %.ml,%,$(wildcard qp_*.ml)) \
|
||||
$(shell grep Input Input_auto_generated.ml | awk '{print $$2 ".ml"}') \
|
||||
Input_*.ml \
|
||||
qp_edit.ml qp_edit qp_edit.native Input_auto_generated.ml;\
|
||||
do \
|
||||
echo $$i ; \
|
||||
|
3
scripts/.gitignore
vendored
3
scripts/.gitignore
vendored
@ -2,3 +2,6 @@
|
||||
*.pyo
|
||||
docopt.py
|
||||
resultsFile/
|
||||
verif_omp/a.out
|
||||
src/*/Makefile
|
||||
src/*/*/
|
||||
|
11
src/.gitignore
vendored
Normal file
11
src/.gitignore
vendored
Normal file
@ -0,0 +1,11 @@
|
||||
*
|
||||
!README.rst
|
||||
!*/
|
||||
*/*
|
||||
!*/*.*
|
||||
*/*.o
|
||||
*/build.ninja
|
||||
*/ezfio_interface.irp.f
|
||||
*/.gitignore
|
||||
*/*.swp
|
||||
|
@ -38,11 +38,6 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
|
||||
|
||||
ao_integrals_n_e = 0.d0
|
||||
|
||||
! _
|
||||
! /| / |_)
|
||||
! | / | \
|
||||
!
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
|
||||
@ -106,7 +101,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
|
||||
endif
|
||||
|
||||
|
||||
IF(DO_PSEUDO) THEN
|
||||
IF(do_pseudo) THEN
|
||||
ao_integrals_n_e += ao_pseudo_integrals
|
||||
ENDIF
|
||||
|
||||
|
@ -292,9 +292,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
|
||||
! call omp_set_max_active_levels(1)
|
||||
|
||||
|
||||
print '(A)', '========== ======================= ===================== ===================== ==========='
|
||||
print '(A)', ' Samples Energy Variance Norm^2 Seconds'
|
||||
print '(A)', '========== ======================= ===================== ===================== ==========='
|
||||
print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
|
||||
print '(A)', ' Samples Energy PT2 Variance Norm^2 Convergence Seconds'
|
||||
print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
|
||||
|
||||
PROVIDE global_selection_buffer
|
||||
|
||||
@ -319,7 +319,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
|
||||
call set_multiple_levels_omp(.True.)
|
||||
! call omp_set_max_active_levels(8)
|
||||
|
||||
print '(A)', '========== ======================= ===================== ===================== ==========='
|
||||
print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
|
||||
|
||||
|
||||
do k=1,N_states
|
||||
pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate)
|
||||
@ -417,6 +418,17 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
|
||||
double precision :: rss
|
||||
double precision, external :: memory_of_double, memory_of_int
|
||||
|
||||
character(len=20) :: format_str1, str_error1, format_str2, str_error2
|
||||
character(len=20) :: format_str3, str_error3, format_str4, str_error4
|
||||
character(len=20) :: format_value1, format_value2, format_value3, format_value4
|
||||
character(len=20) :: str_value1, str_value2, str_value3, str_value4
|
||||
character(len=20) :: str_conv
|
||||
double precision :: value1, value2, value3, value4
|
||||
double precision :: error1, error2, error3, error4
|
||||
integer :: size1,size2,size3,size4
|
||||
|
||||
double precision :: conv_crit
|
||||
|
||||
sending =.False.
|
||||
|
||||
rss = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2)
|
||||
@ -540,14 +552,60 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
|
||||
|
||||
if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
|
||||
time1 = time
|
||||
print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1)', c, &
|
||||
pt2_data % pt2(pt2_stoch_istate) +E, &
|
||||
pt2_data_err % pt2(pt2_stoch_istate), &
|
||||
pt2_data % variance(pt2_stoch_istate), &
|
||||
pt2_data_err % variance(pt2_stoch_istate), &
|
||||
pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), &
|
||||
pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
|
||||
time-time0
|
||||
|
||||
value1 = pt2_data % pt2(pt2_stoch_istate) + E
|
||||
error1 = pt2_data_err % pt2(pt2_stoch_istate)
|
||||
value2 = pt2_data % pt2(pt2_stoch_istate)
|
||||
error2 = pt2_data_err % pt2(pt2_stoch_istate)
|
||||
value3 = pt2_data % variance(pt2_stoch_istate)
|
||||
error3 = pt2_data_err % variance(pt2_stoch_istate)
|
||||
value4 = pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate)
|
||||
error4 = pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate)
|
||||
|
||||
! Max size of the values (FX.Y) with X=size
|
||||
size1 = 15
|
||||
size2 = 12
|
||||
size3 = 12
|
||||
size4 = 12
|
||||
|
||||
! To generate the format: number(error)
|
||||
call format_w_error(value1,error1,size1,8,format_value1,str_error1)
|
||||
call format_w_error(value2,error2,size2,8,format_value2,str_error2)
|
||||
call format_w_error(value3,error3,size3,8,format_value3,str_error3)
|
||||
call format_w_error(value4,error4,size4,8,format_value4,str_error4)
|
||||
|
||||
! value > string with the right format
|
||||
write(str_value1,'('//format_value1//')') value1
|
||||
write(str_value2,'('//format_value2//')') value2
|
||||
write(str_value3,'('//format_value3//')') value3
|
||||
write(str_value4,'('//format_value4//')') value4
|
||||
|
||||
! Convergence criterion
|
||||
conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) / &
|
||||
(1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
|
||||
write(str_conv,'(G10.3)') conv_crit
|
||||
|
||||
write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,&
|
||||
adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),&
|
||||
adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),&
|
||||
adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),&
|
||||
adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),&
|
||||
adjustl(str_conv),&
|
||||
time-time0
|
||||
|
||||
! Old print
|
||||
!print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1,1pE16.6,1pE16.6)', c, &
|
||||
! pt2_data % pt2(pt2_stoch_istate) +E, &
|
||||
! pt2_data_err % pt2(pt2_stoch_istate), &
|
||||
! pt2_data % variance(pt2_stoch_istate), &
|
||||
! pt2_data_err % variance(pt2_stoch_istate), &
|
||||
! pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), &
|
||||
! pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
|
||||
! time-time0, &
|
||||
! pt2_data % pt2(pt2_stoch_istate), &
|
||||
! dabs(pt2_data_err % pt2(pt2_stoch_istate)) / &
|
||||
! (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
|
||||
|
||||
if (stop_now .or. ( &
|
||||
(do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) / &
|
||||
(1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then
|
||||
|
@ -5,7 +5,6 @@ subroutine run_slave_cipsi
|
||||
END_DOC
|
||||
|
||||
call set_multiple_levels_omp(.False.)
|
||||
! call omp_set_max_active_levels(1)
|
||||
distributed_davidson = .False.
|
||||
read_wf = .False.
|
||||
SOFT_TOUCH read_wf distributed_davidson
|
||||
@ -173,10 +172,8 @@ subroutine run_slave_main
|
||||
|
||||
!---
|
||||
call set_multiple_levels_omp(.True.)
|
||||
! call omp_set_max_active_levels(8)
|
||||
call davidson_slave_tcp(0)
|
||||
call set_multiple_levels_omp(.False.)
|
||||
! call omp_set_max_active_levels(1)
|
||||
print *, mpi_rank, ': Davidson done'
|
||||
!---
|
||||
|
||||
|
@ -74,7 +74,7 @@ subroutine run
|
||||
print*,'******************************************************'
|
||||
print*,'Excitation energies (au) (eV)'
|
||||
do i = 2, N_states
|
||||
print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0
|
||||
print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1)) * ha_to_ev
|
||||
enddo
|
||||
print*,''
|
||||
endif
|
||||
|
@ -856,6 +856,7 @@ end subroutine
|
||||
!end subroutine
|
||||
!
|
||||
BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]
|
||||
|
||||
implicit none
|
||||
@ -944,6 +945,29 @@ end subroutine
|
||||
enddo
|
||||
|
||||
deallocate(dets, old_order)
|
||||
integer :: ndet_conf
|
||||
do i = 1, N_configuration
|
||||
ndet_conf = psi_configuration_to_psi_det(2,i) - psi_configuration_to_psi_det(1,i) + 1
|
||||
psi_configuration_n_det(i) = ndet_conf
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, n_elec_alpha_for_psi_configuration, (N_configuration)]
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
integer(bit_kind) :: det_tmp(N_int,2),det_alpha(N_int)
|
||||
n_elec_alpha_for_psi_configuration = 0
|
||||
do i = 1, N_configuration
|
||||
j = psi_configuration_to_psi_det(2,i)
|
||||
det_tmp(:,:) = psi_det(:,:,j)
|
||||
k = 0
|
||||
do l = 1, N_int
|
||||
det_alpha(N_int) = iand(det_tmp(l,1),psi_configuration(l,1,i))
|
||||
k += popcnt(det_alpha(l))
|
||||
enddo
|
||||
n_elec_alpha_for_psi_configuration(i) = k
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
@ -10,7 +10,6 @@ BEGIN_PROVIDER [ double precision, psi_csf_coef, (N_csf, N_states) ]
|
||||
call convertWFfromDETtoCSF(N_states, buffer, psi_csf_coef)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
|
||||
use cfunctions
|
||||
use bitmasks
|
||||
|
@ -509,7 +509,6 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
|
||||
|
||||
|
||||
call set_multiple_levels_omp(.True.)
|
||||
! call omp_set_max_active_levels(5)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
|
||||
ithread = omp_get_thread_num()
|
||||
|
@ -465,7 +465,7 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze)
|
||||
endif
|
||||
|
||||
call set_multiple_levels_omp(.True.)
|
||||
! call omp_set_max_active_levels(4)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
|
||||
ithread = omp_get_thread_num()
|
||||
if (ithread == 0 ) then
|
||||
|
@ -465,7 +465,7 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze)
|
||||
endif
|
||||
|
||||
call set_multiple_levels_omp(.True.)
|
||||
! call omp_set_max_active_levels(4)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
|
||||
ithread = omp_get_thread_num()
|
||||
if (ithread == 0 ) then
|
||||
|
@ -300,7 +300,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
|
||||
shift = N_st_diag*(iter-1)
|
||||
shift2 = N_st_diag*iter
|
||||
|
||||
if ((iter > 1).or.(itertot == 1)) then
|
||||
! if ((iter > 1).or.(itertot == 1)) then
|
||||
! Compute |W_k> = \sum_i |i><i|H|u_k>
|
||||
! -----------------------------------
|
||||
|
||||
@ -310,10 +310,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
|
||||
else
|
||||
call H_u_0_nstates_openmp(W,U,N_st_diag,sze)
|
||||
endif
|
||||
else
|
||||
! Already computed in update below
|
||||
continue
|
||||
endif
|
||||
! else
|
||||
! ! Already computed in update below
|
||||
! continue
|
||||
! endif
|
||||
|
||||
if (dressing_state > 0) then
|
||||
|
||||
|
@ -77,9 +77,9 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
|
||||
!$OMP END PARALLEL
|
||||
|
||||
if (dressing_state > 0) then
|
||||
do k = 1, N_st
|
||||
do k=1,N_st
|
||||
|
||||
do i = 1, sze
|
||||
do i=1,sze
|
||||
H_jj(i) += u_in(i,k) * dressing_column_h(i,k)
|
||||
enddo
|
||||
|
||||
|
@ -304,7 +304,7 @@ BEGIN_PROVIDER [ double precision, c0_weight, (N_states) ]
|
||||
c0_weight(i) = c0_weight(i) * c
|
||||
enddo
|
||||
else
|
||||
c0_weight = 1.d0
|
||||
c0_weight(:) = 1.d0
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
@ -321,7 +321,7 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
|
||||
if (weight_one_e_dm == 0) then
|
||||
state_average_weight(:) = c0_weight(:)
|
||||
else if (weight_one_e_dm == 1) then
|
||||
state_average_weight(:) = 1./N_states
|
||||
state_average_weight(:) = 1.d0/N_states
|
||||
else
|
||||
call ezfio_has_determinants_state_average_weight(exists)
|
||||
if (exists) then
|
||||
@ -384,6 +384,14 @@ END_PROVIDER
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_e_dm_ao, (ao_num, ao_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta
|
||||
END_DOC
|
||||
one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine get_occupation_from_dets(istate,occupation)
|
||||
implicit none
|
||||
|
@ -652,7 +652,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
|
||||
case (1)
|
||||
call get_single_excitation(key_i,key_j,exc,phase,Nint)
|
||||
!DIR$ FORCEINLINE
|
||||
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
|
||||
if (exc(0,1,1) == 1) then
|
||||
! Single alpha
|
||||
m = exc(1,1,1)
|
||||
@ -671,10 +670,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
|
||||
end select
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble,phase)
|
||||
use bitmasks
|
||||
implicit none
|
||||
@ -1009,7 +1004,6 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
@ -282,9 +282,7 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij)
|
||||
double precision :: get_two_e_integral
|
||||
integer :: m,n,p,q
|
||||
integer :: i,j,k
|
||||
integer :: occ(Nint*bit_kind_size,2)
|
||||
double precision :: diag_H_mat_elem, phase,phase_2
|
||||
integer :: n_occ_ab(2)
|
||||
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals ref_bitmask_two_e_energy
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
@ -342,7 +340,6 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij)
|
||||
case (1)
|
||||
call get_single_excitation(key_i,key_j,exc,phase,Nint)
|
||||
!DIR$ FORCEINLINE
|
||||
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
|
||||
if (exc(0,1,1) == 1) then
|
||||
! Mono alpha
|
||||
m = exc(1,1,1)
|
||||
|
@ -73,7 +73,6 @@ subroutine run_dress_slave(thread,iproce,energy)
|
||||
ending = dress_N_cp+1
|
||||
ntask_tbd = 0
|
||||
call set_multiple_levels_omp(.True.)
|
||||
! call omp_set_max_active_levels(8)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(SHARED) &
|
||||
!$OMP PRIVATE(interesting, breve_delta_m, task_id) &
|
||||
@ -86,7 +85,6 @@ subroutine run_dress_slave(thread,iproce,energy)
|
||||
integer, external :: connect_to_taskserver
|
||||
!$OMP CRITICAL
|
||||
call set_multiple_levels_omp(.False.)
|
||||
! call omp_set_max_active_levels(1)
|
||||
if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
|
||||
print *, irp_here, ': Unable to connect to task server'
|
||||
stop -1
|
||||
@ -299,7 +297,6 @@ subroutine run_dress_slave(thread,iproce,energy)
|
||||
|
||||
!$OMP END PARALLEL
|
||||
call set_multiple_levels_omp(.False.)
|
||||
! call omp_set_max_active_levels(1)
|
||||
! do i=0,dress_N_cp+1
|
||||
! call omp_destroy_lock(lck_sto(i))
|
||||
! end do
|
||||
|
51
src/ezfio_files/trexio.irp.f
Normal file
51
src/ezfio_files/trexio.irp.f
Normal file
@ -0,0 +1,51 @@
|
||||
use trexio
|
||||
|
||||
BEGIN_PROVIDER [ logical, use_trexio ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the content of the QP_USE_TREXIO variable
|
||||
END_DOC
|
||||
character(32) :: buffer
|
||||
|
||||
call getenv('QP_USE_TREXIO', buffer)
|
||||
if (trim(buffer) == '0') then
|
||||
print *, 'Using EZFIO storage'
|
||||
use_trexio = .False.
|
||||
else
|
||||
print *, 'Using TREXIO storage'
|
||||
use_trexio = .True.
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ character*(1024), trexio_filename ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Name of the TREXIO file.
|
||||
END_DOC
|
||||
trexio_filename = trim(ezfio_work_dir)//'/trexio.h5'
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(trexio_back_end_t), trexio_backend ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Name of the TREXIO file.
|
||||
END_DOC
|
||||
trexio_backend = TREXIO_HDF5
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(trexio_t), trexio_file ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Name of the TREXIO file.
|
||||
END_DOC
|
||||
integer (trexio_exit_code) :: rc
|
||||
|
||||
trexio_file = 0_trexio_t
|
||||
if (use_trexio) then
|
||||
trexio_file = trexio_open(trexio_filename, 'u', trexio_backend, rc)
|
||||
call trexio_assert(rc, TREXIO_SUCCESS)
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
@ -35,12 +35,13 @@ subroutine print_extrapolated_energy
|
||||
do k=2,min(N_iter,8)
|
||||
write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
|
||||
extrapolated_energy(k,i) - extrapolated_energy(k,1), &
|
||||
(extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
|
||||
(extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * ha_to_ev
|
||||
enddo
|
||||
write(*,*) '=========== ', '=================== ', '=================== ', '==================='
|
||||
enddo
|
||||
|
||||
print *, ''
|
||||
call ezfio_set_fci_energy_extrapolated(extrapolated_energy(min(N_iter,3),1:N_states))
|
||||
|
||||
end subroutine
|
||||
|
||||
|
@ -36,7 +36,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
|
||||
write(*,fmt) '# E ', e_(1:N_states_p)
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
|
||||
write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
|
||||
write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*ha_to_ev
|
||||
endif
|
||||
write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
|
||||
write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p)
|
||||
@ -47,8 +47,8 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), &
|
||||
dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p)
|
||||
write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, &
|
||||
dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p)
|
||||
write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*ha_to_ev, &
|
||||
dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*ha_to_ev, k=1,N_states_p)
|
||||
endif
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
@ -82,19 +82,19 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
|
||||
print *, 'Variational Energy difference (au | eV)'
|
||||
do i=2, N_states_p
|
||||
print*,'Delta E = ', (e_(i) - e_(1)), &
|
||||
(e_(i) - e_(1)) * 27.211396641308d0
|
||||
(e_(i) - e_(1)) * ha_to_ev
|
||||
enddo
|
||||
print *, '-----'
|
||||
print*, 'Variational + perturbative Energy difference (au | eV)'
|
||||
do i=2, N_states_p
|
||||
print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), &
|
||||
(e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0
|
||||
(e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * ha_to_ev
|
||||
enddo
|
||||
print *, '-----'
|
||||
print*, 'Variational + renormalized perturbative Energy difference (au | eV)'
|
||||
do i=2, N_states_p
|
||||
print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), &
|
||||
(e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0
|
||||
(e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * ha_to_ev
|
||||
enddo
|
||||
endif
|
||||
|
||||
|
@ -53,7 +53,11 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
|
||||
! call four_idx_novvvv
|
||||
call four_idx_novvvv_old
|
||||
else
|
||||
call add_integrals_to_map(full_ijkl_bitmask_4)
|
||||
if (32.d-9*dble(ao_num)**4 < dble(qp_max_mem)) then
|
||||
call four_idx_dgemm
|
||||
else
|
||||
call add_integrals_to_map(full_ijkl_bitmask_4)
|
||||
endif
|
||||
endif
|
||||
|
||||
call wall_time(wall_2)
|
||||
@ -77,6 +81,94 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
subroutine four_idx_dgemm
|
||||
implicit none
|
||||
integer :: p,q,r,s,i,j,k,l
|
||||
double precision, allocatable :: a1(:,:,:,:)
|
||||
double precision, allocatable :: a2(:,:,:,:)
|
||||
|
||||
allocate (a1(ao_num,ao_num,ao_num,ao_num))
|
||||
|
||||
print *, 'Getting AOs'
|
||||
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,r,s)
|
||||
do s=1,ao_num
|
||||
do r=1,ao_num
|
||||
do q=1,ao_num
|
||||
call get_ao_two_e_integrals(q,r,s,ao_num,a1(1,q,r,s))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
print *, '1st transformation'
|
||||
! 1st transformation
|
||||
allocate (a2(ao_num,ao_num,ao_num,mo_num))
|
||||
call dgemm('T','N', (ao_num*ao_num*ao_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*ao_num*ao_num))
|
||||
|
||||
! 2nd transformation
|
||||
print *, '2nd transformation'
|
||||
deallocate (a1)
|
||||
allocate (a1(ao_num,ao_num,mo_num,mo_num))
|
||||
call dgemm('T','N', (ao_num*ao_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (ao_num*ao_num*mo_num))
|
||||
|
||||
! 3rd transformation
|
||||
print *, '3rd transformation'
|
||||
deallocate (a2)
|
||||
allocate (a2(ao_num,mo_num,mo_num,mo_num))
|
||||
call dgemm('T','N', (ao_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*mo_num*mo_num))
|
||||
|
||||
! 4th transformation
|
||||
print *, '4th transformation'
|
||||
deallocate (a1)
|
||||
allocate (a1(mo_num,mo_num,mo_num,mo_num))
|
||||
call dgemm('T','N', (mo_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (mo_num*mo_num*mo_num))
|
||||
|
||||
deallocate (a2)
|
||||
|
||||
integer :: n_integrals, size_buffer
|
||||
integer(key_kind) , allocatable :: buffer_i(:)
|
||||
real(integral_kind), allocatable :: buffer_value(:)
|
||||
size_buffer = min(ao_num*ao_num*ao_num,16000000)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals)
|
||||
allocate ( buffer_i(size_buffer), buffer_value(size_buffer) )
|
||||
|
||||
n_integrals = 0
|
||||
!$OMP DO
|
||||
do l=1,mo_num
|
||||
do k=1,mo_num
|
||||
do j=1,l
|
||||
do i=1,k
|
||||
if (abs(a1(i,j,k,l)) < mo_integrals_threshold) then
|
||||
cycle
|
||||
endif
|
||||
n_integrals += 1
|
||||
buffer_value(n_integrals) = a1(i,j,k,l)
|
||||
!DIR$ FORCEINLINE
|
||||
call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals))
|
||||
if (n_integrals == size_buffer) then
|
||||
call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
|
||||
n_integrals = 0
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
|
||||
|
||||
deallocate(buffer_i, buffer_value)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
deallocate (a1)
|
||||
|
||||
call map_unique(mo_integrals_map)
|
||||
|
||||
integer*8 :: get_mo_map_size, mo_map_size
|
||||
mo_map_size = get_mo_map_size()
|
||||
|
||||
end subroutine
|
||||
|
||||
subroutine add_integrals_to_map(mask_ijkl)
|
||||
use bitmasks
|
||||
|
@ -73,11 +73,6 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
|
||||
liwork = -1
|
||||
|
||||
F_save = F
|
||||
!print *, ' Fock matrix'
|
||||
!do i = 1, mo_num
|
||||
! write(*, '(1000(F16.10,X))') F_save(:,i)
|
||||
!enddo
|
||||
|
||||
call dsyevd( 'V', 'U', mo_num, F, &
|
||||
size(F,1), diag, work, lwork, iwork, liwork, info)
|
||||
|
||||
@ -108,16 +103,6 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
|
||||
endif
|
||||
endif
|
||||
|
||||
!print *, ' eigenvalues'
|
||||
!do i = 1, mo_num
|
||||
! write(*, '(1000(F16.10,X))') diag(i)
|
||||
!enddo
|
||||
!print *, ' eigenvectors'
|
||||
!do i = 1, mo_num
|
||||
! write(*, '(1000(F16.10,X))') F(:,i)
|
||||
!enddo
|
||||
|
||||
|
||||
call dgemm('N','N',ao_num,mo_num,mo_num, 1.d0, &
|
||||
mo_coef, size(mo_coef,1), F, size(F,1), &
|
||||
0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
|
||||
|
@ -12,16 +12,6 @@ BEGIN_PROVIDER [double precision, SCF_density_matrix_ao_alpha, (ao_num,ao_num) ]
|
||||
SCF_density_matrix_ao_alpha = 0.d0
|
||||
endif
|
||||
|
||||
! integer :: i, j
|
||||
! double precision :: trace_density
|
||||
! trace_density = 0.d0
|
||||
! do i = 1, ao_num !elec_alpha_num
|
||||
! do j = 1, ao_num !elec_alpha_num
|
||||
! trace_density = trace_density &
|
||||
! + SCF_density_matrix_ao_alpha(j,i) * ao_overlap(j,i)
|
||||
! enddo
|
||||
! enddo
|
||||
! print *, ' trace of SCF_density_matrix_ao_alpha =', trace_density
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -52,8 +52,8 @@ program molden
|
||||
l += 1
|
||||
if (l > ao_num) exit
|
||||
enddo
|
||||
write(i_unit_output,*)''
|
||||
enddo
|
||||
write(i_unit_output,*)''
|
||||
enddo
|
||||
|
||||
|
||||
|
27
src/tools/print_cfg.irp.f
Normal file
27
src/tools/print_cfg.irp.f
Normal file
@ -0,0 +1,27 @@
|
||||
program print_energy
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Prints the energy of the wave function stored in the |EZFIO| directory.
|
||||
END_DOC
|
||||
|
||||
! this has to be done in order to be sure that N_det, psi_det and
|
||||
! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
|
||||
read_wf = .True.
|
||||
touch read_wf
|
||||
PROVIDE N_states
|
||||
call run
|
||||
end
|
||||
|
||||
subroutine run
|
||||
implicit none
|
||||
integer :: i,j
|
||||
|
||||
do i=1,N_configuration
|
||||
print *, i, sum(popcnt(psi_configuration(:,1,i)))
|
||||
enddo
|
||||
|
||||
print *, ''
|
||||
do i=0,elec_num
|
||||
print *, i, cfg_seniority_index(i)
|
||||
enddo
|
||||
end
|
@ -33,8 +33,9 @@ subroutine routine
|
||||
double precision :: norm_mono_a,norm_mono_b
|
||||
double precision :: norm_mono_a_2,norm_mono_b_2
|
||||
double precision :: norm_mono_a_pert_2,norm_mono_b_pert_2
|
||||
double precision :: norm_mono_a_pert,norm_mono_b_pert
|
||||
double precision :: norm_mono_a_pert,norm_mono_b_pert,norm_double_1
|
||||
double precision :: delta_e,coef_2_2
|
||||
|
||||
norm_mono_a = 0.d0
|
||||
norm_mono_b = 0.d0
|
||||
norm_mono_a_2 = 0.d0
|
||||
@ -43,6 +44,7 @@ subroutine routine
|
||||
norm_mono_b_pert = 0.d0
|
||||
norm_mono_a_pert_2 = 0.d0
|
||||
norm_mono_b_pert_2 = 0.d0
|
||||
norm_double_1 = 0.d0
|
||||
do i = 1, min(N_det_print_wf,N_det)
|
||||
print*,''
|
||||
print*,'i = ',i
|
||||
@ -94,6 +96,7 @@ subroutine routine
|
||||
print*,'h1,p1 = ',h1,p1
|
||||
print*,'s2',s2
|
||||
print*,'h2,p2 = ',h2,p2
|
||||
norm_double_1 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1))
|
||||
endif
|
||||
|
||||
print*,'<Ref| H |D_I> = ',hij
|
||||
@ -110,6 +113,7 @@ subroutine routine
|
||||
print*,''
|
||||
print*,'L1 norm of mono alpha = ',norm_mono_a
|
||||
print*,'L1 norm of mono beta = ',norm_mono_b
|
||||
print*,'L1 norm of double exc = ',norm_double_1
|
||||
print*, '---'
|
||||
print*,'L2 norm of mono alpha = ',norm_mono_a_2
|
||||
print*,'L2 norm of mono beta = ',norm_mono_b_2
|
||||
|
@ -55,10 +55,23 @@ subroutine routine_s2
|
||||
integer :: i,j,k
|
||||
double precision :: accu(N_states)
|
||||
|
||||
print *, 'Weights of the CFG'
|
||||
integer :: weights(0:16), ix
|
||||
double precision :: x
|
||||
|
||||
weights(:) = 0
|
||||
do i=1,N_det
|
||||
print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
|
||||
x = -dlog(1.d-32+sum(weight_configuration(det_to_configuration(i),:)))/dlog(10.d0)
|
||||
ix = min(int(x), 16)
|
||||
weights(ix) += 1
|
||||
enddo
|
||||
|
||||
print *, 'Histogram of the weights of the CFG'
|
||||
do i=0,15
|
||||
print *, ' 10^{-', i, '} ', weights(i)
|
||||
end do
|
||||
print *, '< 10^{-', 15, '} ', weights(16)
|
||||
|
||||
|
||||
print*, 'Min weight of the configuration?'
|
||||
read(5,*) wmin
|
||||
|
||||
|
@ -529,10 +529,14 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_
|
||||
c_average += c_1(l) * c_1(l) * state_weights(l)
|
||||
enddo
|
||||
|
||||
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
|
||||
if (nkeys > 0) then
|
||||
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
|
||||
endif
|
||||
nkeys = 0
|
||||
call orb_range_diag_to_all_2_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
|
||||
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
|
||||
if (nkeys > 0) then
|
||||
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
|
||||
endif
|
||||
nkeys = 0
|
||||
|
||||
end do
|
||||
|
71
src/utils/format_w_error.irp.f
Normal file
71
src/utils/format_w_error.irp.f
Normal file
@ -0,0 +1,71 @@
|
||||
subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error)
|
||||
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Format for double precision, value(error)
|
||||
END_DOC
|
||||
|
||||
! in
|
||||
! | value | double precision | value... |
|
||||
! | error | double precision | error... |
|
||||
! | size_nb | integer | X in FX.Y |
|
||||
! | max_nb_digits | integer | Max Y in FX.Y |
|
||||
|
||||
! out
|
||||
! | format_value | character | string FX.Y for the format |
|
||||
! | str_error | character | string of the error |
|
||||
|
||||
! internal
|
||||
! | str_size | character | size in string format |
|
||||
! | nb_digits | integer | number of digits Y in FX.Y depending of the error |
|
||||
! | str_nb_digits | character | nb_digits in string format |
|
||||
! | str_exp | character | string of the value in exponential format |
|
||||
|
||||
! in
|
||||
double precision, intent(in) :: error, value
|
||||
integer, intent(in) :: size_nb, max_nb_digits
|
||||
|
||||
! out
|
||||
character(len=20), intent(out) :: str_error, format_value
|
||||
|
||||
! internal
|
||||
character(len=20) :: str_size, str_nb_digits, str_exp
|
||||
integer :: nb_digits
|
||||
|
||||
! max_nb_digit: Y max
|
||||
! size_nb = Size of the double: X (FX.Y)
|
||||
write(str_size,'(I3)') size_nb
|
||||
|
||||
! Error
|
||||
write(str_exp,'(1pE20.0)') error
|
||||
str_error = trim(adjustl(str_exp))
|
||||
|
||||
! Number of digit: Y (FX.Y) from the exponent
|
||||
str_nb_digits = str_exp(19:20)
|
||||
read(str_nb_digits,*) nb_digits
|
||||
|
||||
! If the error is 0d0
|
||||
if (error <= 1d-16) then
|
||||
write(str_nb_digits,*) max_nb_digits
|
||||
endif
|
||||
|
||||
! If the error is too small
|
||||
if (nb_digits > max_nb_digits) then
|
||||
write(str_nb_digits,*) max_nb_digits
|
||||
str_error(1:1) = '0'
|
||||
endif
|
||||
|
||||
! If the error is too big (>= 0.5)
|
||||
if (error >= 0.5d0) then
|
||||
str_nb_digits = '1'
|
||||
str_error(1:1) = '*'
|
||||
endif
|
||||
|
||||
! FX.Y,A1,A1,A1 for value(str_error)
|
||||
!string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1'
|
||||
|
||||
! FX.Y just for the value
|
||||
format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))
|
||||
|
||||
end
|
22
src/utils/units.irp.f
Normal file
22
src/utils/units.irp.f
Normal file
@ -0,0 +1,22 @@
|
||||
BEGIN_PROVIDER [double precision, ha_to_ev]
|
||||
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Converstion from Hartree to eV
|
||||
END_DOC
|
||||
|
||||
ha_to_ev = 27.211396641308d0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, au_to_D]
|
||||
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Converstion from au to Debye
|
||||
END_DOC
|
||||
|
||||
au_to_D = 2.5415802529d0
|
||||
|
||||
END_PROVIDER
|
||||
|
@ -37,6 +37,10 @@ double precision function binom_func(i,j)
|
||||
else
|
||||
binom_func = dexp( logfact(i)-logfact(j)-logfact(i-j) )
|
||||
endif
|
||||
|
||||
! To avoid .999999 numbers
|
||||
binom_func = floor(binom_func + 0.5d0)
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user