diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 5f0c4771..cbc81032 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -335,7 +335,7 @@ def write_ezfio(res, filename): def get_full_path(file_path): file_path = os.path.expanduser(file_path) file_path = os.path.expandvars(file_path) - file_path = os.path.abspath(file_path) +# file_path = os.path.abspath(file_path) return file_path diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 486ff534..4108ce71 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -3,6 +3,8 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] ! Nucleus-electron interaction, in the |AO| basis set. ! ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` + ! + ! These integrals also contain the pseudopotential integrals. END_DOC implicit none double precision :: alpha, beta, gama, delta @@ -75,11 +77,11 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] !$OMP END DO !$OMP END PARALLEL - endif + IF (DO_PSEUDO) THEN + ao_integrals_n_e += ao_pseudo_integrals + ENDIF - IF (DO_PSEUDO) THEN - ao_integrals_n_e += ao_pseudo_integrals - ENDIF + endif if (write_ao_integrals_n_e) then diff --git a/src/ao_one_e_ints/screening.irp.f b/src/ao_one_e_ints/screening.irp.f index aa23a9b8..1bbe3c73 100644 --- a/src/ao_one_e_ints/screening.irp.f +++ b/src/ao_one_e_ints/screening.irp.f @@ -3,14 +3,11 @@ logical function ao_one_e_integral_zero(i,k) integer, intent(in) :: i,k ao_one_e_integral_zero = .False. - if (.not.(read_ao_one_e_integrals.or.is_periodic)) then + if (.not.((io_ao_integrals_overlap/='None').or.is_periodic)) then if (ao_overlap_abs(i,k) < ao_integrals_threshold) then ao_one_e_integral_zero = .True. return endif endif - if (ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then - ao_one_e_integral_zero = .True. - endif end diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f index 5095deec..d3230370 100644 --- a/src/ao_two_e_ints/screening.irp.f +++ b/src/ao_two_e_ints/screening.irp.f @@ -8,8 +8,8 @@ logical function ao_two_e_integral_zero(i,j,k,l) ao_two_e_integral_zero = .True. return endif - endif - if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then - ao_two_e_integral_zero = .True. + if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then + ao_two_e_integral_zero = .True. + endif endif end diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index acc9ab7a..b6e959d7 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -18,89 +18,89 @@ double precision function ao_two_e_integral(i,j,k,l) if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l) - return - endif + else - dim1 = n_pt_max_integrals + dim1 = n_pt_max_integrals - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral = 0.d0 + num_i = ao_nucl(i) + num_j = ao_nucl(j) + num_k = ao_nucl(k) + num_l = ao_nucl(l) + ao_two_e_integral = 0.d0 - if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) - J_center(p) = nucl_coord(num_j,p) - K_center(p) = nucl_coord(num_k,p) - L_center(p) = nucl_coord(num_l,p) - enddo + if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + I_center(p) = nucl_coord(num_i,p) + J_center(p) = nucl_coord(num_j,p) + K_center(p) = nucl_coord(num_k,p) + L_center(p) = nucl_coord(num_l,p) + enddo - double precision :: coef1, coef2, coef3, coef4 - double precision :: p_inv,q_inv - double precision :: general_primitive_integral + double precision :: coef1, coef2, coef3, coef4 + double precision :: p_inv,q_inv + double precision :: general_primitive_integral - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & - I_power,J_power,I_center,J_center,dim1) - p_inv = 1.d0/pp - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - integral = general_primitive_integral(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_two_e_integral = ao_two_e_integral + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p + do p = 1, ao_prim_num(i) + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) + call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & + I_power,J_power,I_center,J_center,dim1) + p_inv = 1.d0/pp + do r = 1, ao_prim_num(k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) + call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& + ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & + K_power,L_power,K_center,L_center,dim1) + q_inv = 1.d0/qq + integral = general_primitive_integral(dim1, & + P_new,P_center,fact_p,pp,p_inv,iorder_p, & + Q_new,Q_center,fact_q,qq,q_inv,iorder_q) + ao_two_e_integral = ao_two_e_integral + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p - else + else - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - enddo - double precision :: ERI + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + enddo + double precision :: ERI - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - integral = ERI( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - I_power(1),J_power(1),K_power(1),L_power(1), & - I_power(2),J_power(2),K_power(2),L_power(2), & - I_power(3),J_power(3),K_power(3),L_power(3)) - ao_two_e_integral = ao_two_e_integral + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p + do p = 1, ao_prim_num(i) + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) + do r = 1, ao_prim_num(k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) + integral = ERI( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& + I_power(1),J_power(1),K_power(1),L_power(1), & + I_power(2),J_power(2),K_power(2),L_power(2), & + I_power(3),J_power(3),K_power(3),L_power(3)) + ao_two_e_integral = ao_two_e_integral + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif endif - end double precision function ao_two_e_integral_schwartz_accel(i,j,k,l) @@ -343,8 +343,6 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] integer :: kk, m, j1, i1, lmax character*(64) :: fmt - integral = ao_two_e_integral(1,1,1,1) - double precision :: map_mb PROVIDE read_ao_two_e_integrals io_ao_two_e_integrals if (read_ao_two_e_integrals) then @@ -352,66 +350,72 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) print*, 'AO integrals provided' ao_two_e_integrals_in_map = .True. - return - endif + else - print*, 'Providing the AO integrals' - call wall_time(wall_0) - call wall_time(wall_1) - call cpu_time(cpu_1) + print*, 'Providing the AO integrals' + call wall_time(wall_0) + call wall_time(wall_1) + call cpu_time(cpu_1) - integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull - call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals') - - character(len=:), allocatable :: task - allocate(character(len=ao_num*12) :: task) - write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' - do l=1,ao_num - write(task,fmt) (i,l, i=1,l) - integer, external :: add_task_to_taskserver - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then - stop 'Unable to add task to server' + if (.True.) then + ! Avoid openMP + integral = ao_two_e_integral(1,1,1,1) endif - enddo - deallocate(task) - integer, external :: zmq_set_running - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif + integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull + call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals') - PROVIDE nproc - !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1) - i = omp_get_thread_num() - if (i==0) then - call ao_two_e_integrals_in_map_collector(zmq_socket_pull) - else - call ao_two_e_integrals_in_map_slave_inproc(i) + character(len=:), allocatable :: task + allocate(character(len=ao_num*12) :: task) + write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' + do l=1,ao_num + write(task,fmt) (i,l, i=1,l) + integer, external :: add_task_to_taskserver + if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then + stop 'Unable to add task to server' endif - !$OMP END PARALLEL + enddo + deallocate(task) - call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals') + integer, external :: zmq_set_running + if (zmq_set_running(zmq_to_qp_run_socket) == -1) then + print *, irp_here, ': Failed in zmq_set_running' + endif + + PROVIDE nproc + !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call ao_two_e_integrals_in_map_collector(zmq_socket_pull) + else + call ao_two_e_integrals_in_map_slave_inproc(i) + endif + !$OMP END PARALLEL + + call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals') - print*, 'Sorting the map' - call map_sort(ao_integrals_map) - call cpu_time(cpu_2) - call wall_time(wall_2) - integer(map_size_kind) :: get_ao_map_size, ao_map_size - ao_map_size = get_ao_map_size() + print*, 'Sorting the map' + call map_sort(ao_integrals_map) + call cpu_time(cpu_2) + call wall_time(wall_2) + integer(map_size_kind) :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() - print*, 'AO integrals provided:' - print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' - print*, ' Number of AO integrals :', ao_map_size - print*, ' cpu time :',cpu_2 - cpu_1, 's' - print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' + print*, 'AO integrals provided:' + print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' + print*, ' Number of AO integrals :', ao_map_size + print*, ' cpu time :',cpu_2 - cpu_1, 's' + print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' - ao_two_e_integrals_in_map = .True. + ao_two_e_integrals_in_map = .True. + + if (write_ao_two_e_integrals.and.mpi_master) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + endif - if (write_ao_two_e_integrals.and.mpi_master) then - call ezfio_set_work_empty(.False.) - call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) - call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') endif END_PROVIDER diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg index ca2100a1..17e8fd90 100644 --- a/src/becke_numerical_grid/EZFIO.cfg +++ b/src/becke_numerical_grid/EZFIO.cfg @@ -14,3 +14,22 @@ type: double precision doc: threshold on the weight of a given grid point interface: ezfio,provider,ocaml default: 1.e-20 + +[my_grid_becke] +type: logical +doc: if True, the number of angular and radial grid points are read from EZFIO +interface: ezfio,provider,ocaml +default: False + +[my_n_pt_r_grid] +type: integer +doc: Number of radial grid points given from input +interface: ezfio,provider,ocaml +default: 300 + +[my_n_pt_a_grid] +type: integer +doc: Number of angular grid points given from input. Warning, this number cannot be any integer. See file list_angular_grid +interface: ezfio,provider,ocaml +default: 1202 + diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f index e72f6460..bedf0a84 100644 --- a/src/becke_numerical_grid/grid_becke.irp.f +++ b/src/becke_numerical_grid/grid_becke.irp.f @@ -8,7 +8,8 @@ ! ! These numbers are automatically set by setting the grid_type_sgn parameter END_DOC -select case (grid_type_sgn) +if(.not.my_grid_becke)then + select case (grid_type_sgn) case(0) n_points_radial_grid = 23 n_points_integration_angular = 170 @@ -25,6 +26,10 @@ select case (grid_type_sgn) write(*,*) '!!! Quadrature grid not available !!!' stop end select +else + n_points_radial_grid = my_n_pt_r_grid + n_points_integration_angular = my_n_pt_a_grid +endif END_PROVIDER BEGIN_PROVIDER [integer, n_points_grid_per_atom] diff --git a/src/becke_numerical_grid/list_angular_grid b/src/becke_numerical_grid/list_angular_grid new file mode 100644 index 00000000..252c7432 --- /dev/null +++ b/src/becke_numerical_grid/list_angular_grid @@ -0,0 +1,32 @@ +0006 +0014 +0026 +0038 +0050 +0074 +0086 +0110 +0146 +0170 +0194 +0230 +0266 +0302 +0350 +0434 +0590 +0770 +0974 +1202 +1454 +1730 +2030 +2354 +2702 +3074 +3470 +3890 +4334 +4802 +5294 +5810 diff --git a/src/cas_based_on_top/two_body_dens_rout.irp.f b/src/cas_based_on_top/two_body_dens_rout.irp.f new file mode 100644 index 00000000..4a57a868 --- /dev/null +++ b/src/cas_based_on_top/two_body_dens_rout.irp.f @@ -0,0 +1,152 @@ + +subroutine give_n2_ii_val_ab(r1,r2,two_bod_dens) + implicit none + BEGIN_DOC +! contribution from purely inactive orbitals to n2_{\Psi^B}(r_1,r_2) for a CAS wave function + END_DOC + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: two_bod_dens + integer :: i,j,m,n,i_m,i_n + integer :: i_i,i_j + double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + ! You extract the inactive orbitals + allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) + do i_m = 1, n_inact_orb + mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m)) + enddo + do i_m = 1, n_inact_orb + mos_array_inact_r2(i_m) = mos_array_r2(list_inact(i_m)) + enddo + + two_bod_dens = 0.d0 + ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space + do m = 1, n_inact_orb ! electron 1 + ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space + do n = 1, n_inact_orb ! electron 2 + ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2) + two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) + enddo + enddo +end + + +subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate) + BEGIN_DOC +! contribution from inactive and active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function + END_DOC + implicit none + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: two_bod_dens + integer :: i,orb_i,a,orb_a,n,m,b + double precision :: rho + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) + double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) + + two_bod_dens = 0.d0 + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + + ! You extract the inactive orbitals + allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) + do i = 1, n_inact_orb + mos_array_inact_r1(i) = mos_array_r1(list_inact(i)) + enddo + do i= 1, n_inact_orb + mos_array_inact_r2(i) = mos_array_r2(list_inact(i)) + enddo + + ! You extract the active orbitals + allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) ) + do i= 1, n_act_orb + mos_array_act_r1(i) = mos_array_r1(list_act(i)) + enddo + do i= 1, n_act_orb + mos_array_act_r2(i) = mos_array_r2(list_act(i)) + enddo + + ! Contracted density : intermediate quantity + two_bod_dens = 0.d0 + do a = 1, n_act_orb + do i = 1, n_inact_orb + do b = 1, n_act_orb + rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) + two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho + enddo + enddo + enddo +end + + +subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate) + BEGIN_DOC +! contribution from purely active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function + END_DOC + implicit none + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: two_bod_dens + integer :: i,orb_i,a,orb_a,n,m,b,c,d + double precision :: rho + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) + + two_bod_dens = 0.d0 + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + + ! You extract the active orbitals + allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) ) + do i= 1, n_act_orb + mos_array_act_r1(i) = mos_array_r1(list_act(i)) + enddo + do i= 1, n_act_orb + mos_array_act_r2(i) = mos_array_r2(list_act(i)) + enddo + + ! Contracted density : intermediate quantity + two_bod_dens = 0.d0 + do a = 1, n_act_orb ! 1 + do b = 1, n_act_orb ! 2 + do c = 1, n_act_orb ! 1 + do d = 1, n_act_orb ! 2 + rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate) + two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) + enddo + enddo + enddo + enddo + +end + +subroutine give_n2_cas(r1,r2,istate,n2_psi) + implicit none + BEGIN_DOC +! returns mu(r), f_psi, n2_psi for a general cas wave function + END_DOC + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out) :: n2_psi + double precision :: two_bod_dens_ii + double precision :: two_bod_dens_ia + double precision :: two_bod_dens_aa + ! inactive-inactive part of n2_psi(r1,r2) + call give_n2_ii_val_ab(r1,r2,two_bod_dens_ii) + ! inactive-active part of n2_psi(r1,r2) + call give_n2_ia_val_ab(r1,r2,two_bod_dens_ia,istate) + ! active-active part of n2_psi(r1,r2) + call give_n2_aa_val_ab(r1,r2,two_bod_dens_aa,istate) + + n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa + +end diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index d83aa271..950cfd55 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -5,6 +5,7 @@ program casscf END_DOC call reorder_orbitals_for_casscf no_vvvv_integrals = .True. + touch no_vvvv_integrals pt2_max = 0.02 SOFT_TOUCH no_vvvv_integrals pt2_max call run_stochastic_cipsi diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index c0d94b35..aed81063 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -438,6 +438,11 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) ipos=1 do imin=1,N_det,tasksize imax = min(N_det,imin-1+tasksize) + if (imin==1) then + istep = 2 + else + istep = 1 + endif do ishift=0,istep-1 write(task(ipos:ipos+50),'(4(I11,1X),1X,1A)') imin, imax, ishift, istep, '|' ipos = ipos+50 diff --git a/src/ezfio_files/ezfio.irp.f b/src/ezfio_files/ezfio.irp.f index e7a7a159..4f53b173 100644 --- a/src/ezfio_files/ezfio.irp.f +++ b/src/ezfio_files/ezfio.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ character*(128), ezfio_filename ] +BEGIN_PROVIDER [ character*(1024), ezfio_filename ] implicit none BEGIN_DOC ! Name of EZFIO file. It is obtained from the QPACKAGE_INPUT environment @@ -34,7 +34,7 @@ BEGIN_PROVIDER [ character*(128), ezfio_filename ] ! Adjust out-of-memory killer flag such that the current process will be ! killed first by the OOM killer, allowing compute nodes to survive integer :: getpid - character*(64) :: command, pidc + character*(1024) :: command, pidc write(pidc,*) getpid() write(command,*) 'echo 15 > /proc//'//trim(adjustl(pidc))//'/oom_adj' call system(command) @@ -43,7 +43,7 @@ BEGIN_PROVIDER [ character*(128), ezfio_filename ] END_PROVIDER -BEGIN_PROVIDER [ character*(128), ezfio_work_dir ] +BEGIN_PROVIDER [ character*(1024), ezfio_work_dir ] implicit none BEGIN_DOC ! EZFIO/work/ diff --git a/src/ezfio_files/get_unit_and_open.irp.f b/src/ezfio_files/get_unit_and_open.irp.f index 71bbc76c..6440579f 100644 --- a/src/ezfio_files/get_unit_and_open.irp.f +++ b/src/ezfio_files/get_unit_and_open.irp.f @@ -17,7 +17,7 @@ integer function getUnitAndOpen(f,mode) END_DOC character*(*) :: f - character*(128) :: new_f + character*(256) :: new_f integer :: iunit logical :: is_open, exists character :: mode diff --git a/src/ezfio_files/qp_stop.irp.f b/src/ezfio_files/qp_stop.irp.f index b2618246..c1039dc3 100644 --- a/src/ezfio_files/qp_stop.irp.f +++ b/src/ezfio_files/qp_stop.irp.f @@ -1,5 +1,5 @@ - BEGIN_PROVIDER [ character*(128), qp_stop_filename ] -&BEGIN_PROVIDER [ character*(128), qp_kill_filename ] + BEGIN_PROVIDER [ character*(256), qp_stop_filename ] +&BEGIN_PROVIDER [ character*(256), qp_kill_filename ] &BEGIN_PROVIDER [ integer, qp_stop_variable ] implicit none BEGIN_DOC diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index bec74552..ea47c51c 100644 --- a/src/mo_two_e_ints/EZFIO.cfg +++ b/src/mo_two_e_ints/EZFIO.cfg @@ -11,3 +11,9 @@ interface: ezfio,provider,ocaml default: 1.e-15 ezfio_name: threshold_mo +[no_vvvv_integrals] +type: logical +doc: If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices +interface: ezfio,provider,ocaml +default: false + diff --git a/src/mo_two_e_ints/four_idx_novvvv.irp.f b/src/mo_two_e_ints/four_idx_novvvv.irp.f index 054d0a35..a4cdfd52 100644 --- a/src/mo_two_e_ints/four_idx_novvvv.irp.f +++ b/src/mo_two_e_ints/four_idx_novvvv.irp.f @@ -1,11 +1,11 @@ -BEGIN_PROVIDER [ logical, no_vvvv_integrals ] - implicit none - BEGIN_DOC +!BEGIN_PROVIDER [ logical, no_vvvv_integrals ] +! implicit none +! BEGIN_DOC ! If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices - END_DOC - - no_vvvv_integrals = .False. -END_PROVIDER +! END_DOC +! +! no_vvvv_integrals = .False. +!END_PROVIDER BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ] implicit none @@ -56,6 +56,8 @@ subroutine four_idx_novvvv BEGIN_DOC ! Retransform MO integrals for next CAS-SCF step END_DOC + print*,'Using partial transformation' + print*,'It will not transform all integrals with at least 3 indices within the virtuals' integer :: i,j,k,l,n_integrals double precision, allocatable :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:) double precision, external :: get_ao_two_e_integral diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index a9983e51..b926d688 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -189,7 +189,6 @@ subroutine add_integrals_to_map(mask_ijkl) two_e_tmp_2 = 0.d0 do j1 = 1,ao_num call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) - ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) enddo do j1 = 1,ao_num kmax = 0 @@ -747,7 +746,6 @@ subroutine add_integrals_to_map_no_exit_34(mask_ijkl) two_e_tmp_2 = 0.d0 do j1 = 1,ao_num call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) - ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) enddo do j1 = 1,ao_num kmax = 0 diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index 88ee24e8..39e6b399 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -211,7 +211,7 @@ END_PROVIDER END_DOC integer :: iunit, i integer, external :: getUnitAndOpen - character*(128) :: filename + character*(1024) :: filename if (mpi_master) then call getenv('QP_ROOT',filename) filename = trim(filename)//'/data/list_element.txt' diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index 2473d3c8..e3265572 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -2,6 +2,8 @@ BEGIN_PROVIDER [double precision, act_2_rdm_ab_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none BEGIN_DOC +! 12 12 +! 1 2 1 2 == ! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons ! ! diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index 3ff1bb42..c907e425 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -75,7 +75,6 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, P_new(0,1) = 0.d0 P_new(0,2) = 0.d0 P_new(0,3) = 0.d0 - !DIR$ FORCEINLINE call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) if (fact_k < thresh) then diff --git a/src/utils/shank.irp.f b/src/utils/shank.irp.f new file mode 100644 index 00000000..d459399f --- /dev/null +++ b/src/utils/shank.irp.f @@ -0,0 +1,54 @@ +double precision function shank_general(array,n,nmax) + implicit none + integer, intent(in) :: n,nmax + double precision, intent(in) :: array(0:nmax) ! array of the partial sums + integer :: ntmp,i + double precision :: sum(0:nmax),shank1(0:nmax) + if(n.lt.3)then + print*,'You asked to Shank a sum but the order is smaller than 3 ...' + print*,'n = ',n + print*,'stopping ....' + stop + endif + ntmp = n + sum = array + i = 0 + do while(ntmp.ge.2) + i += 1 +! print*,'i = ',i + call shank(sum,ntmp,nmax,shank1) + ntmp = ntmp - 2 + sum = shank1 + shank_general = shank1(ntmp) + enddo +end + + +subroutine shank(array,n,nmax,shank1) + implicit none + integer, intent(in) :: n,nmax + double precision, intent(in) :: array(0:nmax) + double precision, intent(out) :: shank1(0:nmax) + integer :: i,j + double precision :: shank_function + do i = 1, n-1 + shank1(i-1) = shank_function(array,i,nmax) + enddo +end + +double precision function shank_function(array,i,n) + implicit none + integer, intent(in) :: i,n + double precision, intent(in) :: array(0:n) + double precision :: b_n, b_n1 + b_n = array(i) - array(i-1) + b_n1 = array(i+1) - array(i) + if(dabs(b_n1-b_n).lt.1.d-12)then + shank_function = array(i+1) + else + shank_function = array(i+1) - b_n1*b_n1/(b_n1-b_n) + endif + +end + + diff --git a/src/zmq/utils.irp.f b/src/zmq/utils.irp.f index 287a6024..503dedd2 100644 --- a/src/zmq/utils.irp.f +++ b/src/zmq/utils.irp.f @@ -585,7 +585,7 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in) stop 'Wrong end of job' endif - do i=1200,1,-1 + do i=360,1,-1 rc = f77_zmq_send(zmq_to_qp_run_socket, 'end_job '//trim(zmq_state),8+len(trim(zmq_state)),0) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 512, 0) if (trim(message(1:13)) == 'error waiting') then