diff --git a/.travis.yml b/.travis.yml index 18a13949..24687b97 100644 --- a/.travis.yml +++ b/.travis.yml @@ -24,7 +24,7 @@ python: script: - ./configure --production ./config/gfortran.cfg - - source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles + - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD mrcepa0 All_singles - source ./quantum_package.rc ; ninja - source ./quantum_package.rc ; cd ocaml ; make ; cd - - source ./quantum_package.rc ; cd tests ; ./run_tests.sh #-v diff --git a/README.md b/README.md index 5372b7ac..bb63b691 100644 --- a/README.md +++ b/README.md @@ -7,11 +7,14 @@ Set of quantum chemistry programs and libraries. For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>), or below for the installation instructions. + + Demo ==== [![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_295x166.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo") +[![Frozen-core Full-CI energy of Ti](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png)](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png "Convergence of Ti in cc-pv{DTQ}Z") # Installation @@ -155,7 +158,7 @@ Program exited with code 139. #### Why ? -It's caused when we call the DGEM routine of LAPACK. +It's caused when we call the DGEMM routine of LAPACK. ##### Fix diff --git a/data/Titanium.png b/data/Titanium.png new file mode 100644 index 00000000..871babd4 Binary files /dev/null and b/data/Titanium.png differ diff --git a/plugins/FOBOCI/H_apply_dressed_autonom.irp.f b/plugins/FOBOCI/H_apply_dressed_autonom.irp.f index 69929afd..abe6ef2e 100644 --- a/plugins/FOBOCI/H_apply_dressed_autonom.irp.f +++ b/plugins/FOBOCI/H_apply_dressed_autonom.irp.f @@ -273,7 +273,7 @@ subroutine H_apply_dressed_pert_monoexc(key_in, hole_1,particl_1,i_generator,ipr integer,parameter :: size_max = 3072 integer, intent(in) :: Ndet_generators - double precision, intent(in) :: E_ref + double precision, intent(inout) :: E_ref double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators) integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) @@ -438,7 +438,7 @@ subroutine H_apply_dressed_pert(delta_ij_generators_, Ndet_generators,psi_det_g integer, intent(in) :: Ndet_generators - double precision, intent(in) :: E_ref + double precision, intent(inout) :: E_ref double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators) integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) diff --git a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES index e40934be..16fce081 100644 --- a/plugins/FOBOCI/NEEDED_CHILDREN_MODULES +++ b/plugins/FOBOCI/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_no_sorted Hartree_Fock Davidson +Perturbation Selectors_no_sorted Hartree_Fock Davidson CISD diff --git a/plugins/FOBOCI/hcc_1h1p.irp.f b/plugins/FOBOCI/hcc_1h1p.irp.f index 66cf2fd4..bad073db 100644 --- a/plugins/FOBOCI/hcc_1h1p.irp.f +++ b/plugins/FOBOCI/hcc_1h1p.irp.f @@ -15,11 +15,10 @@ subroutine routine call diagonalize_CI call test_hcc call test_mulliken -! call SC2_1h1p(psi_det,psi_coef,energies, & -! diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) allocate(H_matrix(N_det,N_det)) - call SC2_1h1p_full(psi_det,psi_coef,energies, & - H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) + stop 'SC2_1h1p_full is not in the git!' +! call SC2_1h1p_full(psi_det,psi_coef,energies, & +! H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) deallocate(H_matrix) integer :: i,j double precision :: accu,coef_hf diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index 42337258..c81b1266 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/plugins/Full_CI_ZMQ/fci_zmq.irp.f @@ -46,14 +46,16 @@ program fci_zmq PROVIDE psi_det_sorted call diagonalize_CI + call save_wavefunction if (N_det > N_det_max) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted N_det = N_det_max soft_touch N_det psi_det psi_coef + call diagonalize_CI + call save_wavefunction endif - call save_wavefunction print *, 'N_det = ', N_det print *, 'N_states = ', N_states @@ -79,17 +81,13 @@ program fci_zmq E_CI_before(1:N_states) = CI_energy(1:N_states) call ezfio_set_full_ci_energy(CI_energy) enddo - if (N_det > N_det_max) then - N_det = N_det_max - touch N_det psi_det psi_coef - call diagonalize_CI - endif + if(do_pt2_end)then print*,'Last iteration only to compute the PT2' threshold_selectors = 1.d0 threshold_generators = 0.9999d0 E_CI_before(1:N_states) = CI_energy(1:N_states) - call ZMQ_selection(1, pt2) + call ZMQ_selection(0, pt2) print *, 'Final step' print *, 'N_det = ', N_det print *, 'N_states = ', N_states @@ -108,7 +106,7 @@ end -subroutine ZMQ_selection(N, pt2) +subroutine ZMQ_selection(N_in, pt2) use f77_zmq use selection_types @@ -116,13 +114,14 @@ subroutine ZMQ_selection(N, pt2) character*(512) :: task integer(ZMQ_PTR) :: zmq_to_qp_run_socket - integer, intent(in) :: N + integer, intent(in) :: N_in type(selection_buffer) :: b - integer :: i + integer :: i, N integer, external :: omp_get_thread_num double precision, intent(out) :: pt2(N_states) + N = max(N_in,1) provide nproc provide ci_electronic_energy call new_parallel_job(zmq_to_qp_run_socket,"selection") @@ -147,16 +146,18 @@ subroutine ZMQ_selection(N, pt2) if (i==0) then call selection_collector(b, pt2) else - call selection_dressing_slave_inproc(i) + call selection_slave_inproc(i) endif !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, 'selection') - call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN - call copy_H_apply_buffer_to_wf() + if (N_in > 0) then + call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN + call copy_H_apply_buffer_to_wf() + endif end subroutine -subroutine selection_dressing_slave_inproc(i) +subroutine selection_slave_inproc(i) implicit none integer, intent(in) :: i diff --git a/plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f b/plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f new file mode 100644 index 00000000..6e4cf44f --- /dev/null +++ b/plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f @@ -0,0 +1,107 @@ +program selection_slave + implicit none + BEGIN_DOC +! Helper program to compute the PT2 in distributed mode. + END_DOC + + read_wf = .False. + SOFT_TOUCH read_wf + call provide_everything + call switch_qp_run_to_master + call run_wf +end + +subroutine provide_everything + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context mo_mono_elec_integral +! PROVIDE ci_electronic_energy mo_tot_num N_int +end + +subroutine run_wf + use f77_zmq + implicit none + + integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + double precision :: energy(N_states_diag) + character*(64) :: states(2) + integer :: rc, i + + call provide_everything + + zmq_context = f77_zmq_ctx_new () + states(1) = 'selection' + states(2) = 'davidson' + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + + do + + call wait_for_states(states,zmq_state,2) + + if(trim(zmq_state) == 'Stopped') then + + exit + + else if (trim(zmq_state) == 'selection') then + + ! Selection + ! --------- + + print *, 'Selection' + call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states_diag) + + !$OMP PARALLEL PRIVATE(i) + i = omp_get_thread_num() + call selection_slave_tcp(i, energy) + !$OMP END PARALLEL + print *, 'Selection done' + + else if (trim(zmq_state) == 'davidson') then + + ! Davidson + ! -------- + + print *, 'Davidson' + call davidson_miniserver_get() + !$OMP PARALLEL PRIVATE(i) + i = omp_get_thread_num() + call davidson_slave_tcp(i) + !$OMP END PARALLEL + print *, 'Davidson done' + + endif + + end do +end + +subroutine update_energy(energy) + implicit none + double precision, intent(in) :: energy(N_states_diag) + BEGIN_DOC +! Update energy when it is received from ZMQ + END_DOC + integer :: j,k + do j=1,N_states + do k=1,N_det + CI_eigenvectors(k,j) = psi_coef(k,j) + enddo + enddo + call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int) + if (.True.) then + do k=1,size(ci_electronic_energy) + ci_electronic_energy(k) = energy(k) + enddo + TOUCH ci_electronic_energy CI_eigenvectors_s2 CI_eigenvectors + endif + + call write_double(6,ci_energy,'Energy') +end + +subroutine selection_slave_tcp(i,energy) + implicit none + double precision, intent(in) :: energy(N_states_diag) + integer, intent(in) :: i + + call run_selection_slave(0,i,energy) +end + diff --git a/plugins/Full_CI_ZMQ/selection_slave.irp.f b/plugins/Full_CI_ZMQ/selection_slave.irp.f index 31f42e7e..06bcf533 100644 --- a/plugins/Full_CI_ZMQ/selection_slave.irp.f +++ b/plugins/Full_CI_ZMQ/selection_slave.irp.f @@ -23,20 +23,19 @@ subroutine run_wf integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket integer(ZMQ_PTR) :: zmq_to_qp_run_socket double precision :: energy(N_states_diag) - character*(64) :: states(2) + character*(64) :: states(1) integer :: rc, i call provide_everything zmq_context = f77_zmq_ctx_new () states(1) = 'selection' - states(2) = 'davidson' zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() do - call wait_for_states(states,zmq_state,2) + call wait_for_states(states,zmq_state,1) if(trim(zmq_state) == 'Stopped') then @@ -52,24 +51,10 @@ subroutine run_wf !$OMP PARALLEL PRIVATE(i) i = omp_get_thread_num() - call selection_dressing_slave_tcp(i, energy) + call selection_slave_tcp(i, energy) !$OMP END PARALLEL print *, 'Selection done' - else if (trim(zmq_state) == 'davidson') then - - ! Davidson - ! -------- - - print *, 'Davidson' - call davidson_miniserver_get() - - !!$OMP PARALLEL PRIVATE(i) - !i = omp_get_thread_num() - call davidson_slave_tcp(0) - !!$OMP END PARALLEL - print *, 'Davidson done' - endif end do @@ -98,7 +83,7 @@ subroutine update_energy(energy) call write_double(6,ci_energy,'Energy') end -subroutine selection_dressing_slave_tcp(i,energy) +subroutine selection_slave_tcp(i,energy) implicit none double precision, intent(in) :: energy(N_states_diag) integer, intent(in) :: i diff --git a/plugins/MRCC_CASSD/.gitignore b/plugins/MRCC_CASSD/.gitignore deleted file mode 100644 index 97bd070c..00000000 --- a/plugins/MRCC_CASSD/.gitignore +++ /dev/null @@ -1,34 +0,0 @@ -# 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 -ZMQ -ezfio_interface.irp.f -irpf90.make -irpf90_entities -mrcc_cassd -mrcc_noiter -tags \ No newline at end of file diff --git a/plugins/MRCC_CASSD/EZFIO.cfg b/plugins/MRCC_CASSD/EZFIO.cfg deleted file mode 100644 index 17ee7f29..00000000 --- a/plugins/MRCC_CASSD/EZFIO.cfg +++ /dev/null @@ -1,17 +0,0 @@ -[energy] -type: double precision -doc: Calculated energy -interface: ezfio - -[thresh_mrcc] -type: Threshold -doc: Threshold on the convergence of the MRCC energy -interface: ezfio,provider,ocaml -default: 1.e-5 - -[n_it_mrcc_max] -type: Strictly_positive_int -doc: Maximum number of MRCC iterations -interface: ezfio,provider,ocaml -default: 10 - diff --git a/plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES b/plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES deleted file mode 100644 index a8404d62..00000000 --- a/plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils diff --git a/plugins/MRCC_CASSD/README.rst b/plugins/MRCC_CASSD/README.rst deleted file mode 100644 index 0e38fe4c..00000000 --- a/plugins/MRCC_CASSD/README.rst +++ /dev/null @@ -1,80 +0,0 @@ -=========== -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 - -Needed Modules -============== -.. Do not edit this section It was auto-generated -.. 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 -.. by the `update_README.py` script. - - -`mrcc `_ - Undocumented - - -`mrcc_noiter `_ - Undocumented - - -`n_it_mrcc_max `_ - Maximum number of MRCC iterations - - -`print_cas_coefs `_ - Undocumented - - -`run `_ - Undocumented - - -`run_pt2 `_ - Undocumented - - -`thresh_mrcc `_ - Threshold on the convergence of the MRCC energy - diff --git a/plugins/MRCC_CASSD/mrcc_cassd.irp.f b/plugins/MRCC_CASSD/mrcc_cassd.irp.f deleted file mode 100644 index 4fef815d..00000000 --- a/plugins/MRCC_CASSD/mrcc_cassd.irp.f +++ /dev/null @@ -1,120 +0,0 @@ -program mrcc - implicit none - double precision, allocatable :: energy(:) - allocate (energy(N_states)) - - read_wf = .True. - SOFT_TOUCH read_wf - call print_cas_coefs - call set_generators_bitmasks_as_holes_and_particles - call run(N_states,energy) - if(do_pt2_end)then - call run_pt2(N_states,energy) - endif - deallocate(energy) -end - -subroutine run(N_st,energy) - implicit none - - integer, intent(in) :: N_st - double precision, intent(out) :: energy(N_st) - - integer :: i - - double precision :: E_new, E_old, delta_e - integer :: iteration - double precision :: E_past(4), lambda - E_new = 0.d0 - delta_E = 1.d0 - iteration = 0 - lambda = 1.d0 - do while (delta_E > thresh_mrcc) - iteration += 1 - print *, '===========================' - print *, 'MRCC Iteration', iteration - print *, '===========================' - print *, '' - E_old = sum(ci_energy_dressed(1:N_st)) - call write_double(6,ci_energy_dressed(1),"MRCC energy") - call diagonalize_ci_dressed(lambda) - E_new = sum(ci_energy_dressed(1:N_st)) - delta_E = dabs(E_new - E_old) - call save_wavefunction - call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) - if (iteration > n_it_mrcc_max) then - exit - endif - enddo - call write_double(6,ci_energy_dressed(1),"Final MRCC energy") - energy(1:N_st) = ci_energy_dressed(1:N_st) - -end - - -subroutine run_pt2(N_st,energy) - implicit none - integer :: i,j,k - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) - integer, intent(in) :: N_st - double precision, intent(in) :: energy(N_st) - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) - pt2 = 0.d0 - - print*,'Last iteration only to compute the PT2' - threshold_selectors = 1.d0 - threshold_generators = 0.999d0 - - N_det_generators = lambda_mrcc_pt2(0) + N_det_cas - do i=1,N_det_cas - do k=1,N_int - psi_det_generators(k,1,i) = psi_ref(k,1,i) - psi_det_generators(k,2,i) = psi_ref(k,2,i) - enddo - do k=1,N_st - psi_coef_generators(i,k) = psi_ref_coef(i,k) - enddo - enddo - do i=N_det_cas+1,N_det_generators - j = lambda_mrcc_pt2(i) - do k=1,N_int - psi_det_generators(k,1,i) = psi_non_ref(k,1,j) - psi_det_generators(k,2,i) = psi_non_ref(k,2,j) - enddo - do k=1,N_st - psi_coef_generators(i,k) = psi_non_ref_coef(j,k) - enddo - enddo - SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed - - - call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st) - print *, 'Final step' - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', energy - print *, 'E+PT2 = ', energy+pt2 - print *, '-----' - - - call ezfio_set_full_ci_energy_pt2(energy+pt2) - deallocate(pt2,norm_pert) - -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_CASSD/mrcc_noiter.irp.f b/plugins/MRCC_CASSD/mrcc_noiter.irp.f deleted file mode 100644 index 8d95cea9..00000000 --- a/plugins/MRCC_CASSD/mrcc_noiter.irp.f +++ /dev/null @@ -1,91 +0,0 @@ -program mrcc_noiter - implicit none - double precision, allocatable :: energy(:) - allocate (energy(N_states)) - read_wf = .True. - threshold_generators = .9999d0 - SOFT_TOUCH read_wf threshold_generators - call print_cas_coefs - call set_generators_bitmasks_as_holes_and_particles - call run(N_states,energy) - if(do_pt2_end)then - call run_pt2(N_states,energy) - endif - deallocate(energy) -end - -subroutine run(N_st,energy) - implicit none - - integer, intent(in) :: N_st - double precision, intent(out) :: energy(N_st) - integer :: i,j - do j=1,N_states_diag - do i=1,N_det - psi_coef(i,j) = CI_eigenvectors_dressed(i,j) - enddo - enddo - SOFT_TOUCH psi_coef ci_energy_dressed - call write_double(6,ci_energy_dressed(1),"Final MRCC energy") - call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) - call save_wavefunction - energy(:) = ci_energy_dressed(:) -end - - -subroutine run_pt2(N_st,energy) - implicit none - integer :: i,j,k - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) - integer, intent(in) :: N_st - double precision, intent(in) :: energy(N_st) - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) - pt2 = 0.d0 - - print*,'Last iteration only to compute the PT2' - threshold_selectors = 1.d0 - threshold_generators = 0.999d0 - - N_det_generators = lambda_mrcc_pt2(0) - do i=1,N_det_generators - j = lambda_mrcc_pt2(i) - do k=1,N_int - psi_det_generators(k,1,i) = psi_non_ref(k,1,j) - psi_det_generators(k,2,i) = psi_non_ref(k,2,j) - enddo - do k=1,N_st - psi_coef_generators(i,k) = psi_non_ref_coef(j,k) - enddo - enddo - SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed - - - call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st) - print *, 'Final step' - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', energy - print *, 'E+PT2 = ', energy+pt2 - print *, '-----' - - call ezfio_set_full_ci_energy_pt2(energy+pt2) - deallocate(pt2,norm_pert) - -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_CASSD/tree_dependency.png b/plugins/MRCC_CASSD/tree_dependency.png deleted file mode 100644 index e73ff165..00000000 Binary files a/plugins/MRCC_CASSD/tree_dependency.png and /dev/null differ diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 7687b451..79249051 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] integer :: j character*(8) :: st call write_time(output_determinants) - do j=1,min(N_det,N_states_diag) + do j=1,min(N_det,N_states) write(st,'(I4)') j CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion call write_double(output_determinants,CI_energy_dressed(j),'Energy of state '//trim(st)) @@ -349,8 +349,9 @@ integer function searchDet(dets, det, n, Nint) do while(.true.) searchDet = (l+h)/2 c = detCmp(dets(1,1,searchDet), det(1,1), Nint) - if(c == 0) return - if(c == 1) then + if(c == 0) then + return + else if(c == 1) then h = searchDet-1 else l = searchDet+1 @@ -498,12 +499,12 @@ subroutine tamise_exc(key, no, n, N_key) BEGIN_DOC ! Uncodumented : TODO END_DOC - integer,intent(in) :: no, n, N_key + integer,intent(in) :: no, n, N_key integer*2,intent(inout) :: key(4, N_key) - integer :: k,j - integer*2 :: tmp(4) - logical :: exc_inf - integer :: ni + integer :: k,j + integer*2 :: tmp(4) + logical :: exc_inf + integer :: ni k = no j = 2*k @@ -633,7 +634,6 @@ END_PROVIDER allocate(A_ind(0:N_det_ref+1, nex), A_val(N_det_ref+1, nex)) allocate(AtA_ind(N_det_ref * nex), AtA_val(N_det_ref * nex)) !!!!! MAY BE TOO SMALL ? !!!!!!!! allocate(x(nex), AtB(nex)) - allocate(A_val_mwen(nex), A_ind_mwen(nex)) allocate(N_col(nex), col_shortcut(nex)) allocate(x_new(nex)) @@ -645,7 +645,6 @@ END_PROVIDER AtA_val = 0d0 x = 0d0 A_val_mwen = 0d0 - A_ind_mwen = 0 N_col = 0 col_shortcut = 0 @@ -689,7 +688,8 @@ END_PROVIDER end do !$OMP END DO deallocate(lref) - !$OMP END PARALLEL + !$OMP END PARALLEL + print *, 'Done building A_val, A_ind' AtB = 0d0 AtA_size = 0 @@ -698,6 +698,8 @@ END_PROVIDER !$OMP PARALLEL default(none) shared(k, psi_non_ref_coef, A_ind, A_val, x, N_det_ref, nex, N_det_non_ref)& !$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen)& !$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind, s) + allocate(A_val_mwen(nex), A_ind_mwen(nex)) + A_ind_mwen = 0 !$OMP DO schedule(dynamic, 100) do at_row = 1, nex wk = 0 @@ -711,10 +713,10 @@ END_PROVIDER r1 = 1 r2 = 1 do while ((A_ind(r1, at_row) /= 0).and.(A_ind(r2, a_col) /= 0)) - if(A_ind(r1, at_row) < A_ind(r2, a_col)) then - r1 = r1+1 - else if(A_ind(r1, at_row) > A_ind(r2, a_col)) then + if(A_ind(r1, at_row) > A_ind(r2, a_col)) then r2 = r2+1 + else if(A_ind(r1, at_row) < A_ind(r2, a_col)) then + r1 = r1+1 else t = t - A_val(r1, at_row) * A_val(r2, a_col) r1 = r1+1 @@ -748,7 +750,8 @@ END_PROVIDER !$OMP END CRITICAL end if end do - !$OMP END DO + !$OMP END DO NOWAIT + deallocate (A_ind_mwen, A_val_mwen) !$OMP END PARALLEL if(AtA_size > size(AtA_val)) stop "SIZA" @@ -807,14 +810,18 @@ END_PROVIDER if(res < 1d-8) exit end do - + ! rho_mrcc now contains A.X + norm = 0.d0 do i=1,N_det_non_ref norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s) enddo + ! Norm now contains the norm of A.X + do i=1,N_det_ref norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s) enddo + ! Norm now contains the norm of Psi + A.X print *, k, "res : ", res, "norm : ", sqrt(norm) @@ -826,30 +833,46 @@ END_PROVIDER if (rho_mrcc(i,s) == 0.d0) then rho_mrcc(i,s) = 1.d-32 endif + + ! f is such that f.\tilde{c_i} = c_i f = psi_non_ref_coef(i,s) / rho_mrcc(i,s) + ! Avoid numerical instabilities f = min(f,2.d0) f = max(f,-2.d0) + norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s) rho_mrcc(i,s) = f enddo + ! norm now contains the norm of |T.Psi_0> + ! rho_mrcc now contains the f factors f = 1.d0/norm - norm = 0.d0 - do i=1,N_det_non_ref - norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s) - enddo - f = dsqrt(f*norm) + ! f now contains 1/ - print *, 'norm of |T Psi_0> = ', norm*f*f + norm = 1.d0 + do i=1,N_det_ref + norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s) + enddo + ! norm now contains + f = dsqrt(f*norm) + ! f normalises T.Psi_0 such that (1+T)|Psi> is normalized + + norm = norm*f + print *, 'norm of |T Psi_0> = ', dsqrt(norm) + + do i=1,N_det_ref + norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s) + enddo do i=1,N_det_non_ref rho_mrcc(i,s) = rho_mrcc(i,s) * f enddo + ! rho_mrcc now contains the product of the scaling factors and the + ! normalization constant end do - print *, "done" END_PROVIDER @@ -857,13 +880,17 @@ BEGIN_PROVIDER [ double precision, dij, (N_det_ref, N_det_non_ref, N_states) ] integer :: s,i,j double precision, external :: get_dij_index print *, "computing amplitudes..." + !$OMP PARALLEL DEFAULT(shared) PRIVATE(s,i,j) do s=1, N_states + !$OMP DO do i=1, N_det_non_ref do j=1, N_det_ref dij(j, i, s) = get_dij_index(j, i, s, N_int) end do end do + !$OMP END DO end do + !$OMP END PARALLEL print *, "done computing amplitudes" END_PROVIDER @@ -876,10 +903,9 @@ double precision function get_dij_index(II, i, s, Nint) double precision :: HIi, phase if(lambda_type == 0) then - get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) -! get_dij_index = get_dij_index * rho_mrcc(i,s) call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int) - get_dij_index = get_dij_index * rho_mrcc(i,s) * phase + get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase + get_dij_index = get_dij_index * rho_mrcc(i,s) else call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) get_dij_index = HIi * lambda_mrcc(s, i) @@ -1059,18 +1085,20 @@ subroutine apply_hole_local(det, exc, res, ok, Nint) res = det if(h1 /= 0) then - ii = (h1-1)/bit_kind_size + 1 - pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64 - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return - res(ii, s1) = ibclr(res(ii, s1), pos) + ii = (h1-1)/bit_kind_size + 1 + pos = iand(h1-1,bit_kind_size-1) ! mod 64 + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) then + return + endif + res(ii, s1) = ibclr(res(ii, s1), pos) end if - ii = (h2-1)/bit_kind_size + 1 - pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1) - if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return - res(ii, s2) = ibclr(res(ii, s2), pos) - - + ii = (h2-1)/bit_kind_size + 1 + pos = iand(h2-1,bit_kind_size-1) ! mod 64 + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) then + return + endif + res(ii, s2) = ibclr(res(ii, s2), pos) ok = .true. end subroutine @@ -1094,16 +1122,20 @@ subroutine apply_particle_local(det, exc, res, ok, Nint) res = det if(p1 /= 0) then - ii = (p1-1)/bit_kind_size + 1 - pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return - res(ii, s1) = ibset(res(ii, s1), pos) + ii = (p1-1)/bit_kind_size + 1 + pos = iand(p1-1,bit_kind_size-1) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) then + return + endif + res(ii, s1) = ibset(res(ii, s1), pos) end if - ii = (p2-1)/bit_kind_size + 1 - pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) - if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return - res(ii, s2) = ibset(res(ii, s2), pos) + ii = (p2-1)/bit_kind_size + 1 + pos = iand(p2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) then + return + endif + res(ii, s2) = ibset(res(ii, s2), pos) ok = .true. diff --git a/plugins/Psiref_CAS/overwrite_with_cas.irp.f b/plugins/Psiref_CAS/overwrite_with_cas.irp.f index 4d3d217d..d3ced1d1 100644 --- a/plugins/Psiref_CAS/overwrite_with_cas.irp.f +++ b/plugins/Psiref_CAS/overwrite_with_cas.irp.f @@ -1,3 +1,5 @@ program overwrite_w_cas + read_wf = .True. + TOUCH read_wf call extract_ref end diff --git a/plugins/QmcChem/dressed_dmc.irp.f b/plugins/QmcChem/dressed_dmc.irp.f index 1437ed31..803e55dc 100644 --- a/plugins/QmcChem/dressed_dmc.irp.f +++ b/plugins/QmcChem/dressed_dmc.irp.f @@ -57,7 +57,7 @@ program dressed_dmc enddo - call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,,N_int,6) + call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,N_int,6) call save_wavefunction call write_spindeterminants diff --git a/plugins/Selectors_full/zmq.irp.f b/plugins/Selectors_full/zmq.irp.f index dfa94884..8046212b 100644 --- a/plugins/Selectors_full/zmq.irp.f +++ b/plugins/Selectors_full/zmq.irp.f @@ -101,6 +101,7 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id, energy, size_energy) print *, '77_zmq_recv(zmq_to_qp_run_socket,psi_coef,psi_det_size*N_states*8,ZMQ_SNDMORE)' stop 'error' endif + TOUCH psi_det psi_coef rc = f77_zmq_recv(zmq_to_qp_run_socket,energy,size_energy*8,0) if (rc /= size_energy*8) then @@ -110,11 +111,12 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id, energy, size_energy) if (N_det_generators_read > 0) then N_det_generators = N_det_generators_read + TOUCH N_det_generators endif if (N_det_selectors_read > 0) then N_det_selectors = N_det_selectors_read + TOUCH N_det_selectors endif - SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators end diff --git a/plugins/mrcepa0/EZFIO.cfg b/plugins/mrcepa0/EZFIO.cfg index 7580f028..d792390d 100644 --- a/plugins/mrcepa0/EZFIO.cfg +++ b/plugins/mrcepa0/EZFIO.cfg @@ -3,3 +3,31 @@ type: Positive_int doc: lambda type interface: ezfio,provider,ocaml default: 0 + +[energy] +type: double precision +doc: Calculated energy +interface: ezfio + +[energy_pt2] +type: double precision +doc: Calculated energy with PT2 contribution +interface: ezfio + +[energy] +type: double precision +doc: Calculated energy +interface: ezfio + +[thresh_dressed_ci] +type: Threshold +doc: Threshold on the convergence of the dressed CI energy +interface: ezfio,provider,ocaml +default: 1.e-4 + +[n_it_max_dressed_ci] +type: Strictly_positive_int +doc: Maximum number of dressed CI iterations +interface: ezfio,provider,ocaml +default: 10 + diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index b82bc613..1e7ad68d 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -17,8 +17,8 @@ subroutine run(N_st,energy) - thresh_mrcc = 1d-7 - n_it_mrcc_max = 10 + thresh_mrcc = thresh_dressed_ci + n_it_mrcc_max = n_it_max_dressed_ci if(n_it_mrcc_max == 1) then do j=1,N_states_diag @@ -48,8 +48,8 @@ subroutine run(N_st,energy) E_new = sum(ci_energy_dressed) delta_E = dabs(E_new - E_old) call save_wavefunction - call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) - if (iteration > n_it_mrcc_max) then + call ezfio_set_mrcepa0_energy(ci_energy_dressed(1)) + if (iteration >= n_it_mrcc_max) then exit endif enddo @@ -184,7 +184,7 @@ subroutine run_pt2_old(N_st,energy) print *, '-----' -! call ezfio_set_full_ci_energy_pt2(energy+pt2) + call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1)) end @@ -238,5 +238,7 @@ subroutine run_pt2(N_st,energy) print *, 'E+PT2 = ', energy+pt2 print *, '-----' + call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1)) + end diff --git a/src/Davidson/davidson_parallel.irp.f b/src/Davidson/davidson_parallel.irp.f index 0e91ac90..cede52c9 100644 --- a/src/Davidson/davidson_parallel.irp.f +++ b/src/Davidson/davidson_parallel.irp.f @@ -4,12 +4,12 @@ use bitmasks use f77_zmq -subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs) +subroutine davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep) implicit none - integer , intent(in) :: blockb, blocke, bs + integer , intent(in) :: blockb, bs, blockb2, istep integer , intent(inout) :: N integer , intent(inout) :: idx(bs) double precision , intent(inout) :: vt(N_states_diag, bs) @@ -20,11 +20,12 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs) double precision :: s2, hij logical, allocatable :: wrotten(:) - allocate(wrotten(bs)) wrotten = .false. + PROVIDE dav_det - do sh = blockb, blocke + ii=0 + sh = blockb do sh2=1,shortcut_(0,1) exa = 0 do ni=1,N_int @@ -32,7 +33,7 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs) end do if(exa > 2) cycle - do i=shortcut_(sh,1),shortcut_(sh+1,1)-1 + do i=blockb2+shortcut_(sh,1),shortcut_(sh+1,1)-1, istep ii = i - shortcut_(blockb,1) + 1 org_i = sort_idx_(i,1) @@ -48,8 +49,8 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs) ext = ext + popcnt(xor(sorted_i(ni), sorted_(ni,j,1))) end do if(ext <= 4) then - call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij) call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2) + call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij) if(.not. wrotten(ii)) then wrotten(ii) = .true. idx(ii) = org_i @@ -64,39 +65,39 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs) enddo enddo enddo - enddo - do sh=blockb,min(blocke, shortcut_(0,2)) - do sh2=sh, shortcut_(0,2), shortcut_(0,1)*51 - do i=shortcut_(sh2,2),shortcut_(sh2+1,2)-1 - ii += 1 - org_i = sort_idx_(i,2) - do j=shortcut_(sh2,2),shortcut_(sh2+1,2)-1 - if(i == j) cycle - org_j = sort_idx_(j,2) - ext = 0 - do ni=1,N_int - ext = ext + popcnt(xor(sorted_(ni,i,2), sorted_(ni,j,2))) + if (blockb <= shortcut_(0,2)) then + sh=blockb + do sh2=sh, shortcut_(0,2), shortcut_(0,1) + do i=blockb2+shortcut_(sh2,2),shortcut_(sh2+1,2)-1, istep + ii += 1 + org_i = sort_idx_(i,2) + do j=shortcut_(sh2,2),shortcut_(sh2+1,2)-1 + if(i == j) cycle + org_j = sort_idx_(j,2) + ext = 0 + do ni=1,N_int + ext = ext + popcnt(xor(sorted_(ni,i,2), sorted_(ni,j,2))) + end do + if(ext == 4) then + call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij) + call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2) + if(.not. wrotten(ii)) then + wrotten(ii) = .true. + idx(ii) = org_i + vt (:,ii) = 0d0 + st (:,ii) = 0d0 + end if + do istate=1,N_states_diag + vt (istate,ii) += hij*dav_ut(istate,org_j) + st (istate,ii) += s2*dav_ut(istate,org_j) + enddo + end if end do - if(ext == 4) then - call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij) - call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2) - if(.not. wrotten(ii)) then - wrotten(ii) = .true. - idx(ii) = org_i - vt (:,ii) = 0d0 - st (:,ii) = 0d0 - end if - do istate=1,N_states_diag - vt (istate,ii) += hij*dav_ut(istate,org_j) - st (istate,ii) += s2*dav_ut(istate,org_j) - enddo - end if end do - end do - enddo - enddo + enddo + endif N=0 do i=1,bs @@ -107,16 +108,17 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs) st(:,N) = st(:,i) end if end do + + end subroutine -subroutine davidson_collect(blockb, blocke, N, idx, vt, st , v0t, s0t) +subroutine davidson_collect(N, idx, vt, st , v0t, s0t) implicit none - integer , intent(in) :: blockb, blocke integer , intent(in) :: N integer , intent(in) :: idx(N) double precision , intent(in) :: vt(N_states_diag, N) @@ -138,28 +140,48 @@ subroutine davidson_collect(blockb, blocke, N, idx, vt, st , v0t, s0t) end subroutine -subroutine davidson_init(zmq_to_qp_run_socket) +subroutine davidson_init(zmq_to_qp_run_socket,n,n_st_8,ut) use f77_zmq implicit none - integer(ZMQ_PTR), intent(out) :: zmq_to_qp_run_socket + integer(ZMQ_PTR), intent(out) :: zmq_to_qp_run_socket + integer, intent(in) :: n, n_st_8 + double precision, intent(in) :: ut(n_st_8,n) + integer :: i,k + - touch dav_size dav_det dav_ut + dav_size = n + touch dav_size + + do i=1,n + do k=1,N_int + dav_det(k,1,i) = psi_det(k,1,i) + dav_det(k,2,i) = psi_det(k,2,i) + enddo + enddo + do i=1,n + do k=1,N_states_diag + dav_ut(k,i) = ut(k,i) + enddo + enddo + + touch dav_det dav_ut + call new_parallel_job(zmq_to_qp_run_socket,"davidson") end subroutine -subroutine davidson_add_task(zmq_to_qp_run_socket, blockb, blocke) +subroutine davidson_add_task(zmq_to_qp_run_socket, blockb, blockb2, istep) use f77_zmq implicit none integer(ZMQ_PTR) ,intent(in) :: zmq_to_qp_run_socket - integer ,intent(in) :: blockb, blocke + integer ,intent(in) :: blockb, blockb2, istep character*(512) :: task - write(task,*) blockb, blocke + write(task,*) blockb, blockb2, istep call add_task_to_taskserver(zmq_to_qp_run_socket, task) end subroutine @@ -188,7 +210,7 @@ subroutine davidson_run_slave(thread,iproc) integer, intent(in) :: thread, iproc - integer :: worker_id, task_id, blockb, blocke + integer :: worker_id, task_id, blockb character*(512) :: task integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket @@ -228,7 +250,7 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id) character*(512) :: task - integer :: blockb, blocke + integer :: blockb, blockb2, istep integer :: N integer , allocatable :: idx(:) double precision , allocatable :: vt(:,:) @@ -241,10 +263,10 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id) do call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) if(task_id == 0) exit - read (task,*) blockb, blocke - bs = shortcut_(blocke+1,1) - shortcut_(blockb, 1) + read (task,*) blockb, blockb2, istep + bs = shortcut_(blockb+1,1) - shortcut_(blockb, 1) do i=blockb, shortcut_(0,2), shortcut_(0,1) - do j=i, min(i+blocke-blockb, shortcut_(0,2)) + do j=i, min(i, shortcut_(0,2)) bs += shortcut_(j+1,2) - shortcut_(j, 2) end do end do @@ -255,11 +277,11 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id) allocate(st(N_states_diag, bs)) end if - call davidson_process(blockb, blocke, N, idx, vt, st, bs) - + call davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep) call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) - call davidson_push_results(zmq_socket_push, blockb, blocke, N, idx, vt, st, task_id) + call davidson_push_results(zmq_socket_push, blockb, blockb2, N, idx, vt, st, task_id) end do + end subroutine @@ -340,28 +362,28 @@ end subroutine -subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0) +subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0, LDA) use f77_zmq implicit none - integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket - integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - - double precision ,intent(inout) :: v0(dav_size, N_states_diag) - double precision ,intent(inout) :: s0(dav_size, N_states_diag) - - integer :: more, task_id + integer :: LDA + integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull - - integer :: blockb, blocke - integer :: N - integer , allocatable :: idx(:) - double precision , allocatable :: vt(:,:), v0t(:,:), s0t(:,:) - double precision , allocatable :: st(:,:) + double precision ,intent(inout) :: v0(LDA, N_states_diag) + double precision ,intent(inout) :: s0(LDA, N_states_diag) + + integer :: more, task_id, taskn + + integer :: blockb, blocke + integer :: N + integer , allocatable :: idx(:) + double precision , allocatable :: vt(:,:), v0t(:,:), s0t(:,:) + double precision , allocatable :: st(:,:) integer :: msize - msize = (max_workload + max_blocksize)*2 + msize = (1 + max_blocksize)*2 allocate(idx(msize)) allocate(vt(N_states_diag, msize)) allocate(st(N_states_diag, msize)) @@ -376,32 +398,41 @@ subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0) do while (more == 1) call davidson_pull_results(zmq_socket_pull, blockb, blocke, N, idx, vt, st, task_id) !DIR$ FORCEINLINE - call davidson_collect(blockb, blocke, N, idx, vt, st , v0t, s0t) + call davidson_collect(N, idx, vt, st , v0t, s0t) call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) end do deallocate(idx,vt,st) - call dtranspose(v0t,size(v0t,1), v0, size(v0,1), N_states_diag, dav_size) - call dtranspose(s0t,size(s0t,1), s0, size(s0,1), N_states_diag, dav_size) + integer :: i,j + !DIR$ IVDEP + do j=1,N_states_diag + !DIR$ IVDEP + do i=1,dav_size + v0(i,j) = v0t(j,i) + s0(i,j) = s0t(j,i) + enddo + enddo + deallocate(v0t,s0t) end subroutine -subroutine davidson_run(zmq_to_qp_run_socket , v0, s0) +subroutine davidson_run(zmq_to_qp_run_socket , v0, s0, LDA) use f77_zmq implicit none - integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket + integer :: LDA + integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_collector + integer(ZMQ_PTR) :: zmq_collector integer(ZMQ_PTR), external :: new_zmq_pull_socket integer(ZMQ_PTR) :: zmq_socket_pull integer :: i integer, external :: omp_get_thread_num - double precision , intent(inout) :: v0(dav_size, N_states_diag) - double precision , intent(inout) :: s0(dav_size, N_states_diag) + double precision , intent(inout) :: v0(LDA, N_states_diag) + double precision , intent(inout) :: s0(LDA, N_states_diag) call zmq_set_running(zmq_to_qp_run_socket) @@ -412,19 +443,20 @@ subroutine davidson_run(zmq_to_qp_run_socket , v0, s0) PROVIDE nproc - !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+2) - i = omp_get_thread_num() - if (i==0) then - call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0) - call end_zmq_to_qp_run_socket(zmq_collector) - call end_zmq_pull_socket(zmq_socket_pull) - call davidson_miniserver_end() - else if(i==1) then - call davidson_miniserver_run() - else - call davidson_slave_inproc(i) - endif + !$OMP PARALLEL NUM_THREADS(nproc+2) PRIVATE(i) + i = omp_get_thread_num() + if (i == 0 ) then + call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0, LDA) + call end_zmq_to_qp_run_socket(zmq_collector) + call end_zmq_pull_socket(zmq_socket_pull) + call davidson_miniserver_end() + else if (i == 1 ) then + call davidson_miniserver_run () + else + call davidson_slave_inproc(i) + endif !$OMP END PARALLEL + call end_parallel_job(zmq_to_qp_run_socket, 'davidson') end subroutine @@ -505,15 +537,28 @@ end subroutine -BEGIN_PROVIDER [ integer(bit_kind), dav_det, (N_int, 2, dav_size) ] -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ] + BEGIN_PROVIDER [ integer(bit_kind), dav_det, (N_int, 2, dav_size) ] +&BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ] + use bitmasks + implicit none + BEGIN_DOC +! Temporary arrays for parallel davidson +! +! Touched in davidson_miniserver_get + END_DOC + dav_det = 0_bit_kind + dav_ut = -huge(1.d0) END_PROVIDER BEGIN_PROVIDER [ integer, dav_size ] + implicit none + BEGIN_DOC +! Size of the arrays for Davidson +! +! Touched in davidson_miniserver_get + END_DOC + dav_size = 1 END_PROVIDER diff --git a/src/Davidson/davidson_slave.irp.f b/src/Davidson/davidson_slave.irp.f index 9195e46f..e28712e2 100644 --- a/src/Davidson/davidson_slave.irp.f +++ b/src/Davidson/davidson_slave.irp.f @@ -26,10 +26,10 @@ program davidson_slave print *, 'Davidson slave running' - ! !$OMP PARALLEL PRIVATE(i) - !i = omp_get_thread_num() - call davidson_slave_tcp(0) - !!$OMP END PARALLEL + !$OMP PARALLEL PRIVATE(i) + i = omp_get_thread_num() + call davidson_slave_tcp(i) + !$OMP END PARALLEL end do end diff --git a/src/Davidson/diagonalization_hs2.irp.f b/src/Davidson/diagonalization_hs2.irp.f index c44a27d2..c8ac3733 100644 --- a/src/Davidson/diagonalization_hs2.irp.f +++ b/src/Davidson/diagonalization_hs2.irp.f @@ -1,4 +1,4 @@ -subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit) +subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit) use bitmasks implicit none BEGIN_DOC @@ -22,7 +22,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nin integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(inout) :: u_in(dim_in,N_st_diag) - double precision, intent(out) :: energies(N_st) + double precision, intent(out) :: energies(N_st), s2_out(N_st) double precision, allocatable :: H_jj(:), S2_jj(:) double precision :: diag_h_mat_elem @@ -46,6 +46,9 @@ subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nin !$OMP END PARALLEL call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit) + do i=1,N_st_diag + s2_out(i) = S2_jj(i) + enddo deallocate (H_jj,S2_jj) end @@ -79,7 +82,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s END_DOC integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) - double precision, intent(in) :: H_jj(sze), S2_jj(sze) + double precision, intent(in) :: H_jj(sze) + double precision, intent(inout) :: S2_jj(sze) integer, intent(in) :: iunit double precision, intent(inout) :: u_in(dim_in,N_st_diag) double precision, intent(out) :: energies(N_st_diag) @@ -200,7 +204,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s call normalize(u_in(1,k),sze) enddo - + do while (.not.converged) do k=1,N_st_diag @@ -239,12 +243,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! enddo ! enddo - call dgemm('T','N', shift2, N_st_diag, sze, & - 1.d0, U, size(U,1), W(1,shift+1), size(W,1), & + call dgemm('T','N', shift2, N_st_diag, sze, & + 1.d0, U, size(U,1), W(1,shift+1), size(W,1), & 0.d0, h(1,shift+1), size(h,1)) - call dgemm('T','N', shift2, N_st_diag, sze, & - 1.d0, U, size(U,1), S(1,shift+1), size(S,1), & + call dgemm('T','N', shift2, N_st_diag, sze, & + 1.d0, U, size(U,1), S(1,shift+1), size(S,1), & 0.d0, s_(1,shift+1), size(s_,1)) ! Diagonalize h @@ -254,14 +258,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! Compute S2 for each eigenvector ! ------------------------------- - call dgemm('N','N',shift2,shift2,shift2, & - 1.d0, s_, size(s_,1), y, size(y,1), & + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, s_, size(s_,1), y, size(y,1), & 0.d0, s_tmp, size(s_tmp,1)) - - call dgemm('T','N',shift2,shift2,shift2, & - 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & 0.d0, s_, size(s_,1)) - + do k=1,shift2 s2(k) = s_(k,k) + S_z2_Sz enddo @@ -326,8 +330,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! enddo do k=1,N_st_diag do i=1,sze - R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & - * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz) + R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz) enddo if (k <= N_st) then residual_norm(k) = u_dot_u(R(1,k),sze) @@ -367,10 +371,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! enddo ! enddo ! - call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), & - U(1,shift2+k),1,0.d0,c,1) - call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), & - c,1,1.d0,U(1,shift2+k),1) + call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), & + U(1,shift2+k),1,0.d0,c,1) + call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), & + c,1,1.d0,U(1,shift2+k),1) ! ! do l=1,k-1 ! c(1) = u_dot_v(U(1,shift2+k),U(1,shift2+l),sze) @@ -379,7 +383,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! enddo ! enddo ! - call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), & + call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), & U(1,shift2+k),1,0.d0,c,1) call dgemv('N',sze,k-1,-1.d0,U(1,shift2+1),size(U,1), & c,1,1.d0,U(1,shift2+k),1) diff --git a/src/Davidson/diagonalize_CI.irp.f b/src/Davidson/diagonalize_CI.irp.f index ecd2d6b2..3b2c9ed0 100644 --- a/src/Davidson/diagonalize_CI.irp.f +++ b/src/Davidson/diagonalize_CI.irp.f @@ -59,13 +59,13 @@ END_PROVIDER ! size(CI_eigenvectors,1), & ! N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants) ! - call davidson_diag_HS2(psi_det,CI_eigenvectors, & +! call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,& +! min(N_det,N_states_diag),size(CI_eigenvectors,1)) + + call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_eigenvectors_s2, & size(CI_eigenvectors,1),CI_electronic_energy, & N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants) - call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,& - min(N_det,N_states_diag),size(CI_eigenvectors,1)) - else if (diag_algorithm == "Lapack") then diff --git a/src/Davidson/diagonalize_CI_SC2.irp.f b/src/Davidson/diagonalize_CI_SC2.irp.f deleted file mode 100644 index 498792d9..00000000 --- a/src/Davidson/diagonalize_CI_SC2.irp.f +++ /dev/null @@ -1,62 +0,0 @@ -BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states_diag) ] - implicit none - BEGIN_DOC - ! N_states_diag lowest eigenvalues of the CI matrix - END_DOC - - integer :: j - character*(8) :: st - call write_time(output_determinants) - do j=1,N_states_diag - CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion - write(st,'(I4)') j - call write_double(output_determinants,CI_SC2_energy(j),'Energy of state '//trim(st)) - enddo - -END_PROVIDER - - BEGIN_PROVIDER [ double precision, threshold_convergence_SC2] - implicit none - BEGIN_DOC - ! convergence of the correlation energy of SC2 iterations - END_DOC - threshold_convergence_SC2 = 1.d-10 - - END_PROVIDER - - BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ] -&BEGIN_PROVIDER [ double precision, Diag_H_elements_SC2, (N_det) ] - implicit none - BEGIN_DOC - ! Eigenvectors/values of the CI matrix - END_DOC - integer :: i,j - - do j=1,N_states_diag - do i=1,N_det - CI_SC2_eigenvectors(i,j) = psi_coef(i,j) - enddo - CI_SC2_electronic_energy(j) = CI_electronic_energy(j) - enddo - - call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & -! size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) - diag_H_elements_SC2,size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) -END_PROVIDER - -subroutine diagonalize_CI_SC2 - implicit none - BEGIN_DOC -! Replace the coefficients of the CI states_diag by the coefficients of the -! eigenstates of the CI matrix - END_DOC - integer :: i,j - do j=1,N_states_diag - do i=1,N_det - psi_coef(i,j) = CI_SC2_eigenvectors(i,j) - enddo - enddo - SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors diag_h_elements_sc2 -! SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors -end diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index 475487d0..a1a72100 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -24,6 +24,7 @@ subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8) do i=1,N_st e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n) enddo + deallocate (H_jj, v_0) end @@ -56,10 +57,9 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) integer :: N_st_8 integer, external :: align_double - !!!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut - if(N_st /= N_states_diag) stop "H_u_0_nstates N_st /= N_states_diag" - N_st_8 = N_states_diag ! align_double(N_st) + N_st_8 = align_double(N_st) ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -165,7 +165,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) v_0(i,istate) += H_jj(i) * u_0(i,istate) enddo enddo - !deallocate (shortcut, sort_idx, sorted, version, ut) + deallocate (shortcut, sort_idx, sorted, version, ut) end BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ] @@ -200,76 +200,82 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) integer :: i,j,k,l, jj,ii integer :: i0, j0 - integer, allocatable :: shortcut(:,:), sort_idx(:,:) - integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) + integer, allocatable :: shortcut(:,:), sort_idx(:) + integer(bit_kind), allocatable :: sorted(:,:), version(:,:) integer(bit_kind) :: sorted_i(Nint) integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate integer :: N_st_8 integer, external :: align_double - integer :: workload, blockb, blocke + integer :: blockb, blockb2, istep + double precision :: ave_workload, workload integer(ZMQ_PTR) :: handler if(N_st /= N_states_diag .or. sze_8 < N_det) stop "assert fail in H_S2_u_0_nstates" - N_st_8 = N_st !! align_double(N_st) + N_st_8 = N_st ! align_double(N_st) ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) PROVIDE ref_bitmask_energy - allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) + allocate (shortcut(0:n+1,2), sort_idx(n), sorted(Nint,n), version(Nint,n)) allocate(ut(N_st_8,n)) v_0 = 0.d0 s_0 = 0.d0 - if(n /= N_det) stop "n /= N_det" - do i=1,n do istate=1,N_st ut(istate,i) = u_0(i,istate) enddo enddo - call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) - call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) + call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut(0,1), version, n, Nint) + call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut(0,2), version, n, Nint) - dav_size = n - touch dav_size - dav_det = psi_det - dav_ut = ut - - workload = 0 blockb = shortcut(0,1) - blocke = blockb - call davidson_init(handler) + call davidson_init(handler,n,N_st_8,ut) + + + ave_workload = 0.d0 + do sh=1,shortcut(0,1) + ave_workload += shortcut(0,1) + ave_workload += (shortcut(sh+1,1) - shortcut(sh,1))**2 + do i=sh, shortcut(0,2), shortcut(0,1) + do j=i, min(i, shortcut(0,2)) + ave_workload += (shortcut(j+1,2) - shortcut(j, 2))**2 + end do + end do + enddo + ave_workload = ave_workload/dble(shortcut(0,1)) + + do sh=shortcut(0,1),1,-1 - workload += (shortcut(sh+1,1) - shortcut(sh,1))**2 - if(workload > max_workload) then - blocke = sh - call davidson_add_task(handler, blocke, blockb) - blockb = sh-1 - workload = 0 - end if + workload = shortcut(0,1)+dble(shortcut(sh+1,1) - shortcut(sh,1))**2 + do i=sh, shortcut(0,2), shortcut(0,1) + do j=i, min(i, shortcut(0,2)) + workload += (shortcut(j+1,2) - shortcut(j, 2))**2 + end do + end do + istep = 1+ int(0.5d0*workload/ave_workload) + do blockb2=0, istep-1 + call davidson_add_task(handler, sh, blockb2, istep) + enddo enddo - if(blockb > 0) call davidson_add_task(handler, 1, blockb) - call davidson_run(handler, v_0, s_0) + call davidson_run(handler, v_0, s_0, size(v_0,1)) - do istate=1,N_st do i=1,n v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate) s_0(i,istate) = s_0(i,istate) + s2_jj(i)* u_0(i,istate) enddo enddo + deallocate(shortcut, sort_idx, sorted, version) + deallocate(ut) end -BEGIN_PROVIDER [ integer, max_workload ] - max_workload = 1000 -END_PROVIDER - diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 28b7d2e2..afc573aa 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -364,7 +364,7 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) integer(key_kind) :: idx type(map_type), intent(inout) :: map real(integral_kind) :: tmp - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_bielec_integrals_in_map mo_integrals_cache if ( (i >= mo_integrals_cache_min) .and. & (j >= mo_integrals_cache_min) .and. & (k >= mo_integrals_cache_min) .and. & @@ -393,7 +393,7 @@ double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map) integer(key_kind) :: idx type(map_type), intent(inout) :: map real(integral_kind) :: tmp - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_bielec_integrals_in_map mo_integrals_cache if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) > mo_integrals_threshold) then double precision, external :: get_mo_bielec_integral !DIR$ FORCEINLINE @@ -411,7 +411,7 @@ double precision function mo_bielec_integral(i,j,k,l) END_DOC integer, intent(in) :: i,j,k,l double precision :: get_mo_bielec_integral - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_bielec_integrals_in_map mo_integrals_cache !DIR$ FORCEINLINE mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) return diff --git a/tests/bats/cassd.bats b/tests/bats/cassd.bats new file mode 100644 index 00000000..a21b58ac --- /dev/null +++ b/tests/bats/cassd.bats @@ -0,0 +1,17 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh + +@test "CAS_SD H2O cc-pVDZ" { + test_exe cas_sd_selected || skip + INPUT=h2o.ezfio + qp_edit -c $INPUT + ezfio set_file $INPUT + ezfio set perturbation do_pt2_end False + ezfio set determinants n_det_max 1000 + qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" + qp_run cas_sd_selected $INPUT + energy="$(ezfio get cas_sd energy)" + eq $energy -76.22213389282479 1.E-5 +} + diff --git a/tests/bats/common.bats.sh b/tests/bats/common.bats.sh new file mode 100644 index 00000000..2aaff591 --- /dev/null +++ b/tests/bats/common.bats.sh @@ -0,0 +1,44 @@ +#!/usr/bin/env bats + +# floating point number comparison +# Compare two numbers ($1, $2) with a given precision ($3) +# If the numbers are not equal, the exit code is 1 else it is 0 +# So we strip the "-", is the abs value of the poor +function eq() { + declare -a diff + diff=($(awk -v d1=$1 -v d2=$2 -v n1=${1#-} -v n2=${2#-} -v p=$3 'BEGIN{ if ((n1-n2)^2 < p^2) print 0; print 1 " " (d1-d2) " " d1 " " d2 }')) + if [[ "${diff[0]}" == "0" ]] + then + return 0 + else + echo "Test : " ${BATS_TEST_DESCRIPTION} + echo "Error : " ${diff[1]} + echo "Reference : " ${diff[3]} + echo "Computed : " ${diff[2]} + exit 1 + fi +} + + +# ___ +# | ._ o _|_ +# _|_ | | | |_ +# +source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh +TEST_DIR=${QP_ROOT}/tests/work/ + +mkdir -p "${TEST_DIR}" + +cd "${TEST_DIR}" || exit 1 + +function test_exe() { + l_EXE=$(awk "/^$1 / { print \$2 }" < "${QP_ROOT}"/data/executables) + l_EXE=$(echo $l_EXE | sed "s|\$QP_ROOT|$QP_ROOT|") + if [[ -x "$l_EXE" ]] + then + return 0 + else + return 127 + fi +} + diff --git a/tests/bats/convert.bats b/tests/bats/convert.bats new file mode 100644 index 00000000..a1fbd7ce --- /dev/null +++ b/tests/bats/convert.bats @@ -0,0 +1,27 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh + +#=== Convert +@test "gamess convert HBO.out" { + cp ${QP_ROOT}/tests/input/HBO.out . + qp_convert_output_to_ezfio.py HBO.out + qp_edit -c HBO.out.ezfio + ezfio set_file HBO.out.ezfio + qp_run SCF HBO.out.ezfio + # Check energy + energy="$(ezfio get hartree_fock energy)" + eq $energy -100.0185822590964 1.e-10 +} + +@test "g09 convert H2O.log" { + cp ${QP_ROOT}/tests/input/h2o.log . + qp_convert_output_to_ezfio.py h2o.log + qp_edit -c h2o.log.ezfio + ezfio set_file h2o.log.ezfio + qp_run SCF h2o.log.ezfio + # Check energy + energy="$(ezfio get hartree_fock energy)" + eq $energy -76.0270218704265 1E-10 +} + diff --git a/tests/bats/fci.bats b/tests/bats/fci.bats new file mode 100644 index 00000000..174c8f61 --- /dev/null +++ b/tests/bats/fci.bats @@ -0,0 +1,52 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh + +function run_FCI() { + thresh=5.e-5 + test_exe full_ci || skip + qp_edit -c $1 + ezfio set_file $1 + ezfio set perturbation do_pt2_end True + ezfio set determinants n_det_max $2 + ezfio set davidson threshold_davidson 1.e-10 + + qp_run full_ci $1 + energy="$(ezfio get full_ci energy)" + eq $energy $3 $thresh + energy_pt2="$(ezfio get full_ci energy_pt2)" + eq $energy_pt2 $4 $thresh +} + +function run_FCI_ZMQ() { + thresh=5.e-5 + test_exe full_ci || skip + qp_edit -c $1 + ezfio set_file $1 + ezfio set perturbation do_pt2_end True + ezfio set determinants n_det_max $2 + ezfio set davidson threshold_davidson 1.e-10 + + qp_run fci_zmq $1 + energy="$(ezfio get full_ci energy)" + eq $energy $3 $thresh + energy_pt2="$(ezfio get full_ci energy_pt2)" + eq $energy_pt2 $4 $thresh +} + + + +#=== H2O + +@test "qp_set_mo_class H2O cc-pVDZ" { + qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]" +} +@test "FCI H2O cc-pVDZ" { + run_FCI h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02 +} + +@test "FCI-ZMQ H2O cc-pVDZ" { + run_FCI_ZMQ h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02 +} + + diff --git a/tests/bats/foboci.bats b/tests/bats/foboci.bats new file mode 100644 index 00000000..98255969 --- /dev/null +++ b/tests/bats/foboci.bats @@ -0,0 +1,27 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh + +function run_all_1h_1p() { + thresh=1.e-6 + test_exe all_1h_1p || skip + qp_edit -c $1 + ezfio set_file $1 + ezfio set determinants n_det_max $2 + ezfio set perturbation pt2_max $3 + ezfio set davidson threshold_davidson 1.e-10 + + qp_run all_1h_1p $1 | tee $1.F1h1p.out + energy="$(ezfio get all_singles energy)" + eq $energy $4 $thresh +} + + +#=== DHNO + +@test "all_1h_1p DHNO chipman-dzp" { + qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio + run_all_1h_1p dhno.ezfio 10000 0.0000000001 -130.4466283766202 +} + + diff --git a/tests/bats/hf.bats b/tests/bats/hf.bats new file mode 100644 index 00000000..e280c986 --- /dev/null +++ b/tests/bats/hf.bats @@ -0,0 +1,52 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh + +function run_init() { + cp "${QP_ROOT}/tests/input/$1" . + qp_create_ezfio_from_xyz $1 -o $3 $2 + qp_edit -c $3 +} + + +function run_HF() { + thresh=1.e-7 + test_exe SCF || skip + qp_edit -c $1 + ezfio set_file $1 + ezfio set hartree_fock thresh_scf 1.e-11 + qp_run SCF $1 + energy="$(ezfio get hartree_fock energy)" + eq $energy $2 $thresh +} + + + +#=== DHNO +@test "init DHNO chipman-dzp" { + run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio +} + +@test "SCF DHNO chipman-dzp" { + run_HF dhno.ezfio -130.4278777822 +} + +#=== HBO +@test "init HBO STO-3G" { + run_init HBO.xyz "-b STO-3G" hbo.ezfio +} + +@test "SCF HBO STO-3G" { + run_HF hbo.ezfio -98.8251985678084 +} + + +#=== H2O +@test "init H2O cc-pVDZ" { + run_init h2o.xyz "-b cc-pvdz" h2o.ezfio +} + +@test "SCF H2O cc-pVDZ" { + run_HF h2o.ezfio -0.760270218692179E+02 +} + diff --git a/tests/bats/mrcepa0.bats b/tests/bats/mrcepa0.bats new file mode 100644 index 00000000..8b56c606 --- /dev/null +++ b/tests/bats/mrcepa0.bats @@ -0,0 +1,70 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh + +#=== H2O +@test "MRCC-lambda H2O cc-pVDZ" { + INPUT=h2o.ezfio + EXE=mrcc + test_exe $EXE || skip + qp_edit -c $INPUT + ezfio set_file $INPUT + ezfio set determinants threshold_generators 1. + ezfio set determinants threshold_selectors 1. + ezfio set determinants read_wf True + ezfio set mrcepa0 lambda_type 1 + ezfio set mrcepa0 n_it_max_dressed_ci 3 + qp_run $EXE $INPUT + energy="$(ezfio get mrcepa0 energy)" + eq $energy -76.22903276183061 1.e-4 +} + +@test "MRCC H2O cc-pVDZ" { + INPUT=h2o.ezfio + EXE=mrcc + test_exe $EXE || skip + qp_edit -c $INPUT + ezfio set_file $INPUT + ezfio set determinants threshold_generators 1. + ezfio set determinants threshold_selectors 1. + ezfio set determinants read_wf True + ezfio set determinants read_wf True + ezfio set mrcepa0 lambda_type 0 + ezfio set mrcepa0 n_it_max_dressed_ci 3 + qp_run $EXE $INPUT + energy="$(ezfio get mrcepa0 energy)" + eq $energy -76.22899302846875 1.e-4 +} + +@test "MRSC2 H2O cc-pVDZ" { + INPUT=h2o.ezfio + EXE=mrsc2 + test_exe $EXE || skip + qp_edit -c $INPUT + ezfio set_file $INPUT + ezfio set determinants threshold_generators 1. + ezfio set determinants threshold_selectors 1. + ezfio set determinants read_wf True + ezfio set mrcepa0 lambda_type 0 + ezfio set mrcepa0 n_it_max_dressed_ci 3 + qp_run $EXE $INPUT + energy="$(ezfio get mrcepa0 energy)" + eq $energy -76.22647345292708 1.e-4 +} + +@test "MRCEPA0 H2O cc-pVDZ" { + INPUT=h2o.ezfio + EXE=mrcepa0 + test_exe $EXE || skip + qp_edit -c $INPUT + ezfio set_file $INPUT + ezfio set determinants threshold_generators 1. + ezfio set determinants threshold_selectors 1. + ezfio set determinants read_wf True + ezfio set mrcepa0 lambda_type 0 + ezfio set mrcepa0 n_it_max_dressed_ci 3 + qp_run $EXE $INPUT + energy="$(ezfio get mrcepa0 energy)" + eq $energy -76.23199784430074 1.e-4 +} + diff --git a/tests/bats/pseudo.bats b/tests/bats/pseudo.bats new file mode 100644 index 00000000..8cccf229 --- /dev/null +++ b/tests/bats/pseudo.bats @@ -0,0 +1,53 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh + +function run_init() { + cp "${QP_ROOT}/tests/input/$1" . + qp_create_ezfio_from_xyz $1 -o $3 $2 + qp_edit -c $3 +} + + +function run_HF() { + thresh=1.e-7 + test_exe SCF || skip + qp_edit -c $1 + ezfio set_file $1 + ezfio set hartree_fock thresh_scf 1.e-11 + qp_run SCF $1 + energy="$(ezfio get hartree_fock energy)" + eq $energy $2 $thresh +} + + +function run_FCI_ZMQ() { + thresh=5.e-5 + test_exe full_ci || skip + qp_edit -c $1 + ezfio set_file $1 + ezfio set perturbation do_pt2_end True + ezfio set determinants n_det_max $2 + ezfio set davidson threshold_davidson 1.e-10 + + qp_run fci_zmq $1 + energy="$(ezfio get full_ci energy)" + eq $energy $3 $thresh + energy_pt2="$(ezfio get full_ci energy_pt2)" + eq $energy_pt2 $4 $thresh +} + +#=== H2O Pseudo +@test "init H2O VDZ pseudo" { + run_init h2o.xyz "-p bfd -b vdz-bfd" h2o_pseudo.ezfio +} + +@test "SCF H2O VDZ pseudo" { + run_HF h2o_pseudo.ezfio -16.9483703905461 +} + +@test "FCI H2O VDZ pseudo" { + qp_set_mo_class h2o_pseudo.ezfio -core "[1]" -act "[2-12]" -del "[13-23]" + run_FCI_ZMQ h2o_pseudo.ezfio 2000 -0.170399597228904E+02 -0.170400168816800E+02 +} + diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats deleted file mode 100644 index aa283916..00000000 --- a/tests/bats/qp.bats +++ /dev/null @@ -1,205 +0,0 @@ -#!/usr/bin/env bats - -# -# |\/| o _ _ -# | | | _> (_ -# -# floating point number comparison -# Compare two numbers ($1, $2) with a given precision ($3) -# If the numbers are not equal, the exit code is 1 else it is 0 -# So we strip the "-", is the abs value of the poor -function eq() { - declare -a diff - diff=($(awk -v d1=$1 -v d2=$2 -v n1=${1#-} -v n2=${2#-} -v p=$3 'BEGIN{ if ((n1-n2)^2 < p^2) print 0; print 1 " " (d1-d2) " " d1 " " d2 }')) - if [[ "${diff[0]}" == "0" ]] - then - return 0 - else - echo "Test : " ${BATS_TEST_DESCRIPTION} - echo "Error : " ${diff[1]} - echo "Reference : " ${diff[3]} - echo "Computed : " ${diff[2]} - exit 1 - fi -} - - -# ___ -# | ._ o _|_ -# _|_ | | | |_ -# -source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh -TEST_DIR=${QP_ROOT}/tests/work/ - -mkdir -p "${TEST_DIR}" - -cd "${TEST_DIR}" || exit 1 - -function run_init() { - cp "${QP_ROOT}/tests/input/$1" . - qp_create_ezfio_from_xyz $1 -o $3 $2 - qp_edit -c $3 -} - -function test_exe() { - EXE=$(awk "/^$1 / { print \$2 }" < "${QP_ROOT}"/data/executables) - EXE=$(echo $EXE | sed "s|\$QP_ROOT|$QP_ROOT|") - if [[ -x "$EXE" ]] - then - return 0 - else - return 127 - fi -} - -function run_HF() { - thresh=1.e-7 - test_exe SCF || skip - ezfio set_file $1 - ezfio set hartree_fock thresh_scf 1.e-11 - qp_run SCF $1 - energy="$(ezfio get hartree_fock energy)" - eq $energy $2 $thresh -} - -function run_FCI() { - thresh=5.e-5 - test_exe full_ci || skip - ezfio set_file $1 - ezfio set perturbation do_pt2_end True - ezfio set determinants n_det_max $2 - ezfio set davidson threshold_davidson 1.e-10 - - qp_run full_ci $1 - energy="$(ezfio get full_ci energy)" - eq $energy $3 $thresh - energy_pt2="$(ezfio get full_ci energy_pt2)" - eq $energy_pt2 $4 $thresh -} - -function run_all_1h_1p() { - thresh=1.e-6 - test_exe all_1h_1p || skip - ezfio set_file $1 - ezfio set determinants n_det_max $2 - ezfio set perturbation pt2_max $3 - ezfio set davidson threshold_davidson 1.e-10 - - qp_run all_1h_1p $1 | tee $1.F1h1p.out - energy="$(ezfio get all_singles energy)" - eq $energy $4 $thresh -} - -# ___ -# | _ _ _|_ -# | (/_ _> |_ -# - - -#=== DHNO -@test "init DHNO chipman-dzp" { - run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio -} - -@test "SCF DHNO chipman-dzp" { - run_HF dhno.ezfio -130.4278777822 -} - -@test "all_1h_1p DHNO chipman-dzp" { - qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio - run_all_1h_1p dhno.ezfio 10000 0.0000000001 -130.4466283766202 -} - -#=== HBO -@test "init HBO STO-3G" { - run_init HBO.xyz "-b STO-3G" hbo.ezfio -} - -@test "SCF HBO STO-3G" { - run_HF hbo.ezfio -98.8251985678084 -} - - -#=== H2O -@test "init H2O cc-pVDZ" { - run_init h2o.xyz "-b cc-pvdz" h2o.ezfio -} - -@test "SCF H2O cc-pVDZ" { - run_HF h2o.ezfio -0.760270218692179E+02 -} - -@test "FCI H2O cc-pVDZ" { - qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]" - run_FCI h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02 -} - -@test "CAS_SD H2O cc-pVDZ" { - test_exe cas_sd_selected || skip - INPUT=h2o.ezfio - ezfio set_file $INPUT - ezfio set perturbation do_pt2_end False - ezfio set determinants n_det_max 1000 - qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" - qp_run cas_sd_selected $INPUT - energy="$(ezfio get cas_sd energy)" - eq $energy -0.762219854008117E+02 1.E-5 -} - -@test "MRCC H2O cc-pVDZ" { - test_exe mrcc_cassd || skip - INPUT=h2o.ezfio - ezfio set_file $INPUT - ezfio set determinants threshold_generators 1. - ezfio set determinants threshold_selectors 1. - ezfio set determinants read_wf True - qp_run mrcc_cassd $INPUT - energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2288648023833 1.e-4 - -} - - -#=== H2O Pseudo -@test "init H2O VDZ pseudo" { - run_init h2o.xyz "-p bfd -b vdz-bfd" h2o_pseudo.ezfio -} - -@test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -16.9483703905461 -} - -@test "FCI H2O VDZ pseudo" { - qp_set_mo_class h2o_pseudo.ezfio -core "[1]" -act "[2-12]" -del "[13-23]" - run_FCI h2o_pseudo.ezfio 2000 -0.170399597228904E+02 -0.170400168816800E+02 -} - -#=== Convert -@test "gamess convert HBO.out" { - cp ${QP_ROOT}/tests/input/HBO.out . - qp_convert_output_to_ezfio.py HBO.out - ezfio set_file HBO.out.ezfio - qp_run SCF HBO.out.ezfio - # Check energy - energy="$(ezfio get hartree_fock energy)" - eq $energy -100.0185822590964 1.e-10 -} - -@test "g09 convert H2O.log" { - cp ${QP_ROOT}/tests/input/h2o.log . - qp_convert_output_to_ezfio.py h2o.log - ezfio set_file h2o.log.ezfio - qp_run SCF h2o.log.ezfio - # Check energy - energy="$(ezfio get hartree_fock energy)" - eq $energy -76.0270218704265 1E-10 -} - - -# TODO N_int = 1,2,3,4,5 -# TODO mod(64) MOs -# TODO All G2 SCF energies -# TODO Long and short tests -# TODO MP2 -# TODO CISD_selected - diff --git a/tests/bats_to_sh.py b/tests/bats_to_sh.py index 2c6b4a05..8feb9272 100755 --- a/tests/bats_to_sh.py +++ b/tests/bats_to_sh.py @@ -1,6 +1,8 @@ #!/usr/bin/env python -with open('bats/qp.bats','r') as f: +import sys + +with open(sys.argv[1],'r') as f: raw_data = f.read() output = [] diff --git a/tests/run_tests.sh b/tests/run_tests.sh index 2436c60c..4664ce82 100755 --- a/tests/run_tests.sh +++ b/tests/run_tests.sh @@ -1,18 +1,39 @@ #!/bin/bash +LIST=" + +convert.bats +hf.bats +foboci.bats +pseudo.bats +fci.bats +cassd.bats +mrcepa0.bats + +" + + export QP_PREFIX="timeout -s 9 300" export QP_TASK_DEBUG=1 -BATS_FILE=bats/qp.bats - rm -rf work output -if [[ "$1" == "-v" ]] -then - echo "Verbose mode" - ./bats_to_sh.py $BATS_FILE | bash -else - bats $BATS_FILE -fi + +for BATS_FILE in $LIST +do + echo + echo "-~-~-~-~-~-~" + echo + echo "Running tests for ${BATS_FILE%.bats}" + echo + BATS_FILE=bats/$BATS_FILE + if [[ "$1" == "-v" ]] + then + echo "Verbose mode" + ./bats_to_sh.py $BATS_FILE | bash + else + bats $BATS_FILE + fi +done