From 8c4e5a3b159105323da2c5003e07f9fa8efc25d3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 1 Mar 2021 13:15:21 +0100 Subject: [PATCH 01/11] HF --- src/scf_utils/roothaan_hall_scf.irp.f | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 1406d603..5f445ad6 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -23,6 +23,10 @@ END_DOC error_matrix_DIIS(ao_num,ao_num,max_dim_DIIS) & ) + Fock_matrix_DIIS = 0.d0 + error_matrix_DIIS = 0.d0 + mo_coef_save = 0.d0 + call write_time(6) print*,'Energy of the guess = ',SCF_energy @@ -308,6 +312,7 @@ END_DOC do k=1,dim_DIIS if (dabs(X_vector_DIIS(k)) < 1.d-10) cycle do i=1,ao_num + ! FPE here Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + & X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1) enddo From ee9f4516b39d8dd8fa2476d3d5d738bacfd91718 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 2 Mar 2021 14:27:22 +0100 Subject: [PATCH 02/11] Fixing zmq bug --- src/zmq/utils.irp.f | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/zmq/utils.irp.f b/src/zmq/utils.irp.f index 2cdb30f1..058f6bca 100644 --- a/src/zmq/utils.irp.f +++ b/src/zmq/utils.irp.f @@ -963,7 +963,7 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i get_tasks_from_taskserver = 0 - write(message,*) 'get_tasks '//trim(zmq_state), worker_id, n_tasks + write(message,'(A,A,X,I10,I10)') 'get_tasks ', trim(zmq_state), worker_id, n_tasks sze = len(trim(message)) rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0) @@ -974,8 +974,10 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i message = repeat(' ',1024) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 1024, 0) - rc = min(1024,rc) - read(message(1:rc),*, end=10, err=10) reply + if (rc <= 0) then + get_tasks_from_taskserver = -1 + return + endif if (trim(message) == 'get_tasks_reply ok') then continue else if (trim(message) == 'terminate') then @@ -993,6 +995,10 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i do i=1,n_tasks message = repeat(' ',512) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 1024, 0) + if (rc <= 0) then + get_tasks_from_taskserver = -1 + return + endif rc = min(1024,rc) read(message(1:rc),*, end=10, err=10) task_id(i) if (task_id(i) == 0) then @@ -1001,10 +1007,10 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i exit endif rc = 1 - do while (message(rc:rc) == ' ') + do while (rc < 1024 .and. message(rc:rc) == ' ') rc += 1 enddo - do while (message(rc:rc) /= ' ') + do while (rc < 1024 .and. message(rc:rc) /= ' ') rc += 1 enddo rc += 1 From d55acad33fe3d7eee6005173afbb1db3b6bed78b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Mar 2021 14:22:37 +0100 Subject: [PATCH 03/11] libcap without gperf --- configure | 1 + 1 file changed, 1 insertion(+) diff --git a/configure b/configure index 391f5384..31f83002 100755 --- a/configure +++ b/configure @@ -238,6 +238,7 @@ EOF tar --gunzip --extract --file libcap.tar.gz rm libcap.tar.gz cd libcap-*/libcap + prefix=$QP_ROOT make BUILD_GPERF=no prefix=$QP_ROOT make install EOF From 533413277df93882b265dcfee80dd0c879fb5d08 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 14 Mar 2021 01:14:12 +0100 Subject: [PATCH 04/11] Simplified 2_RDM --- src/two_body_rdm/two_e_dm_mo.irp.f | 519 +---------------------------- 1 file changed, 14 insertions(+), 505 deletions(-) diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f index 203c7ff2..47f5983f 100644 --- a/src/two_body_rdm/two_e_dm_mo.irp.f +++ b/src/two_body_rdm/two_e_dm_mo.irp.f @@ -1,317 +1,4 @@ -BEGIN_PROVIDER [double precision, two_e_dm_ab_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] - implicit none - two_e_dm_ab_mo = 0.d0 - integer :: i,j,k,l,iorb,jorb,korb,lorb,istate - BEGIN_DOC - ! two_e_dm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons - ! - ! - ! - ! WHERE ALL ORBITALS (i,j,k,l) BELONG TO ALL OCCUPIED ORBITALS : core, inactive and active - ! - ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} - ! - ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" - ! - ! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA - ! - ! two_e_dm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta - ! - ! Therefore you don't necessary have symmetry between electron 1 and 2 - ! - ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO ARE SET TO ZERO - END_DOC - - two_e_dm_ab_mo = 0.d0 - do istate = 1, N_states - !! PURE ACTIVE PART ALPHA-BETA - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_act_orb - korb = list_act(k) - do l = 1, n_act_orb - lorb = list_act(l) - ! alph beta alph beta - two_e_dm_ab_mo(lorb,korb,jorb,iorb,istate) = act_2_rdm_ab_mo(l,k,j,i,istate) - enddo - enddo - enddo - enddo - !! BETA ACTIVE - ALPHA inactive - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! alph beta alph beta - two_e_dm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) - enddo - enddo - enddo - - !! ALPHA ACTIVE - BETA inactive - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! alph beta alph beta - two_e_dm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - - !! ALPHA INACTIVE - BETA INACTIVE - !! - do j = 1, n_inact_orb - jorb = list_inact(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! alph beta alph beta - two_e_dm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0 - enddo - enddo - - !! BETA ACTIVE - ALPHA CORE - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_core_orb - korb = list_core(k) - ! alph beta alph beta - two_e_dm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) - enddo - enddo - enddo - - !! ALPHA ACTIVE - BETA CORE - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_core_orb - korb = list_core(k) - ! alph beta alph beta - two_e_dm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - - !! ALPHA CORE - BETA CORE - !! - do j = 1, n_core_orb - jorb = list_core(j) - do k = 1, n_core_orb - korb = list_core(k) - ! alph beta alph beta - two_e_dm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0 - enddo - enddo - - enddo -END_PROVIDER - - -BEGIN_PROVIDER [double precision, two_e_dm_aa_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] - implicit none - BEGIN_DOC - ! two_e_dm_aa_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/alpha electrons - ! - ! - ! - ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active - ! - ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2 - ! - ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" - ! - ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero - END_DOC - two_e_dm_aa_mo = 0.d0 - integer :: i,j,k,l,iorb,jorb,korb,lorb,istate - - do istate = 1, N_states - !! PURE ACTIVE PART ALPHA-ALPHA - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_act_orb - korb = list_act(k) - do l = 1, n_act_orb - lorb = list_act(l) - two_e_dm_aa_mo(lorb,korb,jorb,iorb,istate) = & - act_2_rdm_aa_mo(l,k,j,i,istate) - enddo - enddo - enddo - enddo - !! ALPHA ACTIVE - ALPHA inactive - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - - !! ALPHA INACTIVE - ALPHA INACTIVE - do j = 1, n_inact_orb - jorb = list_inact(j) - do k = 1, n_inact_orb - korb = list_inact(k) - two_e_dm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - - - !! ALPHA ACTIVE - ALPHA CORE - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_core_orb - korb = list_core(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - !! ALPHA CORE - ALPHA CORE - - do j = 1, n_core_orb - jorb = list_core(j) - do k = 1, n_core_orb - korb = list_core(k) - two_e_dm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - - enddo - -END_PROVIDER - -BEGIN_PROVIDER [double precision, two_e_dm_bb_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] - implicit none - BEGIN_DOC - ! two_e_dm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons - ! - ! - ! - ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active - ! - ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2 - ! - ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" - ! - ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero - END_DOC - - integer :: i,j,k,l,iorb,jorb,korb,lorb,istate - two_e_dm_bb_mo = 0.d0 - - do istate = 1, N_states - !! PURE ACTIVE PART beta-beta - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_act_orb - korb = list_act(k) - do l = 1, n_act_orb - lorb = list_act(l) - two_e_dm_bb_mo(lorb,korb,jorb,iorb,istate) = & - act_2_rdm_bb_mo(l,k,j,i,istate) - enddo - enddo - enddo - enddo - !! beta ACTIVE - beta inactive - !! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - enddo - enddo - enddo - - !! beta INACTIVE - beta INACTIVE - do j = 1, n_inact_orb - jorb = list_inact(j) - do k = 1, n_inact_orb - korb = list_inact(k) - two_e_dm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - - !! beta ACTIVE - beta CORE - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_core_orb - korb = list_core(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - enddo - enddo - enddo - !! beta CORE - beta CORE - - do j = 1, n_core_orb - jorb = list_core(j) - do k = 1, n_core_orb - korb = list_core(k) - two_e_dm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - - enddo - -END_PROVIDER - -BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] +BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num,1)] implicit none BEGIN_DOC ! two_e_dm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons @@ -334,197 +21,19 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num,N_st two_e_dm_mo = 0.d0 integer :: i,j,k,l,iorb,jorb,korb,lorb,istate - do istate = 1, N_states - !!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!! - !! PURE ACTIVE PART SPIN-TRACE - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_act_orb - korb = list_act(k) - do l = 1, n_act_orb - lorb = list_act(l) - two_e_dm_mo(lorb,korb,jorb,iorb,istate) += & - act_2_rdm_spin_trace_mo(l,k,j,i,istate) - enddo - enddo - enddo + do l=1,mo_num + lorb = list_core_inact_act(l) + do k=1,mo_num + korb = list_core_inact_act(k) + do j=1,mo_num + jorb = list_core_inact_act(j) + do i=1,mo_num + iorb = list_core_inact_act(i) + two_e_dm_mo(iorb,jorb,korb,lorb,1) = state_av_full_occ_2_rdm_spin_trace_mo(i,j,k,l) + enddo enddo - - !!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!! - !!!!! BETA-BETA !!!!! - !! beta ACTIVE - beta inactive - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - enddo - enddo - enddo - !! beta INACTIVE - beta INACTIVE - do j = 1, n_inact_orb - jorb = list_inact(j) - do k = 1, n_inact_orb - korb = list_inact(k) - two_e_dm_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - if (.not.no_core_density)then - !! beta ACTIVE - beta CORE - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_core_orb - korb = list_core(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - two_e_dm_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - enddo - enddo - enddo - !! beta CORE - beta CORE - do j = 1, n_core_orb - jorb = list_core(j) - do k = 1, n_core_orb - korb = list_core(k) - two_e_dm_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - endif - - !!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!! - !!!!! ALPHA-ALPHA !!!!! - !! ALPHA ACTIVE - ALPHA inactive - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - !! ALPHA INACTIVE - ALPHA INACTIVE - do j = 1, n_inact_orb - jorb = list_inact(j) - do k = 1, n_inact_orb - korb = list_inact(k) - two_e_dm_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_core_orb - korb = list_core(k) - ! 1 2 1 2 : DIRECT TERM - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - ! 1 2 1 2 : EXCHANGE TERM - two_e_dm_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - two_e_dm_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - !! ALPHA CORE - ALPHA CORE - do j = 1, n_core_orb - jorb = list_core(j) - do k = 1, n_core_orb - korb = list_core(k) - two_e_dm_mo(korb,jorb,korb,jorb,istate) += 0.5d0 - two_e_dm_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 - enddo - enddo - - !!!!! ALPHA-BETA + BETA-ALPHA !!!!! - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! ALPHA INACTIVE - BETA ACTIVE - ! alph beta alph beta - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - ! beta alph beta alph - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) - ! BETA INACTIVE - ALPHA ACTIVE - ! beta alph beta alpha - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - ! alph beta alph beta - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - !! ALPHA INACTIVE - BETA INACTIVE - do j = 1, n_inact_orb - jorb = list_inact(j) - do k = 1, n_inact_orb - korb = list_inact(k) - ! alph beta alph beta - two_e_dm_mo(korb,jorb,korb,jorb,istate) += 0.5D0 - two_e_dm_mo(jorb,korb,jorb,korb,istate) += 0.5D0 - enddo - enddo - - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_act_orb - jorb = list_act(j) - do k = 1, n_core_orb - korb = list_core(k) - !! BETA ACTIVE - ALPHA CORE - ! alph beta alph beta - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate) - ! beta alph beta alph - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate) - !! ALPHA ACTIVE - BETA CORE - ! alph beta alph beta - two_e_dm_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate) - ! beta alph beta alph - two_e_dm_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate) - enddo - enddo - enddo - !! ALPHA CORE - BETA CORE - do j = 1, n_core_orb - jorb = list_core(j) - do k = 1, n_core_orb - korb = list_core(k) - ! alph beta alph beta - two_e_dm_mo(korb,jorb,korb,jorb,istate) += 0.5D0 - two_e_dm_mo(jorb,korb,jorb,korb,istate) += 0.5D0 - enddo - enddo - + enddo enddo - two_e_dm_mo(:,:,:,:,:) = two_e_dm_mo(:,:,:,:,:) * 2.d0 -END_PROVIDER + END_PROVIDER + From 2465b1b91d843f5588d594a4e3d845cfd74820b6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 21 Mar 2021 10:26:03 +0100 Subject: [PATCH 05/11] Comment to optimize DIIS --- src/scf_utils/roothaan_hall_scf.irp.f | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 5f445ad6..a63b46b6 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -218,7 +218,7 @@ END_DOC scratch(ao_num,ao_num) & ) -! Compute the matrices B and X + ! Compute the matrices B and X B_matrix_DIIS(:,:) = 0.d0 do j=1,dim_DIIS j_DIIS = min(dim_DIIS,mod(iteration_SCF-j,max_dim_DIIS)+1) @@ -226,7 +226,7 @@ END_DOC i_DIIS = min(dim_DIIS,mod(iteration_SCF-i,max_dim_DIIS)+1) -! Compute product of two errors vectors + ! Compute product of two errors vectors call dgemm('N','N',ao_num,ao_num,ao_num, & 1.d0, & @@ -235,7 +235,7 @@ END_DOC 0.d0, & scratch,size(scratch,1)) -! Compute Trace + ! Compute Trace do k=1,ao_num B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + scratch(k,k) @@ -243,6 +243,17 @@ END_DOC enddo enddo +! TODO : Could be simplified (to be checked) +! +! call dgemm('T','N', & +! min(dim_DIIS,iteration_SCF), min(dim_DIIS,iteration_SCF), & +! ao_num*ao_num, & +! 1.d0, & +! error_matrix_DIIS,size(error_matrix_DIIS,1)*size(error_matrix_DIIS,2), & +! error_matrix_DIIS,size(error_matrix_DIIS,1)*size(error_matrix_DIIS,2), & +! 0.d0, & +! B_matrix_DIIS,size(B_matrix_DIIS,1)) + ! Pad B matrix and build the X matrix C_vector_DIIS(:) = 0.d0 From 1d148f84bbf9f5ccdeee8168e4c82bb7722e82e2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 21 Mar 2021 21:06:41 +0100 Subject: [PATCH 06/11] erf -> derf --- src/utils/integration.irp.f | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index c907e425..974417b1 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -447,7 +447,7 @@ double precision function rint(n,rho) else u_inv=1.d0/dsqrt(rho) u=rho*u_inv - rint=0.5d0*u_inv*sqpi*erf(u) + rint=0.5d0*u_inv*sqpi*derf(u) endif return endif @@ -463,7 +463,7 @@ double precision function rint(n,rho) endif u=rho*u_inv two_rho_inv = 0.5d0*u_inv*u_inv - val0=0.5d0*u_inv*sqpi*erf(u) + val0=0.5d0*u_inv*sqpi*derf(u) rint=(val0-v)*two_rho_inv do k=2,n rint=(rint*dfloat(k+k-1)-v)*two_rho_inv @@ -496,7 +496,7 @@ double precision function rint_sum(n_pt_out,rho,d1) else u_inv=1.d0/dsqrt(rho) u=rho*u_inv - rint_sum=0.5d0*u_inv*sqpi*erf(u) *d1(0) + rint_sum=0.5d0*u_inv*sqpi*derf(u) *d1(0) endif do i=2,n_pt_out,2 @@ -515,7 +515,7 @@ double precision function rint_sum(n_pt_out,rho,d1) u_inv=1.d0/dsqrt(rho) u=rho*u_inv two_rho_inv = 0.5d0*u_inv*u_inv - val0=0.5d0*u_inv*sqpi*erf(u) + val0=0.5d0*u_inv*sqpi*derf(u) rint_sum=val0*d1(0) rint_tmp=(val0-v)*two_rho_inv di = 3.d0 From 610304c37ac24b2d374acb89e08592ac4893ab55 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 22 Mar 2021 09:25:36 +0100 Subject: [PATCH 07/11] Accelerate HF --- src/scf_utils/roothaan_hall_scf.irp.f | 33 +++++++-------------------- 1 file changed, 8 insertions(+), 25 deletions(-) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index a63b46b6..3b9eaeb4 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -202,7 +202,7 @@ END_DOC double precision,allocatable :: C_vector_DIIS(:) double precision,allocatable :: scratch(:,:) - integer :: i,j,k,i_DIIS,j_DIIS + integer :: i,j,k,l,i_DIIS,j_DIIS double precision :: rcond, ferr, berr integer, allocatable :: iwork(:) integer :: lwork @@ -222,38 +222,21 @@ END_DOC B_matrix_DIIS(:,:) = 0.d0 do j=1,dim_DIIS j_DIIS = min(dim_DIIS,mod(iteration_SCF-j,max_dim_DIIS)+1) - do i=1,dim_DIIS + do i=1,dim_DIIS i_DIIS = min(dim_DIIS,mod(iteration_SCF-i,max_dim_DIIS)+1) ! Compute product of two errors vectors - - call dgemm('N','N',ao_num,ao_num,ao_num, & - 1.d0, & - error_matrix_DIIS(1,1,i_DIIS),size(error_matrix_DIIS,1), & - error_matrix_DIIS(1,1,j_DIIS),size(error_matrix_DIIS,1), & - 0.d0, & - scratch,size(scratch,1)) - - ! Compute Trace - - do k=1,ao_num - B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + scratch(k,k) + do l=1,ao_num + do k=1,ao_num + B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + & + error_matrix_DIIS(k,l,i_DIIS) * error_matrix_DIIS(k,l,j_DIIS) + enddo enddo + enddo enddo -! TODO : Could be simplified (to be checked) -! -! call dgemm('T','N', & -! min(dim_DIIS,iteration_SCF), min(dim_DIIS,iteration_SCF), & -! ao_num*ao_num, & -! 1.d0, & -! error_matrix_DIIS,size(error_matrix_DIIS,1)*size(error_matrix_DIIS,2), & -! error_matrix_DIIS,size(error_matrix_DIIS,1)*size(error_matrix_DIIS,2), & -! 0.d0, & -! B_matrix_DIIS,size(B_matrix_DIIS,1)) - ! Pad B matrix and build the X matrix C_vector_DIIS(:) = 0.d0 From d217a6b9f1d897aef2375e3cfb9ef7fbc780733b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Mar 2021 00:26:31 +0100 Subject: [PATCH 08/11] Optimized get_all_spin_singles --- src/determinants/spindeterminants.irp.f | 46 +++++++++---------------- 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index 8ba405e0..765ee394 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -1193,16 +1193,10 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu xorvec(k) = xor( spindet(k), buffer(k,i) ) enddo - if (xorvec(1) /= 0_8) then - degree = popcnt(xorvec(1)) - else - degree = 0 - endif + degree = 0 - do k=2,$N_int - if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then + do k=1,$N_int degree = degree + popcnt(xorvec(k)) - endif enddo if ( degree == 4 ) then @@ -1247,23 +1241,19 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single xorvec(k) = xor( spindet(k), buffer(k,i) ) enddo - if (xorvec(1) /= 0_8) then - degree = popcnt(xorvec(1)) - else - degree = 0 - endif + degree = 0 - do k=2,$N_int - if ( (degree <= 2).and.(xorvec(k) /= 0_8) ) then + do k=1,$N_int degree = degree + popcnt(xorvec(k)) - endif enddo - if ( degree == 2 ) then - singles(n_singles) = idx(i) - n_singles = n_singles+1 + if ( degree /= 2 ) then + cycle endif + singles(n_singles) = idx(i) + n_singles = n_singles+1 + enddo n_singles = n_singles-1 @@ -1296,23 +1286,19 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double xorvec(k) = xor( spindet(k), buffer(k,i) ) enddo - if (xorvec(1) /= 0_8) then - degree = popcnt(xorvec(1)) - else - degree = 0 - endif + degree = 0 - do k=2,$N_int - if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then + do k=1,$N_int degree = degree + popcnt(xorvec(k)) - endif enddo - if ( degree == 4 ) then - doubles(n_doubles) = idx(i) - n_doubles = n_doubles+1 + if ( degree /= 4 ) then + cycle endif + doubles(n_doubles) = idx(i) + n_doubles = n_doubles+1 + enddo n_doubles = n_doubles-1 From d79724b7e77357fa075d8da83678f1b511733d7d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Mar 2021 01:04:23 +0100 Subject: [PATCH 09/11] Optimized phasemask --- src/cipsi/selection.irp.f | 77 ++++++++++++++++++++++++++++----------- 1 file changed, 56 insertions(+), 21 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 5afb514e..9b92dc9d 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -8,27 +8,56 @@ subroutine get_mask_phase(det1, pm, Nint) integer(bit_kind), intent(out) :: pm(Nint,2) integer(bit_kind) :: tmp1, tmp2 integer :: i - pm(1:Nint,1:2) = det1(1:Nint,1:2) tmp1 = 0_8 tmp2 = 0_8 - do i=1,Nint - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 1)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 1)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) - pm(i,1) = ieor(pm(i,1), tmp1) - pm(i,2) = ieor(pm(i,2), tmp2) - if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) - if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) - end do + select case (Nint) + +BEGIN_TEMPLATE + case ($Nint) + do i=1,$Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do +SUBST [ Nint ] +1;; +2;; +3;; +4;; +END_TEMPLATE + case default + do i=1,Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do + end select end subroutine @@ -445,11 +474,17 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d endif do i=1,fullinteresting(0) - fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i)) + do k=1,N_int + fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i)) + fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i)) + enddo enddo do i=1,interesting(0) - minilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,interesting(i)) + do k=1,N_int + minilist(k,1,i) = psi_det_sorted(k,1,interesting(i)) + minilist(k,2,i) = psi_det_sorted(k,2,interesting(i)) + enddo enddo do s2=s1,2 From 882d8ee3509a0fe522630b882ad8fa9af9ff7168 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Mar 2021 01:04:23 +0100 Subject: [PATCH 10/11] Improved print_energy --- src/cipsi/selection.irp.f | 77 ++++++++++++++++++++++++++---------- src/tools/print_energy.irp.f | 21 +--------- 2 files changed, 57 insertions(+), 41 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 5afb514e..9b92dc9d 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -8,27 +8,56 @@ subroutine get_mask_phase(det1, pm, Nint) integer(bit_kind), intent(out) :: pm(Nint,2) integer(bit_kind) :: tmp1, tmp2 integer :: i - pm(1:Nint,1:2) = det1(1:Nint,1:2) tmp1 = 0_8 tmp2 = 0_8 - do i=1,Nint - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 1)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 1)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) - pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) - pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) - pm(i,1) = ieor(pm(i,1), tmp1) - pm(i,2) = ieor(pm(i,2), tmp2) - if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) - if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) - end do + select case (Nint) + +BEGIN_TEMPLATE + case ($Nint) + do i=1,$Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do +SUBST [ Nint ] +1;; +2;; +3;; +4;; +END_TEMPLATE + case default + do i=1,Nint + pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1)) + pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16)) + pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32)) + pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32)) + pm(i,1) = ieor(pm(i,1), tmp1) + pm(i,2) = ieor(pm(i,2), tmp2) + if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1) + if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2) + end do + end select end subroutine @@ -445,11 +474,17 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d endif do i=1,fullinteresting(0) - fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i)) + do k=1,N_int + fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i)) + fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i)) + enddo enddo do i=1,interesting(0) - minilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,interesting(i)) + do k=1,N_int + minilist(k,1,i) = psi_det_sorted(k,1,interesting(i)) + minilist(k,2,i) = psi_det_sorted(k,2,interesting(i)) + enddo enddo do s2=s1,2 diff --git a/src/tools/print_energy.irp.f b/src/tools/print_energy.irp.f index f78dffc8..4fe1572c 100644 --- a/src/tools/print_energy.irp.f +++ b/src/tools/print_energy.irp.f @@ -14,24 +14,5 @@ end subroutine run implicit none - integer :: i,j - double precision :: i_H_psi_array(N_states) - double precision :: E(N_states) - double precision :: norm(N_states) - - E(1:N_states) = nuclear_repulsion - norm(1:N_states) = 0.d0 - do i=1,N_det - call i_H_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, & - size(psi_coef,1), N_states, i_H_psi_array) - do j=1,N_states - norm(j) += psi_coef(i,j)*psi_coef(i,j) - E(j) += i_H_psi_array(j) * psi_coef(i,j) - enddo - enddo - - print *, 'Energy:' - do i=1,N_states - print *, E(i)/norm(i) - enddo + print *, psi_energy + nuclear_repulsion end From 697fbddde64144e4fc5b0eed4465be7d35fb8339 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 31 Mar 2021 01:50:53 +0200 Subject: [PATCH 11/11] Optimized get_all_spin_doubles_1 --- src/davidson/davidson_parallel.irp.f | 6 ++++-- src/determinants/spindeterminants.irp.f | 28 +++++++++++-------------- 2 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index 32f89979..10b41f4c 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -438,7 +438,7 @@ 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 + if (imin<=N_det/2) then istep = 2 else istep = 1 @@ -505,7 +505,9 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) print *, irp_here, ': Failed in zmq_set_running' endif - call omp_set_max_active_levels(4) + + call omp_set_max_active_levels(5) + !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread) ithread = omp_get_thread_num() if (ithread == 0 ) then diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index 765ee394..dea4a566 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -1075,19 +1075,17 @@ subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer, integer :: i include 'utils/constants.include.F' integer :: degree - + integer :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /) + integer :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /) n_singles = 1 n_doubles = 1 do i=1,size_buffer degree = popcnt( xor( spindet, buffer(i) ) ) - if ( degree == 4 ) then - doubles(n_doubles) = idx(i) - n_doubles = n_doubles+1 - else if ( degree == 2 ) then - singles(n_singles) = idx(i) - n_singles = n_singles+1 - endif + doubles(n_doubles) = idx(i) + singles(n_singles) = idx(i) + n_doubles = n_doubles+add_double(degree) + n_singles = n_singles+add_single(degree) enddo n_singles = n_singles-1 n_doubles = n_doubles-1 @@ -1113,15 +1111,14 @@ subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_ integer :: i integer(bit_kind) :: v integer :: degree + integer :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /) include 'utils/constants.include.F' n_singles = 1 do i=1,size_buffer degree = popcnt(xor( spindet, buffer(i) )) - if (degree == 2) then - singles(n_singles) = idx(i) - n_singles = n_singles+1 - endif + singles(n_singles) = idx(i) + n_singles = n_singles+add_single(degree) enddo n_singles = n_singles-1 @@ -1145,14 +1142,13 @@ subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_ integer :: i include 'utils/constants.include.F' integer :: degree + integer :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /) n_doubles = 1 do i=1,size_buffer degree = popcnt(xor( spindet, buffer(i) )) - if ( degree == 4 ) then - doubles(n_doubles) = idx(i) - n_doubles = n_doubles+1 - endif + doubles(n_doubles) = idx(i) + n_doubles = n_doubles+add_double(degree) enddo n_doubles = n_doubles-1