From b5637661fa80c3277df196d52a8ceb4cc0a38106 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 15 Nov 2023 13:10:14 +0100 Subject: [PATCH 1/8] Fix shell_index when converting file --- bin/qp_convert_output_to_ezfio | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 091423e4..32721e1e 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -224,14 +224,18 @@ def write_ezfio(res, filename): exponent += [p.expo for p in b.prim] ang_mom.append(str.count(s, "z")) shell_prim_num.append(len(b.prim)) - shell_index += [nshell_tot+1] * len(b.prim) + shell_index += [nshell_tot] * len(b.prim) + + shell_num = len(ang_mom) + assert(shell_index[0] = 1) + assert(shell_index[-1] = shell_num) # ~#~#~#~#~ # # W r i t e # # ~#~#~#~#~ # ezfio.set_basis_basis("Read from ResultsFile") - ezfio.set_basis_shell_num(len(ang_mom)) + ezfio.set_basis_shell_num(shell_num) ezfio.set_basis_basis_nucleus_index(nucl_index) ezfio.set_basis_prim_num(len(coefficient)) From ad1fd55fe9b969f1476e43a2010cc59f57ac9a8c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Feb 2024 13:08:02 +0100 Subject: [PATCH 2/8] Add mo_symmetry --- src/mo_basis/EZFIO.cfg | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg index 4c4f1eca..8349c006 100644 --- a/src/mo_basis/EZFIO.cfg +++ b/src/mo_basis/EZFIO.cfg @@ -32,6 +32,12 @@ doc: |MO| occupation numbers interface: ezfio size: (mo_basis.mo_num) +[mo_symmetry] +type: integer +doc: MOs with the same integer belong to the same irrep. +interface: ezfio +size: (mo_basis.mo_num) + [mo_class] type: MO_class doc: [ Core | Inactive | Active | Virtual | Deleted ], as defined by :ref:`qp_set_mo_class` From 9c49be2f593a53db6eb9dca3ace347559d307a68 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Feb 2024 13:16:55 +0100 Subject: [PATCH 3/8] Added mo_symmetry in qp_convert --- bin/qp_convert_output_to_ezfio | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 0523b6a7..0b8484f6 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -309,10 +309,19 @@ def write_ezfio(res, filename): MoMatrix = [] sym0 = [i.sym for i in res.mo_sets[MO_type]] - sym = [i.sym for i in res.mo_sets[MO_type]] + sym = [i.sym for i in res.mo_sets[MO_type]] for i in range(len(sym)): sym[MOmap[i]] = sym0[i] + irrep = {} + for i in sym: + irrep[i] = 0 + + for i, j in enumerate(irrep.keys()): + irrep[j] = i+1 + + sym = [ irrep[k] for k in sym ] + MoMatrix = [] for i in range(len(MOs)): m = MOs[i] @@ -329,6 +338,7 @@ def write_ezfio(res, filename): ezfio.set_mo_basis_mo_num(mo_num) ezfio.set_mo_basis_mo_coef(MoMatrix) ezfio.set_mo_basis_mo_occ(OccNum) + ezfio.set_mo_basis_mo_symmetry(sym) print("OK") From b9932c0e77b1b07e242e4d55c4679603b7bb964a Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 26 Feb 2024 15:33:36 +0100 Subject: [PATCH 4/8] added swaping between Left/Right MOs when large angles --- plugins/local/tc_scf/routines_rotates.irp.f | 70 +++++++++++++++++++-- src/mo_one_e_ints/spread_dipole_mo.irp.f | 18 ++++++ src/utils/linear_algebra.irp.f | 15 +++-- 3 files changed, 94 insertions(+), 9 deletions(-) diff --git a/plugins/local/tc_scf/routines_rotates.irp.f b/plugins/local/tc_scf/routines_rotates.irp.f index cc825429..c42e846e 100644 --- a/plugins/local/tc_scf/routines_rotates.irp.f +++ b/plugins/local/tc_scf/routines_rotates.irp.f @@ -103,7 +103,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:) double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:) - E_thr = 1d-8 + E_thr = 1d-04 E_old = TC_HF_energy allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num)) mo_r_coef_old = mo_r_coef @@ -164,10 +164,42 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen)) allocate(T(n_degen,n_degen), Snew(n_degen,n_degen)) - do j = 1, n_degen - mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j)) - mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j)) - enddo + print*,'Right orbitals before' + do j = 1, n_degen + write(*,'(100(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j)) + enddo + print*,'Left orbitals before' + do j = 1, n_degen + write(*,'(100(F16.10,X))')mo_l_coef(1:ao_num,list_degen(i,j)) + enddo + if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then + integer :: i_list, j_list + i_list = list_degen(i,1) + j_list = list_degen(i,2) + print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(list_degen(i,2)) + print*,'i_list = ',i_list + print*,'i_list = ',j_list + print*,'Swapping left/right orbitals' + call print_strong_overlap(i_list, j_list) + mo_r_coef_tmp(1:ao_num,1) = mo_r_coef_new(1:ao_num,i_list) + mo_r_coef_tmp(1:ao_num,2) = mo_l_coef(1:ao_num,i_list) + mo_l_coef_tmp(1:ao_num,1) = mo_l_coef(1:ao_num,j_list) + mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list) + else + do j = 1, n_degen + print*,'i_list = ',list_degen(i,j) + mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j)) + mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j)) + enddo + endif + print*,'Right orbitals ' + do j = 1, n_degen + write(*,'(100(F16.10,X))')mo_r_coef_tmp(1:ao_num,j) + enddo + print*,'Left orbitals ' + do j = 1, n_degen + write(*,'(100(F16.10,X))')mo_l_coef_tmp(1:ao_num,j) + enddo ! Orthogonalization of right functions print *, ' Orthogonalization of RIGHT functions' print *, ' ------------------------------------' @@ -445,3 +477,31 @@ subroutine sort_by_tc_fock end + +subroutine print_strong_overlap(i_list, j_list) + implicit none + integer, intent(in) :: i_list,j_list + double precision :: o_i, o_j,o_ij + double precision :: s_mat_r(2,2),s_mat_l(2,2) + o_i = dsqrt(overlap_mo_r(i_list, i_list)) + o_j = dsqrt(overlap_mo_r(j_list, j_list)) + o_ij = overlap_mo_r(j_list, i_list) + s_mat_r(1,1) = o_i*o_i + s_mat_r(2,1) = o_ij/(o_i * o_j) + s_mat_r(2,2) = o_j*o_j + s_mat_r(1,2) = s_mat_r(2,1) + print*,'Right overlap matrix ' + write(*,'(2(F10.5,X))')s_mat_r(1:2,1) + write(*,'(2(F10.5,X))')s_mat_r(1:2,2) + o_i = dsqrt(overlap_mo_l(i_list, i_list)) + o_j = dsqrt(overlap_mo_l(j_list, j_list)) + o_ij = overlap_mo_l(j_list, i_list) + s_mat_l(1,1) = o_i*o_i + s_mat_l(2,1) = o_ij/(o_i * o_j) + s_mat_l(2,2) = o_j*o_j + s_mat_l(1,2) = s_mat_l(2,1) + print*,'Left overlap matrix ' + write(*,'(2(F10.5,X))')s_mat_l(1:2,1) + write(*,'(2(F10.5,X))')s_mat_l(1:2,2) + +end diff --git a/src/mo_one_e_ints/spread_dipole_mo.irp.f b/src/mo_one_e_ints/spread_dipole_mo.irp.f index e4484433..b0a7198b 100644 --- a/src/mo_one_e_ints/spread_dipole_mo.irp.f +++ b/src/mo_one_e_ints/spread_dipole_mo.irp.f @@ -58,3 +58,21 @@ END_PROVIDER ) END_PROVIDER + BEGIN_PROVIDER [double precision, mo_spread_centered_x, (mo_num, mo_num) ] +&BEGIN_PROVIDER [double precision, mo_spread_centered_y, (mo_num, mo_num) ] +&BEGIN_PROVIDER [double precision, mo_spread_centered_z, (mo_num, mo_num) ] + BEGIN_DOC + ! array of the integrals of MO_i * (x^2 - ^2) MO_j = MO_i x^2 MO_j - (MO_i x MO_j)^2 + ! array of the integrals of MO_i * (y^2 - ^2) MO_j = MO_i y^2 MO_j - (MO_i y MO_j)^2 + ! array of the integrals of MO_i * (z^2 - ^2) MO_j = MO_i z^2 MO_j - (MO_i z MO_j)^2 + END_DOC + implicit none + integer :: i,j + do i = 1, mo_num + do j = 1, mo_num + mo_spread_centered_x(j,i) = mo_spread_x(j,i) - mo_dipole_x(j,i)**2 + mo_spread_centered_y(j,i) = mo_spread_y(j,i) - mo_dipole_y(j,i)**2 + mo_spread_centered_z(j,i) = mo_spread_z(j,i) - mo_dipole_z(j,i)**2 + enddo + enddo +END_PROVIDER diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 175beff3..26e390b7 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1920,8 +1920,12 @@ subroutine exp_matrix(X,n,exp_X) call get_A_squared(X,n,r2_mat) call lapack_diagd(eigvalues,eigvectors,r2_mat,n,n) eigvalues=-eigvalues + do i = 1,n + ! t = dsqrt(t^2) where t^2 are eigenvalues of X^2 + eigvalues(i) = dsqrt(eigvalues(i)) + enddo - if(.False.)then + if(.false.)then !!! For debugging and following the book intermediate ! rebuilding the matrix : X^2 = -W t^2 W^T as in 3.1.30 ! matrix_tmp1 = W t^2 @@ -1932,14 +1936,16 @@ subroutine exp_matrix(X,n,exp_X) enddo eigvalues_mat=0.d0 do i = 1,n - ! t = dsqrt(t^2) where t^2 are eigenvalues of X^2 - eigvalues(i) = dsqrt(eigvalues(i)) eigvalues_mat(i,i) = eigvalues(i)*eigvalues(i) enddo call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) + print*,'r2_mat = ' + do i = 1, n + write(*,'(100(F16.10,X))')r2_mat(:,i) + enddo print*,'r2_mat new = ' do i = 1, n write(*,'(100(F16.10,X))')matrix_tmp2(:,i) @@ -1964,7 +1970,8 @@ subroutine exp_matrix(X,n,exp_X) if(dabs(eigvalues(i)).gt.1.d-4)then eigvalues_mat(i,i) = dsin(eigvalues(i))/eigvalues(i) else ! Taylor development of sin(x)/x near x=0 = 1 - x^2/6 - eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0 + eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0 & + + eigvalues(i)*eigvalues(i)*eigvalues(i)*eigvalues(i)*c_1_3*0.025d0 endif enddo ! matrix_tmp1 = W t^-1 sin(t) From ce43b16fc0c7f76023b1744bdfd4dcd23a8aee50 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 28 Feb 2024 14:46:23 +0100 Subject: [PATCH 5/8] Fixed bug in PT2 with fast stochastic convergence --- src/cipsi/run_pt2_slave.irp.f | 8 +++++++- src/tools/diagonalize_h.irp.f | 1 + src/tools/print_energy.irp.f | 3 ++- 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index debae596..cb1dd1f5 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -186,6 +186,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy) type(pt2_type) :: pt2_data integer :: n_tasks, k, N integer :: i_generator, subset + integer :: ifirst integer :: bsize ! Size of selection buffers logical :: sending @@ -202,6 +203,7 @@ subroutine run_pt2_slave_large(thread,iproc,energy) zmq_socket_push = new_zmq_push_socket(thread) + ifirst = 0 b%N = 0 buffer_ready = .False. n_tasks = 1 @@ -250,7 +252,11 @@ subroutine run_pt2_slave_large(thread,iproc,energy) call omp_set_lock(global_selection_buffer_lock) global_selection_buffer%mini = b%mini call merge_selection_buffers(b,global_selection_buffer) - b%cur=0 + if (ifirst /= 0 ) then + b%cur=0 + else + ifirst = 1 + endif call omp_unset_lock(global_selection_buffer_lock) if ( iproc == 1 ) then call omp_set_lock(global_selection_buffer_lock) diff --git a/src/tools/diagonalize_h.irp.f b/src/tools/diagonalize_h.irp.f index c9ae2033..ffc53aa2 100644 --- a/src/tools/diagonalize_h.irp.f +++ b/src/tools/diagonalize_h.irp.f @@ -20,4 +20,5 @@ subroutine routine call diagonalize_CI print*,'N_det = ',N_det call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) + call print_mol_properties end diff --git a/src/tools/print_energy.irp.f b/src/tools/print_energy.irp.f index 4fe1572c..0e67828e 100644 --- a/src/tools/print_energy.irp.f +++ b/src/tools/print_energy.irp.f @@ -14,5 +14,6 @@ end subroutine run implicit none - print *, psi_energy + nuclear_repulsion + call print_mol_properties + print *, psi_energy + nuclear_repulsion end From 30096e07ea33720d6104c13f43f9eb5788967ddc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 28 Feb 2024 15:02:39 +0100 Subject: [PATCH 6/8] Merging TC with cipsi --- .../cipsi_tc_bi_ortho/run_pt2_slave.irp.f | 45 +++++++++---------- .../run_selection_slave.irp.f | 29 ++++++------ .../cipsi_tc_bi_ortho/zmq_selection.irp.f | 3 +- src/cipsi/run_selection_slave.irp.f | 2 +- 4 files changed, 39 insertions(+), 40 deletions(-) diff --git a/plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f b/plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f index aa6546e7..d4f45649 100644 --- a/plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f @@ -31,11 +31,12 @@ subroutine run_pt2_slave(thread,iproc,energy) double precision, intent(in) :: energy(N_states_diag) integer, intent(in) :: thread, iproc - if (N_det > 100000 ) then - call run_pt2_slave_large(thread,iproc,energy) - else - call run_pt2_slave_small(thread,iproc,energy) - endif + call run_pt2_slave_large(thread,iproc,energy) +! if (N_det > 100000 ) then +! call run_pt2_slave_large(thread,iproc,energy) +! else +! call run_pt2_slave_small(thread,iproc,energy) +! endif end subroutine run_pt2_slave_small(thread,iproc,energy) @@ -178,15 +179,12 @@ subroutine run_pt2_slave_large(thread,iproc,energy) type(pt2_type) :: pt2_data integer :: n_tasks, k, N integer :: i_generator, subset + integer :: ifirst integer :: bsize ! Size of selection buffers logical :: sending - double precision :: time_shift - PROVIDE global_selection_buffer global_selection_buffer_lock - call random_number(time_shift) - time_shift = time_shift*15.d0 zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() @@ -198,15 +196,13 @@ subroutine run_pt2_slave_large(thread,iproc,energy) zmq_socket_push = new_zmq_push_socket(thread) + ifirst = 0 b%N = 0 buffer_ready = .False. n_tasks = 1 sending = .False. done = .False. - double precision :: time0, time1 - call wall_time(time0) - time0 = time0+time_shift do while (.not.done) integer, external :: get_tasks_from_taskserver @@ -233,28 +229,29 @@ subroutine run_pt2_slave_large(thread,iproc,energy) ASSERT (b%N == bsize) endif + double precision :: time0, time1 + call wall_time(time0) call pt2_alloc(pt2_data,N_states) b%cur = 0 call select_connected(i_generator,energy,pt2_data,b,subset,pt2_F(i_generator)) + call wall_time(time1) integer, external :: tasks_done_to_taskserver if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then done = .true. endif call sort_selection_buffer(b) - - call wall_time(time1) -! if (time1-time0 > 15.d0) then - call omp_set_lock(global_selection_buffer_lock) - global_selection_buffer%mini = b%mini - call merge_selection_buffers(b,global_selection_buffer) - b%cur=0 - call omp_unset_lock(global_selection_buffer_lock) - call wall_time(time0) -! endif - call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) - if ( iproc == 1 .or. i_generator < 100 .or. done) then + call omp_set_lock(global_selection_buffer_lock) + global_selection_buffer%mini = b%mini + call merge_selection_buffers(b,global_selection_buffer) + if (ifirst /= 0 ) then + b%cur=0 + else + ifirst = 1 + endif + call omp_unset_lock(global_selection_buffer_lock) + if ( iproc == 1 ) then call omp_set_lock(global_selection_buffer_lock) call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), global_selection_buffer, (/task_id/), 1,sending) global_selection_buffer%cur = 0 diff --git a/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f b/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f index d351cc79..39c83c4b 100644 --- a/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f @@ -5,19 +5,22 @@ subroutine run_selection_slave(thread, iproc, energy) implicit none - double precision, intent(in) :: energy(N_states) - integer, intent(in) :: thread, iproc + double precision, intent(in) :: energy(N_states) + integer, intent(in) :: thread, iproc + integer :: rc, i - integer :: rc, i - integer :: worker_id, task_id(1), ctask, ltask - character*(512) :: task - integer(ZMQ_PTR) :: zmq_to_qp_run_socket - integer(ZMQ_PTR) :: zmq_socket_push - integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket - integer(ZMQ_PTR), external :: new_zmq_push_socket - type(selection_buffer) :: buf, buf2 - type(pt2_type) :: pt2_data - logical :: done, buffer_ready + integer :: worker_id, task_id(1), ctask, ltask + character*(512) :: task + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_push_socket + integer(ZMQ_PTR) :: zmq_socket_push + + type(selection_buffer) :: buf, buf2 + logical :: done, buffer_ready + type(pt2_type) :: pt2_data PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order @@ -64,7 +67,7 @@ subroutine run_selection_slave(thread, iproc, energy) stop '-1' end if end if - call select_connected(i_generator, energy, pt2_data, buf,subset, pt2_F(i_generator)) + call select_connected(i_generator, energy, pt2_data, buf, subset, pt2_F(i_generator)) endif integer, external :: task_done_to_taskserver diff --git a/plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f b/plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f index dc3e0f27..22db643f 100644 --- a/plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/zmq_selection.irp.f @@ -11,7 +11,7 @@ subroutine ZMQ_selection(N_in, pt2_data) integer, external :: omp_get_thread_num type(pt2_type), intent(inout) :: pt2_data - PROVIDE psi_det psi_coef N_det qp_max_mem N_states pt2_F s2_eig N_det_generators +! PROVIDE psi_det psi_coef N_det qp_max_mem N_states pt2_F s2_eig N_det_generators N = max(N_in,1) N = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2) @@ -61,7 +61,6 @@ subroutine ZMQ_selection(N_in, pt2_data) ipos=1 task = ' ' - do i= 1, N_det_generators do j=1,pt2_F(i) write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, i, N diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index 91bd3a38..87ebca40 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -65,7 +65,7 @@ subroutine run_selection_slave(thread,iproc,energy) stop '-1' end if end if - call select_connected(i_generator,energy,pt2_data,buf,subset,pt2_F(i_generator)) + call select_connected(i_generator, energy, pt2_data, buf, subset, pt2_F(i_generator)) endif integer, external :: task_done_to_taskserver From b0fdf35983c6458805923c1a440cff72ad7fc150 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Fri, 1 Mar 2024 13:37:46 +0100 Subject: [PATCH 7/8] PROVIDING bug --- plugins/local/tc_bi_ortho/slater_tc_opt.irp.f | 23 +- .../tc_bi_ortho/slater_tc_opt_diag.irp.f | 2 + .../tc_bi_ortho/slater_tc_opt_single.irp.f | 419 ++++++++++-------- plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f | 13 + plugins/local/tc_bi_ortho/tc_hmat.irp.f | 28 +- 5 files changed, 281 insertions(+), 204 deletions(-) diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f index 9901a853..59efc943 100644 --- a/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f +++ b/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f @@ -8,8 +8,13 @@ subroutine provide_all_three_ints_bi_ortho() END_DOC implicit none + double precision :: t1, t2 + PROVIDE ao_two_e_integrals_in_map + print *, ' start provide_all_three_ints_bi_ortho' + call wall_time(t1) + if(three_body_h_tc) then if(three_e_3_idx_term) then @@ -32,6 +37,9 @@ subroutine provide_all_three_ints_bi_ortho() endif + call wall_time(t2) + print *, ' end provide_all_three_ints_bi_ortho after (min) = ', (t2-t1)/60.d0 + return end @@ -83,8 +91,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: degree + PROVIDE pure_three_body_h_tc + hmono = 0.d0 htwoe = 0.d0 htot = 0.d0 @@ -99,7 +110,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, if(degree == 0) then call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) else if (degree == 1) then - call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot) + call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) else if(degree == 2) then call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) endif @@ -111,7 +122,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, if(degree == 0) then call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) else if (degree == 1) then - call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot) + call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) else if(degree == 2) then call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) else @@ -149,16 +160,16 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot) double precision, intent(out) :: htot integer :: degree - htot = 0.d0 + htot = 0.d0 call get_excitation_degree(key_i, key_j, degree, Nint) if(degree.gt.2) return - if(degree == 0)then + if(degree == 0) then call diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_i,htot) - else if (degree == 1)then + else if (degree == 1) then call single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint,key_j, key_i , htot) - else if(degree == 2)then + else if(degree == 2) then call double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) endif diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f index cc1a0603..78f9dc66 100644 --- a/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f +++ b/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f @@ -15,6 +15,8 @@ implicit none double precision :: hmono, htwoe, htot, hthree + PROVIDE N_int + PROVIDE HF_bitmask PROVIDE mo_l_coef mo_r_coef call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot) diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f index 81bf69f4..e57cb05c 100644 --- a/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f +++ b/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -19,6 +19,7 @@ subroutine single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: occ(Nint*bit_kind_size,2) integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk integer :: degree,exc(0:2,2,2) @@ -44,27 +45,28 @@ subroutine single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, call bitstring_to_list_ab(key_i, occ, Ne, Nint) call get_single_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2) - call get_single_excitation_from_fock_tc(key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot) + call get_single_excitation_from_fock_tc(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot) end ! --- -subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot) +subroutine get_single_excitation_from_fock_tc(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot) use bitmasks implicit none + integer, intent(in) :: Nint integer, intent(in) :: h, p, spin double precision, intent(in) :: phase - integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2) + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) double precision, intent(out) :: hmono, htwoe, hthree, htot - integer(bit_kind) :: differences(N_int,2) - integer(bit_kind) :: hole(N_int,2) - integer(bit_kind) :: partcl(N_int,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_partcl(N_int*bit_kind_size,2) + integer(bit_kind) :: differences(Nint,2) + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: partcl(Nint,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer :: occ_partcl(Nint*bit_kind_size,2) integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) integer :: i0,i double precision :: buffer_c(mo_num),buffer_x(mo_num) @@ -74,7 +76,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) enddo - do i = 1, N_int + do i = 1, Nint differences(i,1) = xor(key_i(i,1), ref_closed_shell_bitmask(i,1)) differences(i,2) = xor(key_i(i,2), ref_closed_shell_bitmask(i,2)) hole (i,1) = iand(differences(i,1), ref_closed_shell_bitmask(i,1)) @@ -83,8 +85,8 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h partcl (i,2) = iand(differences(i,2), key_i(i,2)) enddo - call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) - call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint) hmono = mo_bi_ortho_tc_one_e(p,h) htwoe = fock_op_2_e_tc_closed_shell(p,h) @@ -122,7 +124,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h hthree = 0.d0 if (three_body_h_tc .and. elec_num.gt.2 .and. three_e_4_idx_term) then - call three_comp_fock_elem(key_i, h, p, spin, hthree) + call three_comp_fock_elem(Nint, key_i, h, p, spin, hthree) endif htwoe = htwoe * phase @@ -134,24 +136,27 @@ end ! --- -subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) - implicit none - integer,intent(in) :: h_fock,p_fock,ispin_fock - integer(bit_kind), intent(in) :: key_i(N_int,2) - double precision, intent(out) :: hthree - integer :: nexc(2),i,ispin,na,nb - integer(bit_kind) :: hole(N_int,2) - integer(bit_kind) :: particle(N_int,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_particle(N_int*bit_kind_size,2) - integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) - integer(bit_kind) :: det_tmp(N_int,2) +subroutine three_comp_fock_elem(Nint, key_i, h_fock, p_fock, ispin_fock, hthree) + implicit none + integer, intent(in) :: Nint + integer, intent(in) :: h_fock, p_fock, ispin_fock + integer(bit_kind), intent(in) :: key_i(Nint,2) + double precision, intent(out) :: hthree + + integer :: nexc(2),i,ispin,na,nb + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: particle(Nint,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer :: occ_particle(Nint*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) + integer(bit_kind) :: det_tmp(Nint,2) nexc(1) = 0 nexc(2) = 0 + !! Get all the holes and particles of key_i with respect to the ROHF determinant - do i=1,N_int + do i = 1, Nint hole(i,1) = xor(key_i(i,1),ref_bitmask(i,1)) hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2)) particle(i,1) = iand(hole(i,1),key_i(i,1)) @@ -161,13 +166,14 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) nexc(1) = nexc(1) + popcnt(hole(i,1)) nexc(2) = nexc(2) + popcnt(hole(i,2)) enddo + integer :: tmp(2) !DIR$ FORCEINLINE - call bitstring_to_list_ab(particle, occ_particle, tmp, N_int) + call bitstring_to_list_ab(particle, occ_particle, tmp, Nint) ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha ASSERT (tmp(2) == nexc(2)) ! Number of particle beta !DIR$ FORCEINLINE - call bitstring_to_list_ab(hole, occ_hole, tmp, N_int) + call bitstring_to_list_ab(hole, occ_hole, tmp, Nint) ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha ASSERT (tmp(2) == nexc(2)) ! Number of holes beta @@ -181,15 +187,18 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) do ispin=1,2 na = elec_num_tab(ispin) nb = elec_num_tab(iand(ispin,1)+1) - do i=1,nexc(ispin) + do i = 1, nexc(ispin) !DIR$ FORCEINLINE - call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb) + call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb) !DIR$ FORCEINLINE - call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb) + call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb) enddo enddo + end +! --- + subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb) use bitmasks implicit none @@ -365,111 +374,118 @@ subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,N end +! --- -BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ] - implicit none - BEGIN_DOC -! Closed-shell part of the Fock operator for the TC operator - END_DOC - integer :: h0,p0,h,p,k0,k,i - integer :: n_occ_ab(2) - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab_virt(2) - integer :: occ_virt(N_int*bit_kind_size,2) - integer(bit_kind) :: key_test(N_int) - integer(bit_kind) :: key_virt(N_int,2) - double precision :: accu +BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)] - fock_op_2_e_tc_closed_shell = -1000.d0 - call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) - do i = 1, N_int - key_virt(i,1) = full_ijkl_bitmask(i) - key_virt(i,2) = full_ijkl_bitmask(i) - key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) - key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) - enddo - call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) - ! docc ---> virt single excitations - do h0 = 1, n_occ_ab(1) - h=occ(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p = occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu + BEGIN_DOC + ! Closed-shell part of the Fock operator for the TC operator + END_DOC + + implicit none + + PROVIDE N_int + + integer :: h0,p0,h,p,k0,k,i + integer :: n_occ_ab(2) + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab_virt(2) + integer :: occ_virt(N_int*bit_kind_size,2) + integer(bit_kind) :: key_test(N_int) + integer(bit_kind) :: key_virt(N_int,2) + double precision :: accu + + fock_op_2_e_tc_closed_shell = -1000.d0 + call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) + + do i = 1, N_int + key_virt(i,1) = full_ijkl_bitmask(i) + key_virt(i,2) = full_ijkl_bitmask(i) + key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) + key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) enddo - enddo - - do h0 = 1, n_occ_ab_virt(1) - h = occ_virt(h0,1) - do p0 = 1, n_occ_ab(1) - p=occ(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu + call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) + ! docc ---> virt single excitations + do h0 = 1, n_occ_ab(1) + h = occ(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo enddo - enddo - - ! virt ---> virt single excitations - do h0 = 1, n_occ_ab_virt(1) - h=occ_virt(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p = occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu + + do h0 = 1, n_occ_ab_virt(1) + h = occ_virt(h0,1) + do p0 = 1, n_occ_ab(1) + p = occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo enddo - enddo - - do h0 = 1, n_occ_ab_virt(1) - h = occ_virt(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p=occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + + ! virt ---> virt single excitations + do h0 = 1, n_occ_ab_virt(1) + h=occ_virt(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu enddo - fock_op_2_e_tc_closed_shell(p,h) = accu enddo - enddo - - - ! docc ---> docc single excitations - do h0 = 1, n_occ_ab(1) - h=occ(h0,1) - do p0 = 1, n_occ_ab(1) - p = occ(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + + do h0 = 1, n_occ_ab_virt(1) + h = occ_virt(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p=occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu enddo - fock_op_2_e_tc_closed_shell(p,h) = accu enddo - enddo - - do h0 = 1, n_occ_ab(1) - h = occ(h0,1) - do p0 = 1, n_occ_ab(1) - p=occ(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + + + ! docc ---> docc single excitations + do h0 = 1, n_occ_ab(1) + h=occ(h0,1) + do p0 = 1, n_occ_ab(1) + p = occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + do h0 = 1, n_occ_ab(1) + h = occ(h0,1) + do p0 = 1, n_occ_ab(1) + p=occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu enddo - fock_op_2_e_tc_closed_shell(p,h) = accu enddo - enddo ! do i = 1, mo_num ! write(*,'(100(F10.5,X))')fock_op_2_e_tc_closed_shell(:,i) @@ -477,8 +493,10 @@ BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) END_PROVIDER +! --- subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) + BEGIN_DOC ! for single excitation ONLY FOR ONE- AND TWO-BODY TERMS !! @@ -492,8 +510,9 @@ subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) implicit none integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) - double precision, intent(out) :: htot - double precision :: hmono, htwoe + double precision, intent(out) :: htot + + double precision :: hmono, htwoe integer :: occ(Nint*bit_kind_size,2) integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk integer :: degree,exc(0:2,2,2) @@ -517,75 +536,85 @@ subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) call get_single_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) - call get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,htot) -end - - -subroutine get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h,p,spin,phase,hmono,htwoe,htot) - use bitmasks - implicit none - integer,intent(in) :: h,p,spin - double precision, intent(in) :: phase - integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2) - double precision, intent(out) :: hmono,htwoe,htot - integer(bit_kind) :: differences(N_int,2) - integer(bit_kind) :: hole(N_int,2) - integer(bit_kind) :: partcl(N_int,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_partcl(N_int*bit_kind_size,2) - integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) - integer :: i0,i - double precision :: buffer_c(mo_num),buffer_x(mo_num) - do i=1, mo_num - buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h) - buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) - enddo - do i = 1, N_int - differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1)) - differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2)) - hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) - hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) - partcl(i,1) = iand(differences(i,1),key_i(i,1)) - partcl(i,2) = iand(differences(i,2),key_i(i,2)) - enddo - call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) - call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) - hmono = mo_bi_ortho_tc_one_e(p,h) - htwoe = fock_op_2_e_tc_closed_shell(p,h) - ! holes :: direct terms - do i0 = 1, n_occ_ab_hole(1) - i = occ_hole(i0,1) - htwoe -= buffer_c(i) - enddo - do i0 = 1, n_occ_ab_hole(2) - i = occ_hole(i0,2) - htwoe -= buffer_c(i) - enddo - - ! holes :: exchange terms - do i0 = 1, n_occ_ab_hole(spin) - i = occ_hole(i0,spin) - htwoe += buffer_x(i) - enddo - - ! particles :: direct terms - do i0 = 1, n_occ_ab_partcl(1) - i = occ_partcl(i0,1) - htwoe += buffer_c(i) - enddo - do i0 = 1, n_occ_ab_partcl(2) - i = occ_partcl(i0,2) - htwoe += buffer_c(i) - enddo - - ! particles :: exchange terms - do i0 = 1, n_occ_ab_partcl(spin) - i = occ_partcl(i0,spin) - htwoe -= buffer_x(i) - enddo - htwoe = htwoe * phase - hmono = hmono * phase - htot = htwoe + hmono + call get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, htot) + +end + +! --- + +subroutine get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, htot) + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer, intent(in) :: h, p, spin + double precision, intent(in) :: phase + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hmono,htwoe,htot + + integer(bit_kind) :: differences(Nint,2) + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: partcl(Nint,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer :: occ_partcl(Nint*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) + integer :: i0,i + double precision :: buffer_c(mo_num), buffer_x(mo_num) + + do i = 1, mo_num + buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h) + buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) + enddo + + do i = 1, Nint + differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1)) + differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2)) + hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) + hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) + partcl(i,1) = iand(differences(i,1),key_i(i,1)) + partcl(i,2) = iand(differences(i,2),key_i(i,2)) + enddo + + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint) + hmono = mo_bi_ortho_tc_one_e(p,h) + htwoe = fock_op_2_e_tc_closed_shell(p,h) + + ! holes :: direct terms + do i0 = 1, n_occ_ab_hole(1) + i = occ_hole(i0,1) + htwoe -= buffer_c(i) + enddo + do i0 = 1, n_occ_ab_hole(2) + i = occ_hole(i0,2) + htwoe -= buffer_c(i) + enddo + + ! holes :: exchange terms + do i0 = 1, n_occ_ab_hole(spin) + i = occ_hole(i0,spin) + htwoe += buffer_x(i) + enddo + + ! particles :: direct terms + do i0 = 1, n_occ_ab_partcl(1) + i = occ_partcl(i0,1) + htwoe += buffer_c(i) + enddo + do i0 = 1, n_occ_ab_partcl(2) + i = occ_partcl(i0,2) + htwoe += buffer_c(i) + enddo + + ! particles :: exchange terms + do i0 = 1, n_occ_ab_partcl(spin) + i = occ_partcl(i0,spin) + htwoe -= buffer_x(i) + enddo + htwoe = htwoe * phase + hmono = hmono * phase + htot = htwoe + hmono end diff --git a/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f b/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f index 64982ab6..398e96db 100644 --- a/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f +++ b/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f @@ -7,6 +7,10 @@ program tc_bi_ortho ! END_DOC + implicit none + + PROVIDE N_int + my_grid_becke = .True. PROVIDE tc_grid1_a tc_grid1_r my_n_pt_r_grid = tc_grid1_r @@ -66,6 +70,15 @@ subroutine routine_diag() ! provide overlap_bi_ortho ! provide htilde_matrix_elmt_bi_ortho + if(noL_standard) then + PROVIDE noL_0e + PROVIDE noL_1e + PROVIDE noL_2e + endif + + PROVIDE htilde_matrix_elmt_bi_ortho + return + if(N_states .eq. 1) then print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) diff --git a/plugins/local/tc_bi_ortho/tc_hmat.irp.f b/plugins/local/tc_bi_ortho/tc_hmat.irp.f index 88652caa..abec410d 100644 --- a/plugins/local/tc_bi_ortho/tc_hmat.irp.f +++ b/plugins/local/tc_bi_ortho/tc_hmat.irp.f @@ -13,16 +13,34 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)] implicit none integer :: i, j + double precision :: t1, t2 double precision :: htot + + PROVIDE N_int + PROVIDE psi_det + PROVIDE three_e_3_idx_term - call provide_all_three_ints_bi_ortho + if(noL_standard) then + PROVIDE noL_0e + PROVIDE noL_1e + PROVIDE noL_2e + endif + + print *, ' PROVIDING htilde_matrix_elmt_bi_ortho ...' + call wall_time(t1) + + call provide_all_three_ints_bi_ortho() i = 1 j = 1 call htilde_mu_mat_opt_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot) - !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j, htot) & + + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, htot) & !$OMP SHARED (N_det, psi_det, N_int, htilde_matrix_elmt_bi_ortho) + !$OMP DO do i = 1, N_det do j = 1, N_det ! < J |Htilde | I > @@ -31,7 +49,11 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)] htilde_matrix_elmt_bi_ortho(j,i) = htot enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(t2) + print *, ' wall time for htilde_matrix_elmt_bi_ortho (min) =', (t2-t1)/60.d0 END_PROVIDER From eaab1b80648bd2b4ec38aa17722a751fccf42ce0 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 7 Mar 2024 07:34:59 +0100 Subject: [PATCH 8/8] few modif for HHG --- .../multi_s_dipole_moment.irp.f | 23 ++++++++++++++++ src/utils/constants.include.F | 27 +++++++++++++++++++ 2 files changed, 50 insertions(+) diff --git a/src/mol_properties/multi_s_dipole_moment.irp.f b/src/mol_properties/multi_s_dipole_moment.irp.f index 913ae2f3..f21e08cd 100644 --- a/src/mol_properties/multi_s_dipole_moment.irp.f +++ b/src/mol_properties/multi_s_dipole_moment.irp.f @@ -91,3 +91,26 @@ BEGIN_PROVIDER [double precision, multi_s_dipole_moment, (N_states, N_states)] enddo END_PROVIDER + +! --- + + BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment_eigenvec, (N_states, N_states)] +&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment_eigenvec, (N_states, N_states)] +&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment_eigenvec, (N_states, N_states)] +&BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment_eigenval, (N_states)] +&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment_eigenval, (N_states)] +&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment_eigenval, (N_states)] + + implicit none + + PROVIDE multi_s_x_dipole_moment multi_s_y_dipole_moment multi_s_z_dipole_moment + + call lapack_diag(multi_s_x_dipole_moment_eigenval(1), multi_s_x_dipole_moment_eigenvec(1,1), multi_s_x_dipole_moment(1,1), N_states, N_states) + call lapack_diag(multi_s_y_dipole_moment_eigenval(1), multi_s_y_dipole_moment_eigenvec(1,1), multi_s_y_dipole_moment(1,1), N_states, N_states) + call lapack_diag(multi_s_z_dipole_moment_eigenval(1), multi_s_z_dipole_moment_eigenvec(1,1), multi_s_z_dipole_moment(1,1), N_states, N_states) + +END_PROVIDER + +! --- + + diff --git a/src/utils/constants.include.F b/src/utils/constants.include.F index 422eff95..7b01f888 100644 --- a/src/utils/constants.include.F +++ b/src/utils/constants.include.F @@ -18,3 +18,30 @@ double precision, parameter :: c_4_3 = 4.d0/3.d0 double precision, parameter :: c_1_3 = 1.d0/3.d0 double precision, parameter :: sq_op5 = dsqrt(0.5d0) double precision, parameter :: dlog_2pi = dlog(2.d0*dacos(-1.d0)) + +! physical constants and units conversion factors +double precision, parameter :: k_boltzman_si = 1.38066d-23 ! K k^-1 +double precision, parameter :: k_boltzman_au = 3.1667d-6 ! Hartree k^-1 +double precision, parameter :: k_boltzman_m1_au = 315795.26d0 ! Hartree^-1 k +double precision, parameter :: bohr_radius_si = 0.529177d-10 ! m +double precision, parameter :: bohr_radius_cm = 0.529177d-8 ! cm +double precision, parameter :: bohr_radius_angs = 0.529177d0 ! Angstrom +double precision, parameter :: electronmass_si = 9.10953d-31 ! Kg +double precision, parameter :: electronmass_uma = 5.4858d-4 ! uma +double precision, parameter :: electronvolt_si = 1.6021892d-19 ! J +double precision, parameter :: uma_si = 1.66057d-27 ! Kg +double precision, parameter :: debye_si = 3.33564d-30 ! coulomb meter +double precision, parameter :: debye_au = 0.393427228d0 ! e * Bohr +double precision, parameter :: angstrom_to_au = 1.889727d0 ! au +double precision, parameter :: au_to_ohmcmm1 = 46000.0d0 ! (ohm cm)^-1 +double precision, parameter :: au_to_kb = 294210.0d0 ! kbar +double precision, parameter :: au_to_eV = 27.211652d0 +double precision, parameter :: uma_to_au = 1822.89d0 +double precision, parameter :: au_to_terahertz = 2.4189d-5 +double precision, parameter :: au_to_sec = 2.4189d-17 +double precision, parameter :: au_to_fsec = 2.4189d-2 +double precision, parameter :: Wcm2 = 3.5d16 +double precision, parameter :: amconv = 1.66042d-24/9.1095d-28*0.5d0 ! mass conversion: a.m.u to a.u. (ry) +double precision, parameter :: uakbar = 147105.d0 ! pressure conversion from ry/(a.u)^3 to k + +