diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index fd1739f1..be2bf18c 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -150,12 +150,12 @@ end = struct let to_string = function | Huckel -> \"Huckel\" - | HCore -> \"Hcore\" + | HCore -> \"HCore\" let of_string s = - match s with - | \"Huckel\" -> Huckel - | \"Hcore\" -> HCore + match (String.lowercase s) with + | \"huckel\" -> Huckel + | \"hcore\" -> HCore | _ -> failwith (\"Wrong Guess type : \"^s) end @@ -176,10 +176,10 @@ end = struct | Write -> \"Write\" | None -> \"None\" let of_string s = - match s with - | \"Read\" -> Read - | \"Write\" -> Write - | \"None\" -> None + match (String.lowercase s) with + | \"read\" -> Read + | \"write\" -> Write + | \"none\" -> None | _ -> failwith (\"Wrong IO type : \"^s) end diff --git a/plugins/CAS_SD/tree_dependency.png b/plugins/CAS_SD/tree_dependency.png index d3d98e3c..185c2b27 100644 Binary files a/plugins/CAS_SD/tree_dependency.png and b/plugins/CAS_SD/tree_dependency.png differ diff --git a/plugins/Full_CI/.gitignore b/plugins/Full_CI/.gitignore index a806bcbc..2d978fa6 100644 --- a/plugins/Full_CI/.gitignore +++ b/plugins/Full_CI/.gitignore @@ -22,6 +22,7 @@ Properties Pseudo Selectors_full Utils +exc_degree ezfio_interface.irp.f full_ci full_ci_no_skip diff --git a/plugins/Full_CI/tree_dependency.png b/plugins/Full_CI/tree_dependency.png index ad7dcd84..caedb2e0 100644 Binary files a/plugins/Full_CI/tree_dependency.png and b/plugins/Full_CI/tree_dependency.png differ diff --git a/plugins/Generators_CAS/tree_dependency.png b/plugins/Generators_CAS/tree_dependency.png index 7f32349e..5bbc55d0 100644 Binary files a/plugins/Generators_CAS/tree_dependency.png and b/plugins/Generators_CAS/tree_dependency.png differ diff --git a/plugins/Generators_full/tree_dependency.png b/plugins/Generators_full/tree_dependency.png index c2656d20..94ad6358 100644 Binary files a/plugins/Generators_full/tree_dependency.png and b/plugins/Generators_full/tree_dependency.png differ diff --git a/plugins/Hartree_Fock/tree_dependency.png b/plugins/Hartree_Fock/tree_dependency.png index 8a43fb0d..cb1d9738 100644 Binary files a/plugins/Hartree_Fock/tree_dependency.png and b/plugins/Hartree_Fock/tree_dependency.png differ diff --git a/plugins/MRCC/H_apply.irp.f b/plugins/MRCC/H_apply.irp.f deleted file mode 100644 index 7e479529..00000000 --- a/plugins/MRCC/H_apply.irp.f +++ /dev/null @@ -1,28 +0,0 @@ -use bitmasks -BEGIN_SHELL [ /usr/bin/env python ] -from generate_h_apply import * - -s = H_apply("mrcc") -s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_cas, Ndet_non_cas" -s.data["declarations"] += """ - integer, intent(in) :: Ndet_cas,Ndet_non_cas - double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*) - double precision, intent(in) :: delta_ii_(Ndet_cas,*) -""" -s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_cas,Ndet_non_cas,i_generator,key_idx,keys_out,N_int,iproc)" -s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas" -s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas" -s.data["decls_main"] += """ - integer, intent(in) :: Ndet_cas,Ndet_non_cas - double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*) - double precision, intent(in) :: delta_ii_(Ndet_cas,*) -""" -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/NEEDED_CHILDREN_MODULES b/plugins/MRCC/NEEDED_CHILDREN_MODULES deleted file mode 100644 index 04ce9e78..00000000 --- a/plugins/MRCC/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -Perturbation Selectors_full Generators_full diff --git a/plugins/MRCC/mrcc.irp.f b/plugins/MRCC/mrcc.irp.f deleted file mode 100644 index 8a35dde5..00000000 --- a/plugins/MRCC/mrcc.irp.f +++ /dev/null @@ -1,70 +0,0 @@ -program mrcc - implicit none - read_wf = .True. - TOUCH read_wf - call run - call run_mrcc -! call run_mrcc_test -end - -subroutine run - implicit none - - integer :: i,j - print *, 'CAS' - print *, '===' - do i=1,N_det_cas - print *, psi_cas_coef(i,:) - call debug_det(psi_cas(1,1,i),N_int) - enddo - -! print *, 'SD' -! print *, '==' -! do i=1,N_det_non_cas -! print *, psi_non_cas_coef(i,:) -! call debug_det(psi_non_cas(1,1,i),N_int) -! enddo - call write_double(6,ci_energy(1),"Initial CI energy") -end -subroutine run_mrcc_test - implicit none - integer :: i,j - double precision :: pt2 - pt2 = 0.d0 - do j=1,N_det - do i=1,N_det - pt2 += psi_coef(i,1)*psi_coef(j,1) * delta_ij(i,j,1) - enddo - enddo - print *, ci_energy(1) - print *, ci_energy(1)+pt2 -end -subroutine run_mrcc - implicit none - integer :: i,j - - double precision :: E_new, E_old, delta_e - integer :: iteration - E_new = 0.d0 - delta_E = 1.d0 - iteration = 0 - do while (delta_E > 1.d-10) - iteration += 1 - print *, '===========================' - print *, 'MRCC Iteration', iteration - print *, '===========================' - print *, '' - E_old = sum(ci_energy_dressed) - call write_double(6,ci_energy_dressed(1),"MRCC energy") - call diagonalize_ci_dressed - E_new = sum(ci_energy_dressed) - delta_E = dabs(E_new - E_old) - if (iteration > 20) then - exit - endif - enddo - call write_double(6,ci_energy_dressed(1),"Final MRCC energy") - call ezfio_set_mrcc_energy(ci_energy_dressed(1)) -! call save_wavefunction - -end diff --git a/plugins/MRCC/tree_dependency.png b/plugins/MRCC/tree_dependency.png deleted file mode 100644 index db540bc5..00000000 Binary files a/plugins/MRCC/tree_dependency.png and /dev/null differ diff --git a/plugins/MRCC_CASSD/.gitignore b/plugins/MRCC_CASSD/.gitignore new file mode 100644 index 00000000..d81ca7b8 --- /dev/null +++ b/plugins/MRCC_CASSD/.gitignore @@ -0,0 +1,32 @@ +# Automatically created by $QP_ROOT/scripts/module/module_handler.py +.ninja_deps +.ninja_log +AO_Basis +Bitmask +Determinants +Electrons +Ezfio_files +Generators_full +Hartree_Fock +IRPF90_man +IRPF90_temp +Integrals_Bielec +Integrals_Monoelec +MOGuess +MO_Basis +MRCC_Utils +Makefile +Makefile.depend +Nuclei +Perturbation +Properties +Pseudo +Psiref_CAS +Psiref_Utils +Selectors_full +Utils +ezfio_interface.irp.f +irpf90.make +irpf90_entities +mrcc_cassd +tags \ No newline at end of file diff --git a/plugins/MRCC_CASSD/EZFIO.cfg b/plugins/MRCC_CASSD/EZFIO.cfg new file mode 100644 index 00000000..21cc5b98 --- /dev/null +++ b/plugins/MRCC_CASSD/EZFIO.cfg @@ -0,0 +1,4 @@ +[energy] +type: double precision +doc: Calculated energy +interface: ezfio diff --git a/plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES b/plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..a8404d62 --- /dev/null +++ b/plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils diff --git a/plugins/MRCC_CASSD/README.rst b/plugins/MRCC_CASSD/README.rst new file mode 100644 index 00000000..5ef5db62 --- /dev/null +++ b/plugins/MRCC_CASSD/README.rst @@ -0,0 +1,33 @@ +=========== +MRCC Module +=========== + +MRCC as a coupled cluster on a CAS+SD wave function. + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. by the `update_README.py` script. + +.. image:: tree_dependency.png + +* `Perturbation `_ +* `Selectors_full `_ +* `Generators_full `_ +* `Psiref_CAS `_ +* `MRCC_Utils `_ + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. by the `update_README.py` script. + +`mrcc `_ + Undocumented + + +`print_cas_coefs `_ + Undocumented + diff --git a/plugins/MRCC_CASSD/mrcc_cassd.irp.f b/plugins/MRCC_CASSD/mrcc_cassd.irp.f new file mode 100644 index 00000000..e784a167 --- /dev/null +++ b/plugins/MRCC_CASSD/mrcc_cassd.irp.f @@ -0,0 +1,24 @@ +program mrcc + implicit none + if (.not.read_wf) then + print *, 'read_wf has to be true.' + stop 1 + endif + call print_cas_coefs + call run_mrcc +end + +subroutine print_cas_coefs + implicit none + + integer :: i,j + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, psi_cas_coef(i,:) + call debug_det(psi_cas(1,1,i),N_int) + enddo + + call write_double(6,ci_energy(1),"Initial CI energy") +end + diff --git a/plugins/MRCC_Utils/.gitignore b/plugins/MRCC_Utils/.gitignore new file mode 100644 index 00000000..e6279f11 --- /dev/null +++ b/plugins/MRCC_Utils/.gitignore @@ -0,0 +1,31 @@ +# Automatically created by $QP_ROOT/scripts/module/module_handler.py +.ninja_deps +.ninja_log +AO_Basis +Bitmask +Determinants +Electrons +Ezfio_files +Generators_full +Hartree_Fock +IRPF90_man +IRPF90_temp +Integrals_Bielec +Integrals_Monoelec +MOGuess +MO_Basis +Makefile +Makefile.depend +Nuclei +Perturbation +Properties +Pseudo +Psiref_CAS +Psiref_Utils +Selectors_full +Utils +ezfio_interface.irp.f +irpf90.make +irpf90_entities +mrcc_general +tags \ No newline at end of file diff --git a/plugins/MRCC/EZFIO.cfg b/plugins/MRCC_Utils/EZFIO.cfg similarity index 100% rename from plugins/MRCC/EZFIO.cfg rename to plugins/MRCC_Utils/EZFIO.cfg diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f new file mode 100644 index 00000000..455c37da --- /dev/null +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -0,0 +1,28 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("mrcc") +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 mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc)" +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/NEEDED_CHILDREN_MODULES b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..5b16423e --- /dev/null +++ b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Perturbation Selectors_full Generators_full Psiref_Utils diff --git a/plugins/MRCC/README.rst b/plugins/MRCC_Utils/README.rst similarity index 75% rename from plugins/MRCC/README.rst rename to plugins/MRCC_Utils/README.rst index b1a74ccc..71392798 100644 --- a/plugins/MRCC/README.rst +++ b/plugins/MRCC_Utils/README.rst @@ -13,6 +13,7 @@ Needed Modules * `Perturbation `_ * `Selectors_full `_ * `Generators_full `_ +* `Psiref_Utils `_ Documentation ============= @@ -20,23 +21,23 @@ Documentation .. Do not edit this section. It was auto-generated from the .. by the `update_README.py` script. -`ci_eigenvectors_dressed `_ +`ci_eigenvectors_dressed `_ Eigenvectors/values of the CI matrix -`ci_eigenvectors_s2_dressed `_ +`ci_eigenvectors_s2_dressed `_ Eigenvectors/values of the CI matrix -`ci_electronic_energy_dressed `_ +`ci_electronic_energy_dressed `_ Eigenvectors/values of the CI matrix -`ci_energy_dressed `_ +`ci_energy_dressed `_ N_states lowest eigenvalues of the dressed CI matrix -`davidson_diag_hjj_mrcc `_ +`davidson_diag_hjj_mrcc `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -57,7 +58,7 @@ Documentation Initial guess vectors are not necessarily orthonormal -`davidson_diag_mrcc `_ +`davidson_diag_mrcc `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -76,45 +77,45 @@ Documentation Initial guess vectors are not necessarily orthonormal -`delta_ii `_ +`delta_ii `_ Dressing matrix in N_det basis -`delta_ij `_ +`delta_ij `_ Dressing matrix in N_det basis -`diagonalize_ci_dressed `_ +`diagonalize_ci_dressed `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`find_triples_and_quadruples `_ +`find_triples_and_quadruples `_ Undocumented -`h_apply_mrcc `_ +`h_apply_mrcc `_ Calls H_apply on the HF determinant and selects all connected single and double excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. -`h_apply_mrcc_diexc `_ +`h_apply_mrcc_diexc `_ Generate all double excitations of key_in using the bit masks of holes and particles. Assume N_int is already provided. -`h_apply_mrcc_monoexc `_ +`h_apply_mrcc_monoexc `_ Generate all single excitations of key_in using the bit masks of holes and particles. Assume N_int is already provided. -`h_matrix_dressed `_ +`h_matrix_dressed `_ Dressed H with Delta_ij -`h_u_0_mrcc `_ +`h_u_0_mrcc `_ Computes v_0 = H|u_0> .br n : number of determinants @@ -122,38 +123,26 @@ Documentation H_jj : array of -`lambda_mrcc `_ +`lambda_mrcc `_ cm/ or perturbative 1/Delta_E(m) -`lambda_pert `_ +`lambda_pert `_ cm/ or perturbative 1/Delta_E(m) -`mrcc `_ +`mrcc_dress `_ Undocumented -`mrcc_dress `_ +`mrcc_dress_simple `_ Undocumented -`mrcc_dress_simple `_ - Undocumented - - -`psi_cas_lock `_ - Locks on CAS determinants to fill delta_ij - - -`run `_ - Undocumented - - -`run_mrcc `_ - Undocumented - - -`run_mrcc_test `_ +`psi_ref_lock `_ + Locks on ref determinants to fill delta_ij + + +`run_mrcc `_ Undocumented diff --git a/plugins/MRCC/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f similarity index 97% rename from plugins/MRCC/davidson.irp.f rename to plugins/MRCC_Utils/davidson.irp.f index 588e1831..d9e697a1 100644 --- a/plugins/MRCC/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -35,7 +35,7 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i allocate(H_jj(sze)) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(sze,H_jj,N_det_cas,dets_in,Nint,istate,delta_ii,idx_cas) & + !$OMP SHARED(sze,H_jj,N_det_ref,dets_in,Nint,istate,delta_ii,idx_ref) & !$OMP PRIVATE(i) !$OMP DO SCHEDULE(guided) do i=1,sze @@ -43,8 +43,8 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i enddo !$OMP END DO !$OMP DO SCHEDULE(guided) - do i=1,N_det_cas - H_jj(idx_cas(i)) += delta_ii(i,istate) + do i=1,N_det_ref + H_jj(idx_ref(i)) += delta_ii(i,istate) enddo !$OMP END DO !$OMP END PARALLEL @@ -384,7 +384,7 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) integer, parameter :: block_size = 157 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,hij,j,k,idx,jj,ii,vt) & - !$OMP SHARED(n_det_cas,n_det_non_cas,idx_cas,idx_non_cas,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij) + !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij) !$OMP DO SCHEDULE(static) do i=1,n v_0(i) = H_jj(i) * u_0(i) @@ -409,10 +409,10 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) !$OMP END DO !$OMP DO SCHEDULE(guided) - do ii=1,n_det_cas - i = idx_cas(ii) - do jj = 1, n_det_non_cas - j = idx_non_cas(jj) + do ii=1,n_det_ref + i = idx_ref(ii) + do jj = 1, n_det_non_ref + j = idx_non_ref(jj) vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j) vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i) enddo diff --git a/plugins/MRCC/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f similarity index 75% rename from plugins/MRCC/mrcc_dress.irp.f rename to plugins/MRCC_Utils/mrcc_dress.irp.f index 53f94f5b..8c86f7fa 100644 --- a/plugins/MRCC/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -1,25 +1,25 @@ use omp_lib -BEGIN_PROVIDER [ integer(omp_lock_kind), psi_cas_lock, (psi_det_size) ] +BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ] implicit none BEGIN_DOC - ! Locks on CAS determinants to fill delta_ij + ! Locks on ref determinants to fill delta_ij END_DOC integer :: i do i=1,psi_det_size - call omp_init_lock( psi_cas_lock(i) ) + call omp_init_lock( psi_ref_lock(i) ) enddo END_PROVIDER -subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n_selected,det_buffer,Nint,iproc) +subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc - integer, intent(in) :: Ndet_cas, Ndet_non_cas - double precision, intent(inout) :: delta_ij_(Ndet_cas,Ndet_non_cas,*) - double precision, intent(inout) :: delta_ii_(Ndet_cas,*) + 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 @@ -43,18 +43,18 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) - allocate (dIa_hla(N_states,Ndet_non_cas)) + allocate (dIa_hla(N_states,Ndet_non_ref)) ! |I> ! |alpha> do i_alpha=1,N_tq - call get_excitation_degree_vector(psi_non_cas,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_cas,idx_alpha) + call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) ! |I> - do i_I=1,N_det_cas + do i_I=1,N_det_ref ! Find triples and quadruple grand parents - call get_excitation_degree(tq(1,1,i_alpha),psi_cas(1,1,i_I),degree,Nint) + call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint) if (degree > 4) then cycle endif @@ -65,22 +65,22 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n ! |alpha> do k_sd=1,idx_alpha(0) - call get_excitation_degree(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),degree,Nint) + 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 i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),Nint,hIk) + 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_cas(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + 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_cas(k,1,i_I) - tmp_det(k,2) = psi_cas(k,2,i_I) + 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 @@ -108,10 +108,10 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n dka(i_state) = 0.d0 enddo do l_sd=k_sd+1,idx_alpha(0) - call get_excitation_degree(tmp_det,psi_non_cas(1,1,idx_alpha(l_sd)),degree,Nint) + 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_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) - call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hIl) + 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 @@ -124,28 +124,28 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n enddo do i_state=1,N_states - ci_inv(i_state) = 1.d0/psi_cas_coef(i_I,i_state) + ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state) enddo do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - call i_h_j(tq(1,1,i_alpha),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hla) + call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla) do i_state=1,N_states dIa_hla(i_state,k_sd) = dIa(i_state) * hla enddo enddo - call omp_set_lock( psi_cas_lock(i_I) ) + call omp_set_lock( psi_ref_lock(i_I) ) do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) do i_state=1,N_states delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd) - if(dabs(psi_cas_coef(i_I,i_state)).ge.5.d-5)then - delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state) + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + delta_ii_(i_I,i_state) -= dIa_hla(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 enddo - call omp_unset_lock( psi_cas_lock(i_I) ) + call omp_unset_lock( psi_ref_lock(i_I) ) enddo enddo deallocate (dIa_hla) @@ -157,13 +157,13 @@ end -subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,i_generator,n_selected,det_buffer,Nint,iproc) +subroutine mrcc_dress_simple(delta_ij_non_ref_,Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc - integer, intent(in) :: Ndet_non_cas - double precision, intent(inout) :: delta_ij_non_cas_(Ndet_non_cas,Ndet_non_cas,*) + integer, intent(in) :: Ndet_non_ref + double precision, intent(inout) :: delta_ij_non_ref_(Ndet_non_ref,Ndet_non_ref,*) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,m @@ -184,18 +184,18 @@ subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,i_generator,n_select double precision :: f(N_states) do i=1,N_tq - call get_excitation_degree_vector(psi_non_cas,tq(1,1,i),degree,Nint,Ndet_non_cas,idx) + call get_excitation_degree_vector(psi_non_ref,tq(1,1,i),degree,Nint,Ndet_non_ref,idx) call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa) do m=1,N_states f(m) = 1.d0/(ci_electronic_energy(m)-haa) enddo do k=1,idx(0) - call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(k)),Nint,hka) + call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(k)),Nint,hka) do j=k,idx(0) - call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(j)),Nint,haj) + call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(j)),Nint,haj) do m=1,N_states - delta_ij_non_cas_(idx(k), idx(j),m) += haj*hka* f(m) - delta_ij_non_cas_(idx(j), idx(k),m) += haj*hka* f(m) + delta_ij_non_ref_(idx(k), idx(j),m) += haj*hka* f(m) + delta_ij_non_ref_(idx(j), idx(k),m) += haj*hka* f(m) enddo enddo enddo @@ -231,9 +231,9 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq endif ! Select determinants that are triple or quadruple excitations - ! from the CAS + ! from the ref good = .True. - call get_excitation_degree_vector(psi_cas,det_buffer(1,1,i),degree,Nint,N_det_cas,idx) + call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) do k=1,idx(0) if (degree(k) < 3) then good = .False. diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRCC_Utils/mrcc_general.irp.f new file mode 100644 index 00000000..53650164 --- /dev/null +++ b/plugins/MRCC_Utils/mrcc_general.irp.f @@ -0,0 +1,29 @@ +subroutine run_mrcc + implicit none + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + do while (delta_E > 1.d-10) + iteration += 1 + print *, '===========================' + print *, 'MRCC Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCC energy") + call diagonalize_ci_dressed + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + if (iteration > 20) then + exit + endif + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + call ezfio_set_mrcc_energy(ci_energy_dressed(1)) + call save_wavefunction + +end diff --git a/plugins/MRCC/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f similarity index 84% rename from plugins/MRCC/mrcc_utils.irp.f rename to plugins/MRCC_Utils/mrcc_utils.irp.f index 6ae0926c..d8ed5250 100644 --- a/plugins/MRCC/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -7,16 +7,16 @@ integer :: i,k double precision :: ihpsi(N_states), hii - do i=1,N_det_non_cas - call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coef, N_int, N_det_cas,& - size(psi_cas_coef,1), n_states, ihpsi) - call i_h_j(psi_non_cas(1,1,i),psi_non_cas(1,1,i),N_int,hii) + do i=1,N_det_non_ref + 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) + call i_h_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) do k=1,N_states - lambda_pert(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii) + lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) if (dabs(ihpsi(k)).le.1.d-3) then lambda_mrcc(k,i) = lambda_pert(k,i) else - lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k) + lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi(k) endif enddo enddo @@ -26,17 +26,17 @@ END_PROVIDER -!BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ] +!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ] !implicit none !BEGIN_DOC !! Dressing matrix in SD basis !END_DOC -!delta_ij_non_cas = 0.d0 -!call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas) +!delta_ij_non_ref = 0.d0 +!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref) !END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_ij, (N_det_cas,N_det_non_cas,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_cas,N_states) ] + BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] implicit none BEGIN_DOC ! Dressing matrix in N_det basis @@ -44,7 +44,7 @@ END_PROVIDER integer :: i,j,m delta_ij = 0.d0 delta_ii = 0.d0 - call H_apply_mrcc(delta_ij,delta_ii,N_det_cas,N_det_non_cas) + call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref) END_PROVIDER BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] @@ -59,11 +59,11 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) enddo enddo - do ii = 1, N_det_cas - i =idx_cas(ii) + do ii = 1, N_det_ref + i =idx_ref(ii) h_matrix_dressed(i,i,istate) += delta_ii(ii,istate) - do jj = 1, N_det_non_cas - j =idx_cas(jj) + do jj = 1, N_det_non_ref + j =idx_ref(jj) h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate) h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate) enddo diff --git a/plugins/MRCC_Utils/tree_dependency.png b/plugins/MRCC_Utils/tree_dependency.png new file mode 100644 index 00000000..500e5d43 Binary files /dev/null and b/plugins/MRCC_Utils/tree_dependency.png differ diff --git a/plugins/Perturbation/tree_dependency.png b/plugins/Perturbation/tree_dependency.png index 1b361288..dac64544 100644 Binary files a/plugins/Perturbation/tree_dependency.png and b/plugins/Perturbation/tree_dependency.png differ diff --git a/plugins/Properties/tree_dependency.png b/plugins/Properties/tree_dependency.png index 680a08ee..1ba8d487 100644 Binary files a/plugins/Properties/tree_dependency.png and b/plugins/Properties/tree_dependency.png differ diff --git a/plugins/MRCC/.gitignore b/plugins/Psiref_CAS/.gitignore similarity index 95% rename from plugins/MRCC/.gitignore rename to plugins/Psiref_CAS/.gitignore index 48b27784..d98a4abc 100644 --- a/plugins/MRCC/.gitignore +++ b/plugins/Psiref_CAS/.gitignore @@ -25,5 +25,5 @@ Utils ezfio_interface.irp.f irpf90.make irpf90_entities -mrcc +mrcc_general tags \ No newline at end of file diff --git a/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES b/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..7e790003 --- /dev/null +++ b/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Psiref_Utils diff --git a/plugins/Psiref_CAS/README.rst b/plugins/Psiref_CAS/README.rst new file mode 100644 index 00000000..3d4726e1 --- /dev/null +++ b/plugins/Psiref_CAS/README.rst @@ -0,0 +1,24 @@ +======================= +Psiref_threshold Module +======================= + + +Reference wave function is defined as all determinants with coefficients +such that | c_i/c_max | > threshold. + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. by the `update_README.py` script. + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. by the `update_README.py` script. + +.. image:: tree_dependency.png + +* `Psiref_Utils `_ + diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f new file mode 100644 index 00000000..81df60e9 --- /dev/null +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -0,0 +1,28 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_ref ] + implicit none + BEGIN_DOC + ! CAS wave function, defined from the application of the CAS bitmask on the + ! determinants. idx_cas gives the indice of the CAS determinant in psi_det. + END_DOC + integer :: i,j,k + N_det_ref = N_det_cas + do i=1,N_det_ref + do k=1,N_int + psi_ref(k,1,i) = psi_cas(k,1,i) + psi_ref(k,2,i) = psi_cas(k,2,i) + enddo + idx_ref(i) = idx_cas(i) + enddo + do k=1,N_states + do i=1,N_det_ref + psi_ref_coef(i,k) = psi_cas_coef(i,k) + enddo + enddo + +END_PROVIDER + diff --git a/plugins/Psiref_CAS/tree_dependency.png b/plugins/Psiref_CAS/tree_dependency.png new file mode 100644 index 00000000..1a922bdc Binary files /dev/null and b/plugins/Psiref_CAS/tree_dependency.png differ diff --git a/plugins/Psiref_Utils/.gitignore b/plugins/Psiref_Utils/.gitignore new file mode 100644 index 00000000..d98a4abc --- /dev/null +++ b/plugins/Psiref_Utils/.gitignore @@ -0,0 +1,29 @@ +# Automatically created by $QP_ROOT/scripts/module/module_handler.py +.ninja_deps +.ninja_log +AO_Basis +Bitmask +Determinants +Electrons +Ezfio_files +Generators_full +Hartree_Fock +IRPF90_man +IRPF90_temp +Integrals_Bielec +Integrals_Monoelec +MOGuess +MO_Basis +Makefile +Makefile.depend +Nuclei +Perturbation +Properties +Pseudo +Selectors_full +Utils +ezfio_interface.irp.f +irpf90.make +irpf90_entities +mrcc_general +tags \ No newline at end of file diff --git a/plugins/Psiref_Utils/NEEDED_CHILDREN_MODULES b/plugins/Psiref_Utils/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..e69de29b diff --git a/plugins/Psiref_Utils/README.rst b/plugins/Psiref_Utils/README.rst new file mode 100644 index 00000000..75269412 --- /dev/null +++ b/plugins/Psiref_Utils/README.rst @@ -0,0 +1,15 @@ +=================== +Psiref_utils Module +=================== + + +Utilities related to the use of a reference wave function. This module +needs to be loaded with any psi_ref module. + + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. by the `update_README.py` script. + diff --git a/plugins/Psiref_Utils/psi_ref.irp.f b/plugins/Psiref_Utils/psi_ref.irp.f new file mode 100644 index 00000000..852df2d3 --- /dev/null +++ b/plugins/Psiref_Utils/psi_ref.irp.f @@ -0,0 +1,36 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_ref ] + implicit none + BEGIN_DOC + ! Reference wave function, defined as determinants with coefficients > 0.05 + ! idx_ref gives the indice of the ref determinant in psi_det. + END_DOC + integer :: i, k, l + logical :: good + N_det_ref = 0 + do i=1,N_det + good = .False. + do l = 1, N_states + psi_ref_coef(i,l) = 0.d0 + good = good.or.(dabs(psi_coef(i,l)) > 0.05d0) + enddo + if (good) then + N_det_ref = N_det_ref+1 + do k=1,N_int + psi_ref(k,1,N_det_ref) = psi_det(k,1,i) + psi_ref(k,2,N_det_ref) = psi_det(k,2,i) + enddo + idx_ref(N_det_ref) = i + do k=1,N_states + psi_ref_coef(N_det_ref,k) = psi_coef(i,k) + enddo + endif + enddo + call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference') + +END_PROVIDER + diff --git a/plugins/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f new file mode 100644 index 00000000..c5e80199 --- /dev/null +++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f @@ -0,0 +1,123 @@ +use bitmasks + + + BEGIN_PROVIDER [ integer(bit_kind), psi_ref_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_ref_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! Reference determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_ref, psi_ref, psi_ref_coef, & + psi_ref_sorted_bit, psi_ref_coef_sorted_bit) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_non_ref, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_non_ref ] + implicit none + BEGIN_DOC + ! Set of determinants which are not part of the reference, defined from the application + ! of the reference bitmask on the determinants. + ! idx_non_ref gives the indice of the determinant in psi_det. + END_DOC + integer :: i_non_ref,j,k + integer :: degree + logical :: in_ref + i_non_ref =0 + do k=1,N_det + in_ref = .False. + do j=1,N_det_ref + call get_excitation_degree(psi_ref(1,1,j), psi_det(1,1,k), degree, N_int) + if (degree == 0) then + in_ref = .True. + exit + endif + enddo + if (.not.in_ref) then + double precision :: hij + i_non_ref += 1 + do j=1,N_int + psi_non_ref(j,1,i_non_ref) = psi_det(j,1,k) + psi_non_ref(j,2,i_non_ref) = psi_det(j,2,k) + enddo + do j=1,N_states + psi_non_ref_coef(i_non_ref,j) = psi_coef(k,j) + enddo + idx_non_ref(i_non_ref) = k + endif + enddo + N_det_non_ref = i_non_ref +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_ref_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! Reference determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, & + psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit) + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, H_matrix_ref, (N_det_ref,N_det_ref)] + 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) + H_matrix_ref(i,j) = hij + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [double precision, psi_coef_ref_diagonalized, (N_det_ref,N_states)] +&BEGIN_PROVIDER [double precision, psi_ref_energy_diagonalized, (N_states)] + implicit none + integer :: i,j + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + allocate (eigenvectors(size(H_matrix_ref,1),N_det_ref)) + allocate (eigenvalues(N_det_ref)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_ref,size(H_matrix_ref,1),N_det_ref) + 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) + enddo + enddo + + + END_PROVIDER + + BEGIN_PROVIDER [double precision, psi_ref_energy, (N_states)] + implicit none + integer :: i,j,k + double precision :: hij,norm,u_dot_v + psi_ref_energy = 0.d0 + + + do k = 1, N_states + norm = 0.d0 + do i = 1, N_det_ref + norm += psi_ref_coef(i,k) * psi_ref_coef(i,k) + do j = 1, N_det_ref + psi_ref_energy(k) += psi_ref_coef(i,k) * psi_ref_coef(j,k) * H_matrix_ref(i,j) + enddo + enddo + psi_ref_energy(k) = psi_ref_energy(k) /norm + enddo + +END_PROVIDER + + + + diff --git a/plugins/Psiref_Utils/tree_dependency.png b/plugins/Psiref_Utils/tree_dependency.png new file mode 100644 index 00000000..20482ad2 Binary files /dev/null and b/plugins/Psiref_Utils/tree_dependency.png differ diff --git a/plugins/Psiref_threshold/.gitignore b/plugins/Psiref_threshold/.gitignore new file mode 100644 index 00000000..d98a4abc --- /dev/null +++ b/plugins/Psiref_threshold/.gitignore @@ -0,0 +1,29 @@ +# Automatically created by $QP_ROOT/scripts/module/module_handler.py +.ninja_deps +.ninja_log +AO_Basis +Bitmask +Determinants +Electrons +Ezfio_files +Generators_full +Hartree_Fock +IRPF90_man +IRPF90_temp +Integrals_Bielec +Integrals_Monoelec +MOGuess +MO_Basis +Makefile +Makefile.depend +Nuclei +Perturbation +Properties +Pseudo +Selectors_full +Utils +ezfio_interface.irp.f +irpf90.make +irpf90_entities +mrcc_general +tags \ No newline at end of file diff --git a/plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES b/plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..7e790003 --- /dev/null +++ b/plugins/Psiref_threshold/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Psiref_Utils diff --git a/plugins/Psiref_threshold/README.rst b/plugins/Psiref_threshold/README.rst new file mode 100644 index 00000000..3d4726e1 --- /dev/null +++ b/plugins/Psiref_threshold/README.rst @@ -0,0 +1,24 @@ +======================= +Psiref_threshold Module +======================= + + +Reference wave function is defined as all determinants with coefficients +such that | c_i/c_max | > threshold. + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. by the `update_README.py` script. + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. by the `update_README.py` script. + +.. image:: tree_dependency.png + +* `Psiref_Utils `_ + diff --git a/plugins/Psiref_threshold/psi_ref.irp.f b/plugins/Psiref_threshold/psi_ref.irp.f new file mode 100644 index 00000000..5e722822 --- /dev/null +++ b/plugins/Psiref_threshold/psi_ref.irp.f @@ -0,0 +1,38 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_ref ] + implicit none + BEGIN_DOC + ! Reference wave function, defined as determinants with coefficients > 0.05 + ! idx_ref gives the indice of the ref determinant in psi_det. + END_DOC + integer :: i, k, l + logical :: good + double precision, parameter :: threshold=0.05d0 + N_det_ref = 0 + t = threshold * abs_psi_coef_max + do i=1,N_det + good = .False. + do l = 1, N_states + psi_ref_coef(i,l) = 0.d0 + good = good.or.(dabs(psi_coef(i,l)) > t) + enddo + if (good) then + N_det_ref = N_det_ref+1 + do k=1,N_int + psi_ref(k,1,N_det_ref) = psi_det(k,1,i) + psi_ref(k,2,N_det_ref) = psi_det(k,2,i) + enddo + idx_ref(N_det_ref) = i + do k=1,N_states + psi_ref_coef(N_det_ref,k) = psi_coef(i,k) + enddo + endif + enddo + call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference') + +END_PROVIDER + diff --git a/plugins/Psiref_threshold/tree_dependency.png b/plugins/Psiref_threshold/tree_dependency.png new file mode 100644 index 00000000..9c2088e1 Binary files /dev/null and b/plugins/Psiref_threshold/tree_dependency.png differ diff --git a/plugins/QmcChem/tree_dependency.png b/plugins/QmcChem/tree_dependency.png index 9160ab75..3b844e5c 100644 Binary files a/plugins/QmcChem/tree_dependency.png and b/plugins/QmcChem/tree_dependency.png differ diff --git a/plugins/Selectors_full/tree_dependency.png b/plugins/Selectors_full/tree_dependency.png index b6378b39..f49b2e9a 100644 Binary files a/plugins/Selectors_full/tree_dependency.png and b/plugins/Selectors_full/tree_dependency.png differ diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index 52f06f7b..f2b261ad 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -132,7 +132,7 @@ def dict_module_genelogy_path(d_module_genelogy): def get_l_module_with_ezfio_cfg(): """ - Return all the module who have a EZFIO.cfg + Return all the modules that have a EZFIO.cfg """ from os import listdir from os.path import isfile @@ -739,7 +739,7 @@ def create_build_ninja_module(path_module): l_string += ["rule make_all_binaries", " command = ninja -f {0}".format(ROOT_BUILD_NINJA), - " pool = console", " description = Compile all the module", + " pool = console", " description = Compiling all modules", ""] l_string += ["rule make_clean", " command = module_handler.py clean {0}".format(path_module.rel), @@ -768,7 +768,7 @@ def create_build_ninja_global(): l_string += ["rule make_all_binaries", " command = ninja -f {0}".format(ROOT_BUILD_NINJA), - " pool = console", " description = Compile all the module", + " pool = console", " description = Compiling all modules", ""] l_string += ["rule make_clean", diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index ff63a18d..f44fb097 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -30,7 +30,7 @@ p = re.compile(ur'\|(\d+)><\d+\|') def get_pseudo_str(l_atom): """ - Run EMSL_local for geting the str of the speudo potential + Run EMSL_local for getting the str of the pseudo potential str_ele : Element Symbol: Na diff --git a/scripts/qp_set_frozen_core.py b/scripts/qp_set_frozen_core.py index 843df4b3..3f95a9e6 100755 --- a/scripts/qp_set_frozen_core.py +++ b/scripts/qp_set_frozen_core.py @@ -3,7 +3,7 @@ import os import sys -sys.path = [ os.environ["QP_ROOT"]+"/EZFIO/Python" ] + sys.path +sys.path = [ os.environ["QP_ROOT"]+"/install/EZFIO/Python" ] + sys.path from ezfio import ezfio ezfio.set_filename(sys.argv[1]) diff --git a/src/AO_Basis/tree_dependency.png b/src/AO_Basis/tree_dependency.png index d69b7130..acaeb7af 100644 Binary files a/src/AO_Basis/tree_dependency.png and b/src/AO_Basis/tree_dependency.png differ diff --git a/src/Bitmask/tree_dependency.png b/src/Bitmask/tree_dependency.png index dfa495c3..24ce3397 100644 Binary files a/src/Bitmask/tree_dependency.png and b/src/Bitmask/tree_dependency.png differ diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index 966dd14b..823cdd39 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -47,6 +47,14 @@ Documentation Needed for diag_H_mat_elem +`abs_psi_coef_max `_ + Max and min values of the coefficients + + +`abs_psi_coef_min `_ + Max and min values of the coefficients + + `ac_operator `_ Needed for diag_H_mat_elem @@ -120,11 +128,11 @@ Documentation Initial guess vectors are not necessarily orthonormal -`connected_to_ref `_ +`connected_to_ref `_ Undocumented -`connected_to_ref_by_mono `_ +`connected_to_ref_by_mono `_ Undocumented @@ -133,7 +141,7 @@ Documentation After calling this subroutine, N_det, psi_det and psi_coef need to be touched -`create_wf_of_psi_svd_matrix `_ +`create_wf_of_psi_svd_matrix `_ Matrix of wf coefficients. Outer product of alpha and beta determinants @@ -209,7 +217,7 @@ Documentation det_coef -`det_connections `_ +`det_connections `_ Build connection proxy between determinants @@ -339,7 +347,7 @@ Documentation Determinants are taken from the psi_det_sorted_ab array -`generate_all_alpha_beta_det_products `_ +`generate_all_alpha_beta_det_products `_ Create a wave function from all possible alpha x beta determinants @@ -359,15 +367,15 @@ Documentation Applies get_excitation_degree to an array of determinants -`get_index_in_psi_det_alpha_unique `_ +`get_index_in_psi_det_alpha_unique `_ Returns the index of the determinant in the ``psi_det_alpha_unique`` array -`get_index_in_psi_det_beta_unique `_ +`get_index_in_psi_det_beta_unique `_ Returns the index of the determinant in the ``psi_det_beta_unique`` array -`get_index_in_psi_det_sorted_bit `_ +`get_index_in_psi_det_sorted_bit `_ Returns the index of the determinant in the ``psi_det_sorted_bit`` array @@ -383,7 +391,11 @@ Documentation Returns -`get_s2_u0 `_ +`get_s2_u0 `_ + Undocumented + + +`get_s2_u0_old `_ Undocumented @@ -480,7 +492,7 @@ Documentation .br -`is_in_wavefunction `_ +`is_in_wavefunction `_ True if the determinant ``det`` is in the wave function @@ -500,7 +512,7 @@ Documentation Energy of the reference bitmask used in Slater rules -`n_con_int `_ +`n_con_int `_ Number of integers to represent the connections between determinants @@ -508,11 +520,11 @@ Documentation Number of determinants in the wave function -`n_det_alpha_unique `_ +`n_det_alpha_unique `_ Unique alpha determinants -`n_det_beta_unique `_ +`n_det_beta_unique `_ Unique beta determinants @@ -565,7 +577,7 @@ Documentation Energy of the reference bitmask used in Slater rules -`occ_pattern_search_key `_ +`occ_pattern_search_key `_ Return an integer*8 corresponding to a determinant index for searching @@ -654,11 +666,19 @@ Documentation Undocumented +`psi_coef_max `_ + Max and min values of the coefficients + + +`psi_coef_min `_ + Max and min values of the coefficients + + `psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_ab `_ +`psi_coef_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate @@ -677,19 +697,19 @@ Documentation is empty -`psi_det_alpha `_ +`psi_det_alpha `_ List of alpha determinants of psi_det -`psi_det_alpha_unique `_ +`psi_det_alpha_unique `_ Unique alpha determinants -`psi_det_beta `_ +`psi_det_beta `_ List of beta determinants of psi_det -`psi_det_beta_unique `_ +`psi_det_beta_unique `_ Unique beta determinants @@ -701,7 +721,7 @@ Documentation Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_ab `_ +`psi_det_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate @@ -715,7 +735,7 @@ Documentation function. -`psi_det_sorted_next_ab `_ +`psi_det_sorted_next_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate @@ -750,31 +770,31 @@ Documentation psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation -`psi_svd_alpha `_ +`psi_svd_alpha `_ SVD wave function -`psi_svd_beta `_ +`psi_svd_beta `_ SVD wave function -`psi_svd_coefs `_ +`psi_svd_coefs `_ SVD wave function -`psi_svd_matrix `_ +`psi_svd_matrix `_ Matrix of wf coefficients. Outer product of alpha and beta determinants -`psi_svd_matrix_columns `_ +`psi_svd_matrix_columns `_ Matrix of wf coefficients. Outer product of alpha and beta determinants -`psi_svd_matrix_rows `_ +`psi_svd_matrix_rows `_ Matrix of wf coefficients. Outer product of alpha and beta determinants -`psi_svd_matrix_values `_ +`psi_svd_matrix_values `_ Matrix of wf coefficients. Outer product of alpha and beta determinants @@ -782,7 +802,7 @@ Documentation Undocumented -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file @@ -835,15 +855,15 @@ Documentation Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`save_wavefunction_general `_ +`save_wavefunction_general `_ Save the wave function into the EZFIO file -`save_wavefunction_unsorted `_ +`save_wavefunction_unsorted `_ Save the wave function into the EZFIO file @@ -857,7 +877,7 @@ Documentation for a given couple of hole/particle excitations i. -`sort_dets_by_3_highest_electrons `_ +`sort_dets_by_3_highest_electrons `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate @@ -890,6 +910,6 @@ Documentation Thresholds on selectors (fraction of the norm) -`write_spindeterminants `_ +`write_spindeterminants `_ Undocumented diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index 2d40b621..3f4b62d2 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -11,6 +11,7 @@ integer*8 function det_search_key(det,Nint) do i=2,Nint det_search_key = ieor(det_search_key,iand(det(i,1),det(i,2))) enddo + det_search_key = iand(huge(det(1,1)),det_search_key) end @@ -27,6 +28,7 @@ integer*8 function occ_pattern_search_key(det,Nint) do i=2,Nint occ_pattern_search_key = ieor(occ_pattern_search_key,iand(det(i,1),det(i,2))) enddo + occ_pattern_search_key = iand(huge(det(1,1)),occ_pattern_search_key) end diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 8cc545f5..6834a745 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -446,6 +446,24 @@ subroutine filter_3_highest_electrons( det_in, det_out, Nint ) enddo end + BEGIN_PROVIDER [ double precision, psi_coef_max, (N_states) ] +&BEGIN_PROVIDER [ double precision, psi_coef_min, (N_states) ] +&BEGIN_PROVIDER [ double precision, abs_psi_coef_max, (N_states) ] +&BEGIN_PROVIDER [ double precision, abs_psi_coef_min, (N_states) ] + implicit none + BEGIN_DOC + ! Max and min values of the coefficients + END_DOC + integer:: i + do i=1,N_states + psi_coef_min(i) = minval(psi_coef(:,i)) + psi_coef_max(i) = maxval(psi_coef(:,i)) + abs_psi_coef_min(i) = dabs(psi_coef_min(i)) + abs_psi_coef_max(i) = dabs(psi_coef_max(i)) + enddo + +END_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_ab, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ double precision, psi_coef_sorted_ab, (N_det,N_states) ] &BEGIN_PROVIDER [ integer, psi_det_sorted_next_ab, (2,psi_det_size) ] diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index fd42fe51..72c1b9aa 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -79,7 +79,7 @@ BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] END_PROVIDER -subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) +subroutine get_s2_u0_old(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) implicit none use bitmasks integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax) @@ -89,32 +89,60 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) integer :: i,j,l double precision :: s2_tmp s2 = 0.d0 -!print*,'IN get_s2_u0' -!print*,'n,nmax = ',n,nmax - double precision :: accu - accu = 0.d0 - do i = 1,n - accu += psi_coefs_tmp(i) * psi_coefs_tmp(i) -! print*,'psi_coef = ',psi_coefs_tmp(i) - enddo -!print*,'accu = ',accu -!print*,'' !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) REDUCTION(+:s2) SCHEDULE(dynamic) do i=1,n do j=i+1,n call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int) - s2 += (psi_coefs_tmp(i)+psi_coefs_tmp(i))*psi_coefs_tmp(j)*s2_tmp -! s2 = s2 + 2.d0 * psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp + s2 += psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp enddo enddo !$OMP END PARALLEL DO + s2 = s2+s2 do i=1,n call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) s2 += psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp enddo s2 += S_z2_Sz -!print*,'s2 = ',s2 -!print*,'' +end + +subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax) + integer, intent(in) :: n,nmax + double precision, intent(in) :: psi_coefs_tmp(nmax) + double precision, intent(out) :: s2 + double precision :: s2_tmp + integer :: i,j,l,jj + integer, allocatable :: idx(:) + s2 = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,s2_tmp,idx) & + !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold)& + !$OMP REDUCTION(+:s2) + allocate(idx(0:n)) + !$OMP DO SCHEDULE(dynamic) + do i=1,n + idx(0) = i + call filter_connected_davidson(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx) + do jj=1,idx(0) + j = idx(jj) + if ( dabs(psi_coefs_tmp(j)) + dabs(psi_coefs_tmp(i)) & + > davidson_threshold ) then + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int) + s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp + endif + enddo + enddo + !$OMP END DO + deallocate(idx) + !$OMP END PARALLEL + s2 = s2+s2 + do i=1,n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) + s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp + enddo + s2 = s2 + S_z2_Sz end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 7d431879..a8b394d9 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1092,10 +1092,9 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) ASSERT (Nint == N_int) ASSERT (n>0) PROVIDE ref_bitmask_energy - integer, parameter :: block_size = 157 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,hij,j,k,idx,jj,vt) & - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0) + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold) !$OMP DO SCHEDULE(static) do i=1,n v_0(i) = H_jj(i) * u_0(i) @@ -1109,7 +1108,7 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) do jj=1,idx(0) j = idx(jj) - if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then + if ( dabs(u_0(j)) + dabs(u_0(i)) > davidson_threshold ) then call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) vt (i) = vt (i) + hij*u_0(j) vt (j) = vt (j) + hij*u_0(i) diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index 5daf4e2c..ad9d9960 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -14,11 +14,13 @@ integer*8 function spin_det_search_key(det,Nint) END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: det(Nint) + integer(bit_kind), parameter :: unsigned_shift = not(huge(1_bit_kind)) ! 100...00 integer :: i spin_det_search_key = det(1) do i=2,Nint spin_det_search_key = ieor(spin_det_search_key,det(i)) enddo + spin_det_search_key = spin_det_search_key-unsigned_shift end @@ -346,53 +348,6 @@ subroutine write_spindeterminants call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_svd_matrix_rows) call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_svd_matrix_columns) -! integer :: n_svd_coefs -! double precision :: norm, f -! f = 1.d0/dble(N_states) -! norm = 1.d0 -! do n_svd_coefs=1,N_det_alpha_unique -! do k=1,N_states -! norm -= psi_svd_coefs(n_svd_coefs,k)*psi_svd_coefs(n_svd_coefs,k) -! enddo -! if (norm < 1.d-4) then -! exit -! endif -! enddo -! n_svd_coefs -= 1 -! call ezfio_set_spindeterminants_n_svd_coefs(n_svd_coefs) -! -! double precision, allocatable :: dtmp(:,:,:) -! allocate(dtmp(N_det_alpha_unique,n_svd_coefs,N_states)) -! do k=1,N_states -! do j=1,n_svd_coefs -! do i=1,N_det_alpha_unique -! dtmp(i,j,k) = psi_svd_alpha(i,j,k) -! enddo -! enddo -! enddo -! call ezfio_set_spindeterminants_psi_svd_alpha(dtmp) -! deallocate(dtmp) -! -! allocate(dtmp(N_det_beta_unique,n_svd_coefs,N_states)) -! do k=1,N_states -! do j=1,n_svd_coefs -! do i=1,N_det_beta_unique -! dtmp(i,j,k) = psi_svd_beta(i,j,k) -! enddo -! enddo -! enddo -! call ezfio_set_spindeterminants_psi_svd_beta(dtmp) -! deallocate(dtmp) -! -! allocate(dtmp(n_svd_coefs,N_states,1)) -! do k=1,N_states -! do j=1,n_svd_coefs -! dtmp(j,k,1) = psi_svd_coefs(j,k) -! enddo -! enddo -! call ezfio_set_spindeterminants_psi_svd_coefs(dtmp) -! deallocate(dtmp) - end @@ -419,28 +374,6 @@ BEGIN_PROVIDER [ double precision, psi_svd_matrix_values, (N_det,N_states) ] PROVIDE psi_coef_sorted_bit -! l=0 -! do j=1,N_det_beta_unique -! do k=1,N_int -! tmp_det(k,2) = psi_det_beta_unique(k,j) -! enddo -! do i=1,N_det_alpha_unique -! do k=1,N_int -! tmp_det(k,1) = psi_det_alpha_unique(k,i) -! enddo -! idx = get_index_in_psi_det_sorted_bit(tmp_det,N_int) -! if (idx > 0) then -! l += 1 -! psi_svd_matrix_rows(l) = i -! psi_svd_matrix_columns(l) = j -! do k=1,N_states -! psi_svd_matrix_values(l,k) = psi_coef_sorted_bit(idx,k) -! enddo -! endif -! enddo -! enddo -! ASSERT (l == N_det) - integer, allocatable :: iorder(:), to_sort(:) integer, external :: get_index_in_psi_det_alpha_unique integer, external :: get_index_in_psi_det_beta_unique diff --git a/src/Determinants/tree_dependency.png b/src/Determinants/tree_dependency.png index 8731537a..f9eb10c3 100644 Binary files a/src/Determinants/tree_dependency.png and b/src/Determinants/tree_dependency.png differ diff --git a/src/Electrons/tree_dependency.png b/src/Electrons/tree_dependency.png index 8ac27f56..b90a1f83 100644 Binary files a/src/Electrons/tree_dependency.png and b/src/Electrons/tree_dependency.png differ diff --git a/src/Ezfio_files/tree_dependency.png b/src/Ezfio_files/tree_dependency.png index d54e695c..48f53991 100644 Binary files a/src/Ezfio_files/tree_dependency.png and b/src/Ezfio_files/tree_dependency.png differ diff --git a/src/Integrals_Bielec/tree_dependency.png b/src/Integrals_Bielec/tree_dependency.png index 30d4bbdb..4161fd0a 100644 Binary files a/src/Integrals_Bielec/tree_dependency.png and b/src/Integrals_Bielec/tree_dependency.png differ diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90 index 66a13419..bd00dc51 100644 --- a/src/Integrals_Monoelec/pseudopot.f90 +++ b/src/Integrals_Monoelec/pseudopot.f90 @@ -702,7 +702,7 @@ integer n_k(klocmax_max) double precision a(3),g_a,b(3),g_b,c(3),d(3) integer n_a(3),n_b(3),ntotA,ntotB,ntot,m integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p -double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0 +double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0,coef double precision,allocatable :: array_R_loc(:,:,:) double precision,allocatable :: array_coefs(:,:,:,:,:,:) double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg @@ -713,100 +713,119 @@ double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) arg=g_a*ac**2+g_b*bc**2 if(arg.gt.-dlog(10.d-20))then - Vloc=0.d0 - return + Vloc=0.d0 + return endif - + ntotA=n_a(1)+n_a(2)+n_a(3) ntotB=n_b(1)+n_b(2)+n_b(3) ntot=ntotA+ntotB - + if(ac.eq.0.d0.and.bc.eq.0.d0)then - accu=0.d0 - - do k=1,klocmax - accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k)) - enddo - Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3)) - !bigI frequently is null - return + accu=0.d0 + + do k=1,klocmax + accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k)) + enddo + Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3)) + !bigI frequently is null + return endif - + freal=dexp(-g_a*ac**2-g_b*bc**2) d2 = 0.d0 do i=1,3 - d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i)) - d2=d2+d(i)*d(i) + d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i)) + d2=d2+d(i)*d(i) enddo d2=dsqrt(d2) dreal=2.d0*d2 - - - theta_DC0=dacos(d(3)/d2) - phi_DC0=datan2(d(2)/d2,d(1)/d2) - - if (isnan(theta_DC0).or.isnan(phi_DC0)) then - print *, 'NaN in /src/Integrals_Monoelec/pseudopot.f90 at line 449.' - print *, 'Try to break symmetry in your molecule (1.-16 is OK).' - stop 1 - endif - -allocate (array_R_loc(-2:ntot_max+klocmax_max,klocmax_max,0:ntot_max)) -allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) - + + + allocate (array_R_loc(-2:ntot_max+klocmax_max,klocmax_max,0:ntot_max)) + allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) + do ktot=-2,ntotA+ntotB+klocmax - do l=0,ntot - do k=1,klocmax - array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) - enddo - enddo - enddo - - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3) & - *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) & - *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p) & - *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) - enddo - enddo - enddo - enddo - enddo - enddo - - accu=0.d0 - do k=1,klocmax - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - - do l=0,ntot - do m=-l,l - prod=ylm(l,m,theta_DC0,phi_DC0)*array_coefs(k1,k2,k3,k1p,k2p,k3p) & - *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) - ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) - accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) + do l=0,ntot + do k=1,klocmax + array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) + enddo enddo - enddo - - enddo - enddo - enddo - enddo - enddo - enddo enddo + + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)& + *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)& + *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)& + *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) + enddo + enddo + enddo + enddo + enddo + enddo + + + accu=0.d0 + if(d2 == 0.d0)then + l=0 + m=0 + coef=1.d0/dsqrt(4.d0*dacos(-1.d0)) + do k=1,klocmax + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) & + *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) + ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) + accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + else + theta_DC0=dacos(d(3)/d2) + phi_DC0=datan2(d(2)/d2,d(1)/d2) + + do k=1,klocmax + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + do l=0,ntot + do m=-l,l + coef=ylm(l,m,theta_DC0,phi_DC0) + prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) & + *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) + ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) + accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + endif Vloc=f*accu - + deallocate (array_R_loc) deallocate (array_coefs) end diff --git a/src/Integrals_Monoelec/tree_dependency.png b/src/Integrals_Monoelec/tree_dependency.png index 7431ba2d..f56c1e77 100644 Binary files a/src/Integrals_Monoelec/tree_dependency.png and b/src/Integrals_Monoelec/tree_dependency.png differ diff --git a/src/MOGuess/tree_dependency.png b/src/MOGuess/tree_dependency.png index 994660a3..f33b4bb3 100644 Binary files a/src/MOGuess/tree_dependency.png and b/src/MOGuess/tree_dependency.png differ diff --git a/src/MO_Basis/tree_dependency.png b/src/MO_Basis/tree_dependency.png index 9c3e1be6..c8086369 100644 Binary files a/src/MO_Basis/tree_dependency.png and b/src/MO_Basis/tree_dependency.png differ diff --git a/src/Nuclei/tree_dependency.png b/src/Nuclei/tree_dependency.png index 7b3ed613..72cfaeee 100644 Binary files a/src/Nuclei/tree_dependency.png and b/src/Nuclei/tree_dependency.png differ diff --git a/src/Pseudo/tree_dependency.png b/src/Pseudo/tree_dependency.png index a2dc7c27..5a9af5ae 100644 Binary files a/src/Pseudo/tree_dependency.png and b/src/Pseudo/tree_dependency.png differ diff --git a/src/Utils/tree_dependency.png b/src/Utils/tree_dependency.png index f0aa66a7..38b04785 100644 Binary files a/src/Utils/tree_dependency.png and b/src/Utils/tree_dependency.png differ