From 9fd26ca1c8a6b88a33847673952d26f6cb3b68c3 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 28 Sep 2021 00:30:10 +0200 Subject: [PATCH 01/20] added dav_dressed_ext_rout.irp.f --- INSTALL.rst | 2 +- src/ao_basis/aos_in_r.irp.f | 26 +- .../dav_dressed_ext_rout.irp.f | 485 ++++++++++++++++++ src/dav_general_mat/dav_ext_rout.irp.f | 24 +- src/dav_general_mat/test_dav.irp.f | 39 +- src/mo_basis/mos_in_r.irp.f | 3 + 6 files changed, 543 insertions(+), 36 deletions(-) create mode 100644 src/dav_general_mat/dav_dressed_ext_rout.irp.f diff --git a/INSTALL.rst b/INSTALL.rst index 229bf40a..336d350c 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -209,7 +209,7 @@ ZeroMQ and its Fortran binding .. code:: bash - cp f77_zmq_free.h ${QP_ROOT}/src/ZMQ/f77_zmq.h + cp f77_zmq_free.h ${QP_ROOT}/src/zmq/f77_zmq.h Zlib diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f index 7fcb980a..902827eb 100644 --- a/src/ao_basis/aos_in_r.irp.f +++ b/src/ao_basis/aos_in_r.irp.f @@ -12,21 +12,21 @@ double precision function ao_value(i,r) integer :: power_ao(3) double precision :: accu,dx,dy,dz,r2 num_ao = ao_nucl(i) - power_ao(1:3)= ao_power(i,1:3) - center_ao(1:3) = nucl_coord(num_ao,1:3) - dx = (r(1) - center_ao(1)) - dy = (r(2) - center_ao(2)) - dz = (r(3) - center_ao(3)) - r2 = dx*dx + dy*dy + dz*dz - dx = dx**power_ao(1) - dy = dy**power_ao(2) - dz = dz**power_ao(3) +! power_ao(1:3)= ao_power(i,1:3) +! center_ao(1:3) = nucl_coord(num_ao,1:3) +! dx = (r(1) - center_ao(1)) +! dy = (r(2) - center_ao(2)) +! dz = (r(3) - center_ao(3)) +! r2 = dx*dx + dy*dy + dz*dz +! dx = dx**power_ao(1) +! dy = dy**power_ao(2) +! dz = dz**power_ao(3) accu = 0.d0 - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) - enddo +! do m=1,ao_prim_num(i) +! beta = ao_expo_ordered_transp(m,i) +! accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) +! enddo ao_value = accu * dx * dy * dz end diff --git a/src/dav_general_mat/dav_dressed_ext_rout.irp.f b/src/dav_general_mat/dav_dressed_ext_rout.irp.f new file mode 100644 index 00000000..27873bfc --- /dev/null +++ b/src/dav_general_mat/dav_dressed_ext_rout.irp.f @@ -0,0 +1,485 @@ +subroutine davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_diag,dressing_state,dressing_vec,idress,converged,hcalc) + use mmap_module + implicit none + BEGIN_DOC + ! Davidson diagonalization. + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! sze : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: sze, N_st, N_st_diag,idress + double precision, intent(inout) :: u_in(sze,N_st_diag) + double precision, intent(inout) :: H_jj(sze) + double precision, intent(out) :: energies(N_st_diag) + double precision, intent(in) :: dressing_vec(sze,N_st) + integer, intent(in) :: dressing_state + logical, intent(out) :: converged + external hcalc + + double precision :: f + + integer :: iter + integer :: i,j,k,l,m + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: s_tmp(:,:) + double precision :: diag_h_mat_elem + double precision, allocatable :: residual_norm(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st) + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + logical :: state_ok(N_st_diag*davidson_sze_max) + integer :: nproc_target + integer :: order(N_st_diag) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag*davidson_sze_max) + !!!! TO CHANGE !!!! + integer :: idx_dress(1) + idx_dress = idress + + + if (dressing_state > 0) then + do k=1,N_st + do i=1,sze + H_jj(i) += u_in(i,k) * dressing_vec(i,k) + enddo + enddo + endif + + l = idx_dress(1) + f = 1.0d0/u_in(l,1) + + include 'constants.include.F' + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag*itermax, N_st_diag*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + + provide threshold_davidson !nthreads_davidson + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab +! maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + maxab = sze + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag*itermax) &! U + + 1.0d0*dble(sze*m)*(N_st_diag*itermax) &! W + + 3.0d0*(N_st_diag*itermax)**2 &! h,y,s_tmp + + 1.d0*(N_st_diag*itermax) &! lambda + + 1.d0*(N_st_diag) &! residual_norm + ! In H_u_0_nstates_zmq + + 2.d0*(N_st_diag*N_det) &! u_t, v_t, on collector + + 2.d0*(N_st_diag*N_det) &! u_t, v_t, on slave + + 0.5d0*maxab &! idx0 in H_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag,'Number of states in diagonalization') + call write_int(6,sze,'Number of determinants') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag*itermax), & + + ! Small + h(N_st_diag*itermax,N_st_diag*itermax), & + y(N_st_diag*itermax,N_st_diag*itermax), & + s_tmp(N_st_diag*itermax,N_st_diag*itermax), & + residual_norm(N_st_diag), & + lambda(N_st_diag*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + s_tmp = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag >= N_st) + ASSERT (sze > 0) + + ! Davidson iterations + ! =================== + + converged = .False. + + do k=N_st+1,N_st_diag + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) + enddo + u_in(k,k) = u_in(k,k) + 10.d0 + enddo + ! Normalize all states + do k=1,N_st_diag + call normalize(u_in(1,k),sze) + enddo + ! Copy from the guess input "u_in" to the working vectors "U" + + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 8) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag*(iter-1) + shift2 = N_st_diag*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + ! Gram-Schmidt to orthogonalize all new guess with the previous vectors + call ortho_qr(U,size(U,1),sze,shift2) + call ortho_qr(U,size(U,1),sze,shift2) + ! it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag) + ! where sze is the size of the vector, N_st_diag is the number of states + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze) + else + ! Already computed in update below + continue + endif + + if (dressing_state > 0) then + + if (N_st == 1) then + + + do istate=1,N_st_diag + do i=1,sze + W(i,shift+istate) += dressing_vec(i,1) *f * U(l,shift+istate) + W(l,shift+istate) += dressing_vec(i,1) *f * U(i,shift+istate) + enddo + + enddo + + else + print*,'Not implemented yet for multi state ...' + stop +! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & +! psi_coef, size(psi_coef,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, & +! dressing_vec, size(dressing_vec,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! +! +! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & +! dressing_vec, size(dressing_vec,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, & +! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) + + endif + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), U, size(U,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + ! Diagonalize h + ! --------------- + + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif + + ! Compute Energy for each eigenvector + ! ----------------------------------- + + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, h, size(h,1), y, size(y,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + 0.d0, h, size(h,1)) + + do k=1,shift2 + lambda(k) = h(k,k) + enddo + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + nuclear_repulsion + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + do k=1,N_st + if (residual_norm(k) > 1.d8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + ! Re-contract U and update W + ! -------------------------------- + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + enddo + + + call nullify_small_elements(sze,N_st_diag,U,size(U,1),threshold_davidson_pt2) + do k=1,N_st_diag + do i=1,sze + u_in(i,k) = U(i,k) + enddo + enddo + + do k=1,N_st_diag + energies(k) = lambda(k) + enddo + write_buffer = '======' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + call write_time(6) + + deallocate(W) + + deallocate ( & + residual_norm, & + U, overlap, & + h, y, s_tmp, & + lambda & + ) + FREE nthreads_davidson +end + + + + + + + diff --git a/src/dav_general_mat/dav_ext_rout.irp.f b/src/dav_general_mat/dav_ext_rout.irp.f index a2aad413..ff2167e4 100644 --- a/src/dav_general_mat/dav_ext_rout.irp.f +++ b/src/dav_general_mat/dav_ext_rout.irp.f @@ -1,5 +1,5 @@ -subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,hcalc) +subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) use mmap_module implicit none BEGIN_DOC @@ -9,7 +9,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia ! ! u_in : guess coefficients on the various states. Overwritten on exit ! - ! dim_in : leftmost dimension of u_in + ! sze : leftmost dimension of u_in ! ! sze : Number of determinants ! @@ -21,9 +21,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia ! ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) END_DOC - integer, intent(in) :: dim_in, sze, N_st, N_st_diag_in + integer, intent(in) :: sze, N_st, N_st_diag_in double precision, intent(in) :: H_jj(sze) - double precision, intent(inout) :: u_in(dim_in,N_st_diag_in) + double precision, intent(inout) :: u_in(sze,N_st_diag_in) double precision, intent(out) :: energies(N_st) external hcalc @@ -157,19 +157,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia write(6,'(A)') write_buffer(1:6+41*N_st) -! if (disk_based) then -! ! Create memory-mapped files for W and S -! type(c_ptr) :: ptr_w, ptr_s -! integer :: fd_s, fd_w -! call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& -! 8, fd_w, .False., ptr_w) -! call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),& -! 4, fd_s, .False., ptr_s) -! call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/)) -! call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/)) -! else - allocate(W(sze,N_st_diag*itermax)) -! endif + allocate(W(sze,N_st_diag*itermax)) allocate( & ! Large @@ -398,7 +386,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia write(6,'(A)') '' call write_time(6) - deallocate(W) + deallocate(W) deallocate ( & residual_norm, & diff --git a/src/dav_general_mat/test_dav.irp.f b/src/dav_general_mat/test_dav.irp.f index 0a7faaf0..4a65ba7c 100644 --- a/src/dav_general_mat/test_dav.irp.f +++ b/src/dav_general_mat/test_dav.irp.f @@ -8,12 +8,13 @@ program test_dav touch read_wf PROVIDE threshold_davidson nthreads_davidson call routine + call test_dav_dress end subroutine routine implicit none double precision, allocatable :: u_in(:,:), H_jj(:), energies(:),h_mat(:,:) - integer :: dim_in,sze,N_st,N_st_diag_in,dressing_state + integer :: dim_in,sze,N_st,N_st_diag_in logical :: converged integer :: i,j external hcalc_template @@ -21,9 +22,8 @@ subroutine routine N_st_diag_in = N_states_diag sze = N_det dim_in = sze - dressing_state = 0 !!!! MARK THAT u_in mut dimensioned with "N_st_diag_in" as a second dimension - allocate(u_in(dim_in,N_st_diag_in),H_jj(sze),h_mat(sze,sze),energies(N_st)) + allocate(u_in(dim_in,N_st_diag_in),H_jj(sze),h_mat(sze,sze),energies(N_st_diag_in)) u_in = 0.d0 do i = 1, N_st u_in(1,i) = 1.d0 @@ -42,7 +42,38 @@ subroutine routine print*,'energies = ',energies !!! hcalc_template is the routine that computes v = H u !!! and you can use the routine "davidson_general_ext_rout" - call davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,hcalc_template) + call davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc_template) print*,'energies = ',energies end + +subroutine test_dav_dress + implicit none + double precision, allocatable :: u_in(:,:), H_jj(:), energies(:) + integer :: sze,N_st,N_st_diag_in,dressing_state + logical :: converged + integer :: i,j + external hcalc_template + double precision, allocatable :: dressing_vec(:) + integer :: idress + N_st = N_states + N_st_diag_in = N_states_diag + sze = N_det + dressing_state = 0 + idress = 1 + !!!! MARK THAT u_in mut dimensioned with "N_st_diag_in" as a second dimension + allocate(u_in(sze,N_st_diag_in),H_jj(sze),energies(N_st_diag_in)) + allocate(dressing_vec(sze)) + dressing_vec = 0.d0 + u_in = 0.d0 + do i = 1, N_st + u_in(1,i) = 1.d0 + enddo + do i = 1, sze + H_jj(i) = H_matrix_all_dets(i,i) + nuclear_repulsion + enddo + print*,'dressing davidson ' + call davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_diag_in,dressing_state,dressing_vec,idress,converged,hcalc_template) + print*,'energies(1) = ',energies(1) + +end diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f index ee2795d0..e5d3b243 100644 --- a/src/mo_basis/mos_in_r.irp.f +++ b/src/mo_basis/mos_in_r.irp.f @@ -1,6 +1,9 @@ subroutine give_all_mos_at_r(r,mos_array) implicit none + BEGIN_DOC +! mos_array(i) = ith MO function evaluated at "r" + END_DOC double precision, intent(in) :: r(3) double precision, intent(out) :: mos_array(mo_num) double precision :: aos_array(ao_num) From 1d53e6fda2f7667ec83219df5febb6b5249f96ae Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 28 Sep 2021 11:53:34 +0200 Subject: [PATCH 02/20] added the possibility to introduce a threshold for saving the wave function in cis.irp.f --- src/cis/cis.irp.f | 2 +- src/determinants/EZFIO.cfg | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f index acec29c2..c4100105 100644 --- a/src/cis/cis.irp.f +++ b/src/cis/cis.irp.f @@ -79,6 +79,6 @@ subroutine run call ezfio_set_cis_energy(CI_energy) psi_coef = ci_eigenvectors SOFT_TOUCH psi_coef - call save_wavefunction_truncated(1.d-12) + call save_wavefunction_truncated(threshold_save_wf) end diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index 662c6fbb..1e85693b 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -136,3 +136,9 @@ doc: If |true|, discard any Slater determinants with an interaction smaller than interface: ezfio,provider,ocaml default: False + +[threshold_save_wf] +type: Threshold +doc: Threshold on the coefficients of the wave function when saving it into the ezfio +interface: ezfio,provider,ocaml +default: 1.e-14 From 0c7c8513b1c5ccb1a6bba46dcc75a5b4c6b59580 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 7 Oct 2021 17:37:24 +0200 Subject: [PATCH 03/20] added the possibility to have a mu(r) in the functionals --- src/dft_keywords/EZFIO.cfg | 5 +++++ src/dft_one_e/NEED | 1 + src/dft_one_e/mu_erf_dft.irp.f | 12 ++++++++++++ src/functionals/sr_lda.irp.f | 18 ++++++++++++------ src/functionals/sr_pbe.irp.f | 14 ++++++++++---- 5 files changed, 40 insertions(+), 10 deletions(-) diff --git a/src/dft_keywords/EZFIO.cfg b/src/dft_keywords/EZFIO.cfg index b452c863..a13cbac9 100644 --- a/src/dft_keywords/EZFIO.cfg +++ b/src/dft_keywords/EZFIO.cfg @@ -16,3 +16,8 @@ doc: Percentage of HF exchange in the DFT model interface: ezfio,provider,ocaml default: 0. +[mu_dft_type] +type: character*(32) +doc: type of mu(r) for rsdft [ cst ] +interface: ezfio, provider, ocaml +default: cst diff --git a/src/dft_one_e/NEED b/src/dft_one_e/NEED index 3a942f28..615ee97e 100644 --- a/src/dft_one_e/NEED +++ b/src/dft_one_e/NEED @@ -6,3 +6,4 @@ ao_one_e_ints ao_two_e_ints mo_two_e_erf_ints ao_two_e_erf_ints +mu_of_r diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index 53effcb6..8161324b 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -8,3 +8,15 @@ BEGIN_PROVIDER [double precision, mu_erf_dft] mu_erf_dft = mu_erf END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_of_r_dft, (n_points_final_grid)] + implicit none + integer :: i + do i = 1, n_points_final_grid + if(mu_dft_type == "cst")then + mu_of_r_dft(i) = mu_erf_dft + else + mu_of_r_dft(i) = mu_of_r_hf(i) + endif + enddo +END_PROVIDER diff --git a/src/functionals/sr_lda.irp.f b/src/functionals/sr_lda.irp.f index 965a744c..bd062a02 100644 --- a/src/functionals/sr_lda.irp.f +++ b/src/functionals/sr_lda.irp.f @@ -21,7 +21,9 @@ weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - call ex_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,vx_a,vx_b) + double precision :: mu_local + mu_local = mu_of_r_dft(i) + call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,vx_a,vx_b) energy_x_sr_lda(istate) += weight * e_x enddo enddo @@ -48,7 +50,9 @@ weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - call ec_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_c,vc_a,vc_b) + double precision :: mu_local + mu_local = mu_of_r_dft(i) + call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,vc_a,vc_b) energy_c_sr_lda(istate) += weight * e_c enddo enddo @@ -122,8 +126,10 @@ END_PROVIDER weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - call ec_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b) - call ex_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b) + double precision :: mu_local + mu_local = mu_of_r_dft(i) + call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b) + call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b) do j =1, ao_num aos_sr_vc_alpha_lda_w(j,i,istate) = sr_vc_a * aos_in_r_array(j,i)*weight aos_sr_vc_beta_lda_w(j,i,istate) = sr_vc_b * aos_in_r_array(j,i)*weight @@ -147,8 +153,6 @@ END_PROVIDER double precision :: mu,weight double precision :: e_c,sr_vc_a,sr_vc_b,e_x,sr_vx_a,sr_vx_b double precision, allocatable :: rhoa(:),rhob(:) - double precision :: mu_local - mu_local = mu_erf_dft allocate(rhoa(N_states), rhob(N_states)) do istate = 1, N_states do i = 1, n_points_final_grid @@ -158,6 +162,8 @@ END_PROVIDER weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) + double precision :: mu_local + mu_local = mu_of_r_dft(i) call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b) call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b) do j =1, ao_num diff --git a/src/functionals/sr_pbe.irp.f b/src/functionals/sr_pbe.irp.f index 93c51067..7053cfb6 100644 --- a/src/functionals/sr_pbe.irp.f +++ b/src/functionals/sr_pbe.irp.f @@ -35,9 +35,11 @@ grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m) grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo - + + double precision :: mu_local + mu_local = mu_of_r_dft(i) ! inputs - call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) energy_x_sr_pbe(istate) += ex * weight @@ -135,8 +137,10 @@ END_PROVIDER grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo + double precision :: mu_local + mu_local = mu_of_r_dft(i) ! inputs - call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) vx_rho_a *= weight @@ -292,8 +296,10 @@ END_PROVIDER grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo + double precision :: mu_local + mu_local = mu_of_r_dft(i) ! inputs - call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) vx_rho_a *= weight From 2125cd69aba5f560cdc4e697cc2b871f8f7fe0ad Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 7 Oct 2021 23:05:43 +0200 Subject: [PATCH 04/20] added mu(rsc) --- src/dft_one_e/mu_erf_dft.irp.f | 57 +++++++++++++++++++++++++++++---- src/dft_utils_func/mu_rsc.irp.f | 13 ++++++++ 2 files changed, 63 insertions(+), 7 deletions(-) create mode 100644 src/dft_utils_func/mu_rsc.irp.f diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index 8161324b..a1d7a5c1 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -12,11 +12,54 @@ END_PROVIDER BEGIN_PROVIDER [double precision, mu_of_r_dft, (n_points_final_grid)] implicit none integer :: i - do i = 1, n_points_final_grid - if(mu_dft_type == "cst")then - mu_of_r_dft(i) = mu_erf_dft - else - mu_of_r_dft(i) = mu_of_r_hf(i) - endif - enddo + if(mu_dft_type == "Read")then + call ezfio_get_mu_of_r_mu_of_r_disk(mu_of_r_dft) + else + do i = 1, n_points_final_grid + if(mu_dft_type == "cst")then + mu_of_r_dft(i) = mu_erf_dft + else if(mu_dft_type == "hf")then + mu_of_r_dft(i) = mu_of_r_hf(i) + else if(mu_dft_type == "rsc")then + mu_of_r_dft(i) = mu_rsc_of_r(i) + else + print*,'mu_dft_type is not of good type = ',mu_dft_type + print*,'it must be of type Read, cst, hf, rsc' + print*,'Stopping ...' + stop + endif + enddo + endif +END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_rsc_of_r, (n_points_final_grid)] + implicit none + integer :: i + double precision :: mu_rs_c,rho,r(3), dm_a, dm_b + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) + rho = dm_a + dm_b + mu_rsc_of_r(i) = mu_rs_c(rho) + enddo + +END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_of_r_dft_average] + implicit none + integer :: i + double precision :: mu_rs_c,rho,r(3), dm_a, dm_b + mu_of_r_dft_average = 0.d0 + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) + rho = dm_a + dm_b + mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i) + enddo + mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num) + print*,'mu_of_r_dft_average = ',mu_of_r_dft_average END_PROVIDER diff --git a/src/dft_utils_func/mu_rsc.irp.f b/src/dft_utils_func/mu_rsc.irp.f new file mode 100644 index 00000000..386a1c2d --- /dev/null +++ b/src/dft_utils_func/mu_rsc.irp.f @@ -0,0 +1,13 @@ +double precision function mu_rs_c(rho) + implicit none + double precision, intent(in) :: rho + include 'constants.include.F' + double precision :: cst_rs,alpha_rs,rs + cst_rs = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0) + alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi + + rs = cst_rs * rho**(-1.d0/3.d0) + mu_rs_c = alpha_rs/dsqrt(rs) + +end + From 8a7d6444abe7d837b33863f0e59ce85b6a95e23b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 18 Oct 2021 08:56:25 +0200 Subject: [PATCH 05/20] removed some stuffs to clean the CIPSI module --- src/cipsi/EZFIO.cfg | 6 - src/cipsi/NEED | 1 - src/cipsi/pert_rdm_providers.irp.f | 183 -------------- src/cipsi/pt2_stoch_routines.irp.f | 2 +- src/cipsi/selection.irp.f | 30 +-- src/cipsi/update_2rdm.irp.f | 223 ------------------ src/cipsi/zmq_selection.irp.f | 2 +- src/davidson_dressed/diagonalize_ci.irp.f | 2 +- src/dft_one_e/mu_erf_dft.irp.f | 1 + src/iterations/print_summary.irp.f | 2 +- src/tools/NEED | 1 + .../print_e_components.irp.f | 0 12 files changed, 21 insertions(+), 432 deletions(-) delete mode 100644 src/cipsi/pert_rdm_providers.irp.f delete mode 100644 src/cipsi/update_2rdm.irp.f rename src/{davidson => two_body_rdm}/print_e_components.irp.f (100%) diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg index 19b45ac1..7fcf19eb 100644 --- a/src/cipsi/EZFIO.cfg +++ b/src/cipsi/EZFIO.cfg @@ -1,9 +1,3 @@ -[pert_2rdm] -type: logical -doc: If true, computes the one- and two-body rdms with perturbation theory -interface: ezfio,provider,ocaml -default: False - [save_wf_after_selection] type: logical doc: If true, saves the wave function after the selection, before the diagonalization diff --git a/src/cipsi/NEED b/src/cipsi/NEED index bfbc559a..85d01f79 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -2,5 +2,4 @@ perturbation zmq mpi iterations -two_body_rdm csf diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f deleted file mode 100644 index eca8decc..00000000 --- a/src/cipsi/pert_rdm_providers.irp.f +++ /dev/null @@ -1,183 +0,0 @@ - -use bitmasks -use omp_lib - -BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock] - use f77_zmq - implicit none - call omp_init_lock(pert_2rdm_lock) -END_PROVIDER - -BEGIN_PROVIDER [integer, n_orb_pert_rdm] - implicit none - n_orb_pert_rdm = n_act_orb -END_PROVIDER - -BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)] - implicit none - list_orb_reverse_pert_rdm = list_act_reverse - -END_PROVIDER - -BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)] - implicit none - list_orb_pert_rdm = list_act - -END_PROVIDER - -BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)] - implicit none - pert_2rdm_provider = 0.d0 - -END_PROVIDER - -subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection) - use bitmasks - use selection_types - implicit none - - integer, intent(in) :: n_det_connection - double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) - integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - integer, intent(in) :: i_generator, sp, h1, h2 - double precision, intent(in) :: mat(N_states, mo_num, mo_num) - logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) - double precision, intent(in) :: fock_diag_tmp(mo_num) - double precision, intent(in) :: E0(N_states) - type(pt2_type), intent(inout) :: pt2_data - type(selection_buffer), intent(inout) :: buf - logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate, jstate - integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states) - double precision, external :: diag_H_mat_elem_fock - double precision :: E_shift - - logical, external :: detEq - double precision, allocatable :: values(:) - integer, allocatable :: keys(:,:) - integer :: nkeys - integer :: sze_buff - sze_buff = 5 * mo_num ** 2 - allocate(keys(4,sze_buff),values(sze_buff)) - nkeys = 0 - if(sp == 3) then - s1 = 1 - s2 = 2 - else - s1 = sp - s2 = sp - end if - call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) - E_shift = 0.d0 - - if (h0_type == 'CFG') then - j = det_to_configuration(i_generator) - E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j) - endif - - do p1=1,mo_num - if(bannedOrb(p1, s1)) cycle - ib = 1 - if(sp /= 3) ib = p1+1 - - do p2=ib,mo_num - -! ----- -! /!\ Generating only single excited determinants doesn't work because a -! determinant generated by a single excitation may be doubly excited wrt -! to a determinant of the future. In that case, the determinant will be -! detected as already generated when generating in the future with a -! double excitation. -! -! if (.not.do_singles) then -! if ((h1 == p1) .or. (h2 == p2)) then -! cycle -! endif -! endif -! -! if (.not.do_doubles) then -! if ((h1 /= p1).and.(h2 /= p2)) then -! cycle -! endif -! endif -! ----- - - if(bannedOrb(p2, s2)) cycle - if(banned(p1,p2)) cycle - - - if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle - call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) - - if (do_only_cas) then - integer, external :: number_of_holes, number_of_particles - if (number_of_particles(det)>0) then - cycle - endif - if (number_of_holes(det)>0) then - cycle - endif - endif - - if (do_ddci) then - logical, external :: is_a_two_holes_two_particles - if (is_a_two_holes_two_particles(det)) then - cycle - endif - endif - - if (do_only_1h1p) then - logical, external :: is_a_1h1p - if (.not.is_a_1h1p(det)) cycle - endif - - - Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) - - sum_e_pert = 0d0 - integer :: degree - call get_excitation_degree(det,HF_bitmask,degree,N_int) - if(degree == 2)cycle - do istate=1,N_states - delta_E = E0(istate) - Hii + E_shift - alpha_h_psi = mat(istate, p1, p2) - val = alpha_h_psi + alpha_h_psi - tmp = dsqrt(delta_E * delta_E + val * val) - if (delta_E < 0.d0) then - tmp = -tmp - endif - e_pert = 0.5d0 * (tmp - delta_E) - coef(istate) = e_pert / alpha_h_psi - print*,e_pert,coef,alpha_h_psi - pt2_data % pt2(istate) += e_pert - pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - enddo - - do istate=1,N_states - alpha_h_psi = mat(istate, p1, p2) - e_pert = coef(istate) * alpha_h_psi - do jstate=1,N_states - pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate) - enddo - - if (weight_selection /= 5) then - ! Energy selection - sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) - - else - ! Variance selection - sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) - endif - end do - call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff) - - if(sum_e_pert <= buf%mini) then - call add_to_selection_buffer(buf, det, sum_e_pert) - end if - end do - end do - call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) -end - - diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index b366a268..3594aaf2 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -133,7 +133,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_det_hii selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max if (h0_type == 'CFG') then PROVIDE psi_configuration_hii det_to_configuration diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index eda9642c..f1ec6ff6 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -464,14 +464,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d allocate (fullminilist (N_int, 2, fullinteresting(0)), & minilist (N_int, 2, interesting(0)) ) - if(pert_2rdm)then - allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) - do i=1,fullinteresting(0) - do j = 1, N_states - coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) - enddo - enddo - endif +! if(pert_2rdm)then +! allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) +! do i=1,fullinteresting(0) +! do j = 1, N_states +! coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) +! enddo +! enddo +! endif do i=1,fullinteresting(0) do k=1,N_int @@ -531,19 +531,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) - if(.not.pert_2rdm)then +! if(.not.pert_2rdm)then call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) - else - call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) - endif +! else +! call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) +! endif end if enddo if(s1 /= s2) monoBdo = .false. enddo deallocate(fullminilist,minilist) - if(pert_2rdm)then - deallocate(coef_fullminilist_rev) - endif +! if(pert_2rdm)then +! deallocate(coef_fullminilist_rev) +! endif enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f deleted file mode 100644 index 260c48fd..00000000 --- a/src/cipsi/update_2rdm.irp.f +++ /dev/null @@ -1,223 +0,0 @@ -use bitmasks - -subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: n_det_connection,sze_buff - double precision, intent(in) :: coef(N_states) - integer(bit_kind), intent(in) :: det(N_int,2) - integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) - integer, intent(inout) :: keys(4,sze_buff),nkeys - double precision, intent(inout) :: values(sze_buff) - integer :: i,j - integer :: exc(0:2,2,2) - integer :: degree - double precision :: phase, contrib - do i = 1, n_det_connection - call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) - if(degree.gt.2)cycle - contrib = 0.d0 - do j = 1, N_states - contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * phase * coef(j) - enddo - ! case of single excitations - if(degree == 1)then - if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then - call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) - nkeys = 0 - endif - call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff) - else - !! case of double excitations - ! if (nkeys + 4 .ge. sze_buff)then - ! call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) - ! nkeys = 0 - ! endif - ! call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) - endif - enddo -!call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) -!nkeys = 0 - -end - -subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: sze_buff - integer(bit_kind), intent(in) :: det1(N_int,2) - integer(bit_kind), intent(in) :: det2(N_int,2) - integer,intent(in) :: exc(0:2,2,2) - double precision,intent(in) :: phase, contrib - integer, intent(inout) :: nkeys, keys(4,sze_buff) - double precision, intent(inout):: values(sze_buff) - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2),ispin,other_spin - integer :: h1,h2,p1,p2,i - call bitstring_to_list_ab(det1, occ, n_occ_ab, N_int) - - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - p1 = exc(1,2,1) - ispin = 1 - other_spin = 2 - else - ! Mono beta - h1 = exc(1,1,2) - p1 = exc(1,2,2) - ispin = 2 - other_spin = 1 - endif - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - !update the alpha/beta part - do i = 1, n_occ_ab(other_spin) - h2 = occ(i,other_spin) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - !update the same spin part -!do i = 1, n_occ_ab(ispin) -! h2 = occ(i,ispin) -! if(list_orb_reverse_pert_rdm(h2).lt.0)return -! h2 = list_orb_reverse_pert_rdm(h2) - -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 - -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 -! -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 - -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 -!enddo - -end - -subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: sze_buff - integer,intent(in) :: exc(0:2,2,2) - double precision,intent(in) :: phase, contrib - integer, intent(inout) :: nkeys, keys(4,sze_buff) - double precision, intent(inout):: values(sze_buff) - integer :: h1,h2,p1,p2 - - if (exc(0,1,1) == 1) then - ! Double alpha/beta - h1 = exc(1,1,1) - h2 = exc(1,1,2) - p1 = exc(1,2,1) - p2 = exc(1,2,2) - ! check if the orbitals involved are within the orbital range - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - if(list_orb_reverse_pert_rdm(p2).lt.0)return - p2 = list_orb_reverse_pert_rdm(p2) - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - - else - if (exc(0,1,1) == 2) then - ! Double alpha/alpha - h1 = exc(1,1,1) - h2 = exc(2,1,1) - p1 = exc(1,2,1) - p2 = exc(2,2,1) - else if (exc(0,1,2) == 2) then - ! Double beta - h1 = exc(1,1,2) - h2 = exc(2,1,2) - p1 = exc(1,2,2) - p2 = exc(2,2,2) - endif - ! check if the orbitals involved are within the orbital range - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - if(list_orb_reverse_pert_rdm(p2).lt.0)return - p2 = list_orb_reverse_pert_rdm(p2) - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - - nkeys += 1 - values(nkeys) = - 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = - 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - endif - -end - - diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 58630709..1bfe87c0 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data) PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') diff --git a/src/davidson_dressed/diagonalize_ci.irp.f b/src/davidson_dressed/diagonalize_ci.irp.f index 7619532a..709ee0e6 100644 --- a/src/davidson_dressed/diagonalize_ci.irp.f +++ b/src/davidson_dressed/diagonalize_ci.irp.f @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] enddo do j=1,min(N_det,N_states) write(st,'(I4)') j - call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st)) + call write_double(6,CI_energy_dressed(j),'Energy dressed of state '//trim(st)) call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st)) enddo diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index a1d7a5c1..21f6dc21 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -58,6 +58,7 @@ BEGIN_PROVIDER [double precision, mu_of_r_dft_average] r(3) = final_grid_points(3,i) call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) rho = dm_a + dm_b + if(mu_of_r_dft(i).gt.1.d+3)cycle mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i) enddo mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num) diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 641ee209..8e6285e2 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -98,7 +98,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s enddo endif - call print_energy_components() +! call print_energy_components() end subroutine diff --git a/src/tools/NEED b/src/tools/NEED index c07c9109..0f4e17b0 100644 --- a/src/tools/NEED +++ b/src/tools/NEED @@ -2,3 +2,4 @@ fci mo_two_e_erf_ints aux_quantities hartree_fock +two_body_rdm diff --git a/src/davidson/print_e_components.irp.f b/src/two_body_rdm/print_e_components.irp.f similarity index 100% rename from src/davidson/print_e_components.irp.f rename to src/two_body_rdm/print_e_components.irp.f From 416cd45d7a469961db5ff03a8d27639fec1f1a20 Mon Sep 17 00:00:00 2001 From: Fileto Rodriguez Date: Fri, 22 Oct 2021 12:18:01 +0200 Subject: [PATCH 06/20] added mu of r grad_rho --- src/dft_one_e/mu_erf_dft.irp.f | 12 +++++++++ src/dft_utils_func/mu_of_r_dft.irp.f | 37 ++++++++++++++++++++++++++++ src/dft_utils_func/mu_rsc.irp.f | 13 ---------- 3 files changed, 49 insertions(+), 13 deletions(-) create mode 100644 src/dft_utils_func/mu_of_r_dft.irp.f delete mode 100644 src/dft_utils_func/mu_rsc.irp.f diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index a1d7a5c1..f9bff24d 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -44,9 +44,21 @@ BEGIN_PROVIDER [double precision, mu_rsc_of_r, (n_points_final_grid)] rho = dm_a + dm_b mu_rsc_of_r(i) = mu_rs_c(rho) enddo +END_PROVIDER +BEGIN_PROVIDER [double precision, mu_grad_rho, (n_points_final_grid)] + implicit none + integer :: i + double precision :: mu_grad_rho_func, r(3) + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + mu_grad_rho(i) = mu_grad_rho_func(r) + enddo END_PROVIDER + BEGIN_PROVIDER [double precision, mu_of_r_dft_average] implicit none integer :: i diff --git a/src/dft_utils_func/mu_of_r_dft.irp.f b/src/dft_utils_func/mu_of_r_dft.irp.f new file mode 100644 index 00000000..7cba0a60 --- /dev/null +++ b/src/dft_utils_func/mu_of_r_dft.irp.f @@ -0,0 +1,37 @@ +double precision function mu_rs_c(rho) + implicit none + double precision, intent(in) :: rho + include 'constants.include.F' + double precision :: cst_rs,alpha_rs,rs + cst_rs = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0) + alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi + + rs = cst_rs * rho**(-1.d0/3.d0) + mu_rs_c = alpha_rs/dsqrt(rs) + +end + +double precision function mu_grad_rho_func(r) + implicit none + double precision , intent(in) :: r(3) + integer :: m + double precision :: rho, dm_a, dm_b, grad_dm_a(3), grad_dm_b(3) + double precision :: eta, grad_rho(3), grad_sqr + eta = 0.135d0 + call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b) + rho = dm_a + dm_b + do m = 1,3 + grad_rho(m) = grad_dm_a(m) + grad_dm_b(m) + enddo + grad_sqr=0.d0 + do m = 1,3 + grad_sqr=grad_sqr+grad_rho(m)*grad_rho(m) + enddo + grad_sqr = dsqrt(grad_sqr) + if (rho<1.d-12) then + mu_grad_rho_func = 1.d-10 + else + mu_grad_rho_func = eta * grad_sqr / rho + endif + +end diff --git a/src/dft_utils_func/mu_rsc.irp.f b/src/dft_utils_func/mu_rsc.irp.f deleted file mode 100644 index 386a1c2d..00000000 --- a/src/dft_utils_func/mu_rsc.irp.f +++ /dev/null @@ -1,13 +0,0 @@ -double precision function mu_rs_c(rho) - implicit none - double precision, intent(in) :: rho - include 'constants.include.F' - double precision :: cst_rs,alpha_rs,rs - cst_rs = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0) - alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi - - rs = cst_rs * rho**(-1.d0/3.d0) - mu_rs_c = alpha_rs/dsqrt(rs) - -end - From 699f655b8943b5aa7486be845294d8dec8d86305 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 25 Oct 2021 10:35:22 +0200 Subject: [PATCH 07/20] added int grad --- src/dft_utils_in_r/ints_grad.irp.f | 39 ++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 src/dft_utils_in_r/ints_grad.irp.f diff --git a/src/dft_utils_in_r/ints_grad.irp.f b/src/dft_utils_in_r/ints_grad.irp.f new file mode 100644 index 00000000..239fe554 --- /dev/null +++ b/src/dft_utils_in_r/ints_grad.irp.f @@ -0,0 +1,39 @@ + BEGIN_PROVIDER [ double precision, mo_grad_ints, (mo_num, mo_num,3)] + implicit none + BEGIN_DOC +! mo_grad_ints(i,j,m) = + END_DOC + integer :: i,j,ipoint,m + double precision :: weight + mo_grad_ints = 0.d0 + do m = 1, 3 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do j = 1, mo_num + do i = 1, mo_num + mo_grad_ints(i,j,m) += mos_grad_in_r_array(j,ipoint,m) * mos_in_r_array(i,ipoint) * weight + enddo + enddo + enddo + enddo + + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, mo_grad_ints_transp, (3,mo_num, mo_num)] + implicit none + BEGIN_DOC +! mo_grad_ints(i,j,m) = + END_DOC + integer :: i,j,ipoint,m + double precision :: weight + do m = 1, 3 + do j = 1, mo_num + do i = 1, mo_num + mo_grad_ints_transp(m,i,j) = mo_grad_ints(i,j,m) + enddo + enddo + enddo + + +END_PROVIDER From 1f0c48023d20540a3101f5c46b9317479a980eb6 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 26 Oct 2021 19:11:04 +0200 Subject: [PATCH 08/20] dav_double_dress_ext_rout.irp.f --- .../dav_double_dress_ext_rout.irp.f | 520 ++++++++++++++++++ src/dav_general_mat/dav_ext_rout.irp.f | 2 +- src/utils/sort.irp.f | 209 ++----- 3 files changed, 557 insertions(+), 174 deletions(-) create mode 100644 src/dav_general_mat/dav_double_dress_ext_rout.irp.f diff --git a/src/dav_general_mat/dav_double_dress_ext_rout.irp.f b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f new file mode 100644 index 00000000..9bcfc5ab --- /dev/null +++ b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f @@ -0,0 +1,520 @@ +subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies,sze,N_st,N_st_diag,converged,hcalc) + use mmap_module + BEGIN_DOC + ! Generic Davidson diagonalization with TWO DRESSING VECTORS + ! + ! Dress_jj : DIAGONAL DRESSING of the Hamiltonian + ! + ! Dressing_vec : COLUMN / LINE DRESSING VECTOR + ! + ! idx_dress : position of the basis function used to use the Dressing_vec (usually the largest coeff) + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! u_in : guess coefficients on the various states. Overwritten on exit + ! + ! sze : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! N_st_diag : Number of states in which H is diagonalized. Assumed > sze + ! + ! Initial guess vectors are not necessarily orthonormal + ! + ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) + END_DOC + implicit none + integer, intent(in) :: sze, N_st, N_st_diag, idx_dress + double precision, intent(in) :: H_jj(sze),Dress_jj(sze),Dressing_vec(sze,N_st) + double precision, intent(inout) :: u_in(sze,N_st_diag) + double precision, intent(out) :: energies(N_st_diag) + logical, intent(out) :: converged + external hcalc + + double precision, allocatable :: H_jj_tmp(:) + ASSERT (N_st > 0) + ASSERT (sze > 0) + allocate(H_jj_tmp(sze)) + + do i=1,sze + H_jj_tmp(i) = H_jj(i) + Dress_jj(i) + enddo + do k=1,N_st + do i=1,sze + H_jj_tmp(i) += u_in(i,k) * Dressing_vec(i,k) + enddo + enddo + + integer :: iter + integer :: i,j,k,l,m + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: s_tmp(:,:) + double precision, allocatable :: residual_norm(:),inv_c_idx_dress_vec(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st),inv_c_idx_dress + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + logical :: state_ok(N_st_diag*davidson_sze_max) + integer :: nproc_target + integer :: order(N_st_diag) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag*davidson_sze_max) + + + allocate(inv_c_idx_dress_vec(N_st)) + inv_c_idx_dress = 1.d0/u_in(idx_dress,1) + do i = 1, N_st + inv_c_idx_dress_vec(i) = 1.d0/u_in(idx_dress,i) + enddo + include 'constants.include.F' + + integer :: N_st_diag_in + N_st_diag_in = N_st_diag + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag_in*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag_in*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag_in))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag_in*itermax, N_st_diag_in*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag_in*itermax) &! U + + 1.0d0*dble(sze*m)*(N_st_diag_in*itermax) &! W + + 3.0d0*(N_st_diag_in*itermax)**2 &! h,y,s_tmp + + 1.d0*(N_st_diag_in*itermax) &! lambda + + 1.d0*(N_st_diag_in) &! residual_norm + ! In H_u_0_nstates_zmq + + 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on collector + + 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on slave + + 0.5d0*maxab &! idx0 in H_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag_in,'Number of states in diagonalization') + call write_int(6,sze,'Number of basis functions ') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag_in*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag_in*itermax), & + + ! Small + h(N_st_diag_in*itermax,N_st_diag_in*itermax), & + y(N_st_diag_in*itermax,N_st_diag_in*itermax), & + s_tmp(N_st_diag_in*itermax,N_st_diag_in*itermax), & + residual_norm(N_st_diag_in), & + lambda(N_st_diag_in*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + s_tmp = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag_in >= N_st) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + ! Davidson iterations + ! =================== + + converged = .False. + + do k=N_st+1,N_st_diag_in + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) + enddo + u_in(k,k) = u_in(k,k) + 10.d0 + enddo + do k=1,N_st_diag_in + call normalize(u_in(1,k),sze) + enddo + + do k=1,N_st_diag_in + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 8) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag_in*(iter-1) + shift2 = N_st_diag_in*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag_in,sze) + ! Compute then the DIAGONAL PART OF THE DRESSING + ! += Dress_jj(i) * + call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze) + else + ! Already computed in update below + continue + endif + + + if (N_st == 1) then + + l = idx_dress + double precision :: f + f = inv_c_idx_dress + do istate=1,N_st_diag_in + do i=1,sze + W(i,shift+istate) += Dressing_vec(i,1) *f * U(l,shift+istate) + W(l,shift+istate) += Dressing_vec(i,1) *f * U(i,shift+istate) + enddo + enddo + + else + print*,'dav_double_dressed routine not yet implemented for N_st > 1' +! +! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, & +! psi_coef, size(psi_coef,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, & +! Dressing_vec, size(Dressing_vec,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! +! +! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, & +! Dressing_vec, size(Dressing_vec,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, & +! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), U, size(U,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + ! Diagonalize h + ! --------------- + + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif + + ! Compute Energy for each eigenvector + ! ----------------------------------- + + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, h, size(h,1), y, size(y,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + 0.d0, h, size(h,1)) + + do k=1,shift2 + lambda(k) = h(k,k) + enddo + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag_in + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag_in, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag_in, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag_in + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj_tmp(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + do k=1,N_st + if (residual_norm(k) > 1.d8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + ! Re-contract U and update W + ! -------------------------------- + + call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag_in + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + + do k=1,N_st_diag_in + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + enddo + + + call nullify_small_elements(sze,N_st_diag_in,U,size(U,1),threshold_davidson_pt2) + do k=1,N_st_diag_in + do i=1,sze + u_in(i,k) = U(i,k) + enddo + enddo + + do k=1,N_st_diag_in + energies(k) = lambda(k) + enddo + write_buffer = '======' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + call write_time(6) + + deallocate(W) + + deallocate ( & + residual_norm, & + U, overlap, & + h, y, s_tmp, & + lambda & + ) + FREE nthreads_davidson +end + + +subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze) + implicit none + BEGIN_DOC + ! Routine that computes the diagonal part of the dressing + ! + ! v(i) += u(i) * dress_diag(i) + ! + ! !!!!!!!! WARNING !!!!!!!! the vector v is not initialized + ! + ! !!!!!!!! SO MAKE SURE THERE ARE SOME MEANINGFUL VALUES IN THERE + END_DOC + integer, intent(in) :: N_st,sze + double precision, intent(in) :: u(sze,N_st),dress_diag(sze) + double precision, intent(inout) :: v(sze,N_st) + integer :: i,istate + do istate = 1, N_st + do i = 1, sze + v(i,istate) += dress_diag(i) * u(i,istate) + enddo + enddo +end + + + + + + + + + + + diff --git a/src/dav_general_mat/dav_ext_rout.irp.f b/src/dav_general_mat/dav_ext_rout.irp.f index ff2167e4..0f2d7680 100644 --- a/src/dav_general_mat/dav_ext_rout.irp.f +++ b/src/dav_general_mat/dav_ext_rout.irp.f @@ -3,7 +3,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co use mmap_module implicit none BEGIN_DOC - ! Davidson diagonalization with specific diagonal elements of the H matrix + ! Generic Davidson diagonalization ! ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson ! diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index a63eb4a3..2a655eed 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -38,7 +38,15 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer, external :: omp_get_num_threads - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) + if (omp_get_num_threads() == 1) then + !$OMP PARALLEL DEFAULT(SHARED) + !$OMP SINGLE + call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) + !$OMP END SINGLE + !$OMP END PARALLEL + else + call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) + endif end recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level) @@ -49,7 +57,7 @@ BEGIN_TEMPLATE $type :: c, tmp integer :: itmp integer :: i, j - + if(isize<2)return c = x( shiftr(first+last,1) ) @@ -81,11 +89,16 @@ BEGIN_TEMPLATE endif else if (first < i-1) then + !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,first,i,level) call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2) + !$OMP END TASK endif if (j+1 < last) then + !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,last,j,level) call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2) + !$OMP END TASK endif + !$OMP TASKWAIT endif end @@ -249,60 +262,7 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE - -!---------------------- INTEL -IRP_IF INTEL - BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - character, allocatable :: tmp(:) - if (isize < 2) return - call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) - allocate(tmp(n)) - call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp) - deallocate(tmp) - iorder(1:isize) = iorder(1:isize)+1 - call $Xset_order(x,iorder,isize) - end - - subroutine $Xsort_noidx(x,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer :: n - character, allocatable :: tmp(:) - if (isize < 2) return - call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) - allocate(tmp(n)) - call ippsSortRadixAscend_$ityp_I(x, isize, tmp) - deallocate(tmp) - end - -SUBST [ X, type, ityp, n, ippsz ] - ; real ; 32f ; 4 ; 13 ;; - i ; integer ; 32s ; 4 ; 11 ;; - i2 ; integer*2 ; 16s ; 2 ; 7 ;; -END_TEMPLATE - -BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -329,12 +289,12 @@ BEGIN_TEMPLATE endif end subroutine $Xsort -SUBST [ X, type ] - d ; double precision ;; +SUBST [ X, type, Y ] + ; real ; i ;; + d ; double precision ; i8 ;; END_TEMPLATE BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -346,112 +306,8 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer :: n - if (isize < 2) then - return - endif - call sorted_$Xnumber(x,isize,n) - if (isize == n) then - return - endif - if ( isize < 32) then - call insertion_$Xsort(x,iorder,isize) - else - call $Xradix_sort(x,iorder,isize,-1) - endif - end subroutine $Xsort - -SUBST [ X, type ] - i8 ; integer*8 ;; -END_TEMPLATE - -!---------------------- END INTEL -IRP_ELSE -!---------------------- NON-INTEL -BEGIN_TEMPLATE - - subroutine $Xsort_noidx(x,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer, allocatable :: iorder(:) - integer :: i - allocate(iorder(isize)) - do i=1,isize - iorder(i)=i - enddo - call $Xsort(x,iorder,isize) - deallocate(iorder) - end subroutine $Xsort_noidx - -SUBST [ X, type ] - ; real ;; - d ; double precision ;; - i ; integer ;; - i8 ; integer*8 ;; - i2 ; integer*2 ;; -END_TEMPLATE - -BEGIN_TEMPLATE - - subroutine $Xsort(x,iorder,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - if (isize < 2) then - return - endif -! call sorted_$Xnumber(x,isize,n) -! if (isize == n) then -! return -! endif - if ( isize < 32) then - call insertion_$Xsort(x,iorder,isize) - else -! call heap_$Xsort(x,iorder,isize) - call quick_$Xsort(x,iorder,isize) - endif - end subroutine $Xsort - -SUBST [ X, type ] - ; real ;; - d ; double precision ;; -END_TEMPLATE - -BEGIN_TEMPLATE - - subroutine $Xsort(x,iorder,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - if (isize < 2) then - return - endif - call sorted_$Xnumber(x,isize,n) - if (isize == n) then - return - endif - if ( isize < 32) then - call insertion_$Xsort(x,iorder,isize) - else - call $Xradix_sort(x,iorder,isize,-1) - endif +! call $Xradix_sort(x,iorder,isize,-1) + call quick_$Xsort(x,iorder,isize) end subroutine $Xsort SUBST [ X, type ] @@ -460,11 +316,6 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE -IRP_ENDIF -!---------------------- END NON-INTEL - - - BEGIN_TEMPLATE subroutine $Xset_order(x,iorder,isize) implicit none @@ -562,12 +413,10 @@ SUBST [ X, type ] i2; integer*2 ;; END_TEMPLATE - BEGIN_TEMPLATE -recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) + recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) implicit none - BEGIN_DOC ! Sort integer array x(isize) using the radix sort algorithm. ! iorder in input should be (1,2,3,...,isize), and in output @@ -703,14 +552,24 @@ recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) endif +! !$OMP PARALLEL DEFAULT(SHARED) if (isize > 1000000) +! !$OMP SINGLE if (i3>1_$int_type) then +! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(i3 > 1000000) call $Xradix_sort$big(x,iorder,i3,iradix_new-1) +! !$OMP END TASK endif if (isize-i3>1_$int_type) then +! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(isize-i3 > 1000000) call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1) +! !$OMP END TASK endif +! !$OMP TASKWAIT +! !$OMP END SINGLE +! !$OMP END PARALLEL + return endif @@ -765,11 +624,16 @@ recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) if (i1>1_$int_type) then + !$OMP TASK FIRSTPRIVATE(i0,iradix,i1) SHARED(x,iorder) if(i1 >1000000) call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1) + !$OMP END TASK endif if (i0>1) then + !$OMP TASK FIRSTPRIVATE(i0,iradix) SHARED(x,iorder) if(i0 >1000000) call $Xradix_sort$big(x,iorder,i0,iradix-1) + !$OMP END TASK endif + !$OMP TASKWAIT end @@ -782,4 +646,3 @@ SUBST [ X, type, integer_size, is_big, big, int_type ] END_TEMPLATE - From fccd7e2d1a22c3e1599311cd431bed602314cc71 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 26 Oct 2021 20:49:48 +0200 Subject: [PATCH 09/20] reput the sort.irp.f --- src/utils/sort.irp.f | 209 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 173 insertions(+), 36 deletions(-) diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index 2a655eed..a63eb4a3 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -38,15 +38,7 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer, external :: omp_get_num_threads - if (omp_get_num_threads() == 1) then - !$OMP PARALLEL DEFAULT(SHARED) - !$OMP SINGLE - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) - !$OMP END SINGLE - !$OMP END PARALLEL - else - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) - endif + call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) end recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level) @@ -57,7 +49,7 @@ BEGIN_TEMPLATE $type :: c, tmp integer :: itmp integer :: i, j - + if(isize<2)return c = x( shiftr(first+last,1) ) @@ -89,16 +81,11 @@ BEGIN_TEMPLATE endif else if (first < i-1) then - !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,first,i,level) call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2) - !$OMP END TASK endif if (j+1 < last) then - !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,last,j,level) call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2) - !$OMP END TASK endif - !$OMP TASKWAIT endif end @@ -262,7 +249,60 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE + +!---------------------- INTEL +IRP_IF INTEL + BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp) + deallocate(tmp) + iorder(1:isize) = iorder(1:isize)+1 + call $Xset_order(x,iorder,isize) + end + + subroutine $Xsort_noidx(x,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer :: n + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixAscend_$ityp_I(x, isize, tmp) + deallocate(tmp) + end + +SUBST [ X, type, ityp, n, ippsz ] + ; real ; 32f ; 4 ; 13 ;; + i ; integer ; 32s ; 4 ; 11 ;; + i2 ; integer*2 ; 16s ; 2 ; 7 ;; +END_TEMPLATE + +BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -289,12 +329,12 @@ BEGIN_TEMPLATE endif end subroutine $Xsort -SUBST [ X, type, Y ] - ; real ; i ;; - d ; double precision ; i8 ;; +SUBST [ X, type ] + d ; double precision ;; END_TEMPLATE BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -306,8 +346,112 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer :: n -! call $Xradix_sort(x,iorder,isize,-1) - call quick_$Xsort(x,iorder,isize) + if (isize < 2) then + return + endif + call sorted_$Xnumber(x,isize,n) + if (isize == n) then + return + endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else + call $Xradix_sort(x,iorder,isize,-1) + endif + end subroutine $Xsort + +SUBST [ X, type ] + i8 ; integer*8 ;; +END_TEMPLATE + +!---------------------- END INTEL +IRP_ELSE +!---------------------- NON-INTEL +BEGIN_TEMPLATE + + subroutine $Xsort_noidx(x,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer, allocatable :: iorder(:) + integer :: i + allocate(iorder(isize)) + do i=1,isize + iorder(i)=i + enddo + call $Xsort(x,iorder,isize) + deallocate(iorder) + end subroutine $Xsort_noidx + +SUBST [ X, type ] + ; real ;; + d ; double precision ;; + i ; integer ;; + i8 ; integer*8 ;; + i2 ; integer*2 ;; +END_TEMPLATE + +BEGIN_TEMPLATE + + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + if (isize < 2) then + return + endif +! call sorted_$Xnumber(x,isize,n) +! if (isize == n) then +! return +! endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else +! call heap_$Xsort(x,iorder,isize) + call quick_$Xsort(x,iorder,isize) + endif + end subroutine $Xsort + +SUBST [ X, type ] + ; real ;; + d ; double precision ;; +END_TEMPLATE + +BEGIN_TEMPLATE + + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + if (isize < 2) then + return + endif + call sorted_$Xnumber(x,isize,n) + if (isize == n) then + return + endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else + call $Xradix_sort(x,iorder,isize,-1) + endif end subroutine $Xsort SUBST [ X, type ] @@ -316,6 +460,11 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE +IRP_ENDIF +!---------------------- END NON-INTEL + + + BEGIN_TEMPLATE subroutine $Xset_order(x,iorder,isize) implicit none @@ -413,10 +562,12 @@ SUBST [ X, type ] i2; integer*2 ;; END_TEMPLATE + BEGIN_TEMPLATE - recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) +recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) implicit none + BEGIN_DOC ! Sort integer array x(isize) using the radix sort algorithm. ! iorder in input should be (1,2,3,...,isize), and in output @@ -552,24 +703,14 @@ BEGIN_TEMPLATE endif -! !$OMP PARALLEL DEFAULT(SHARED) if (isize > 1000000) -! !$OMP SINGLE if (i3>1_$int_type) then -! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(i3 > 1000000) call $Xradix_sort$big(x,iorder,i3,iradix_new-1) -! !$OMP END TASK endif if (isize-i3>1_$int_type) then -! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(isize-i3 > 1000000) call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1) -! !$OMP END TASK endif -! !$OMP TASKWAIT -! !$OMP END SINGLE -! !$OMP END PARALLEL - return endif @@ -624,16 +765,11 @@ BEGIN_TEMPLATE if (i1>1_$int_type) then - !$OMP TASK FIRSTPRIVATE(i0,iradix,i1) SHARED(x,iorder) if(i1 >1000000) call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1) - !$OMP END TASK endif if (i0>1) then - !$OMP TASK FIRSTPRIVATE(i0,iradix) SHARED(x,iorder) if(i0 >1000000) call $Xradix_sort$big(x,iorder,i0,iradix-1) - !$OMP END TASK endif - !$OMP TASKWAIT end @@ -646,3 +782,4 @@ SUBST [ X, type, integer_size, is_big, big, int_type ] END_TEMPLATE + From b19569905057f36aedbd478620d57cf4f8b4542b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 28 Oct 2021 15:19:45 +0200 Subject: [PATCH 10/20] added H_matrix_diag_all_dets --- src/determinants/utils.irp.f | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/determinants/utils.irp.f b/src/determinants/utils.irp.f index 3aec16f9..bd2ee293 100644 --- a/src/determinants/utils.irp.f +++ b/src/determinants/utils.irp.f @@ -20,6 +20,25 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] !$OMP END PARALLEL DO END_PROVIDER +BEGIN_PROVIDER [ double precision, H_matrix_diag_all_dets,(N_det) ] + use bitmasks + implicit none + BEGIN_DOC + ! |H| matrix on the basis of the Slater determinants defined by psi_det + END_DOC + integer :: i + double precision :: hij + integer :: degree(N_det) + call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,hij,degree) & + !$OMP SHARED (N_det, psi_det, N_int,H_matrix_diag_all_dets) + do i =1,N_det + call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij) + H_matrix_diag_all_dets(i) = hij + enddo + !$OMP END PARALLEL DO +END_PROVIDER + BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ] use bitmasks From 18bb5ee9172841bd4c8a33ca144f65be2b469a22 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 29 Oct 2021 11:44:37 +0200 Subject: [PATCH 11/20] added davidson for diagonal dress --- .../dav_diag_dressed_ext_rout.irp.f | 481 ++++++++++++++++++ .../dav_double_dress_ext_rout.irp.f | 2 +- .../dav_dressed_ext_rout.irp.f | 4 +- src/dav_general_mat/dav_ext_rout.irp.f | 8 +- src/determinants/utils.irp.f | 2 + 5 files changed, 491 insertions(+), 6 deletions(-) create mode 100644 src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f diff --git a/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f new file mode 100644 index 00000000..2f3d7f80 --- /dev/null +++ b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f @@ -0,0 +1,481 @@ + +subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) + use mmap_module + implicit none + BEGIN_DOC + ! Generic Davidson diagonalization with ONE DIAGONAL DRESSING OPERATOR + ! + ! Dress_jj : DIAGONAL DRESSING of the Hamiltonian + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! u_in : guess coefficients on the various states. Overwritten on exit + ! + ! sze : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze + ! + ! Initial guess vectors are not necessarily orthonormal + ! + ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) + END_DOC + integer, intent(in) :: sze, N_st, N_st_diag_in + double precision, intent(in) :: H_jj(sze),Dress_jj(sze) + double precision, intent(inout) :: u_in(sze,N_st_diag_in) + double precision, intent(out) :: energies(N_st) + external hcalc + + integer :: iter, N_st_diag + integer :: i,j,k,l,m + logical, intent(inout) :: converged + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: residual_norm(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st) + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + integer :: nproc_target + integer :: order(N_st_diag_in) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:)!, S_d(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag_in*davidson_sze_max) + + include 'constants.include.F' + + N_st_diag = N_st_diag_in + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag*itermax, N_st_diag*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + provide threshold_davidson !nthreads_davidson + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab + maxab = sze + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag*itermax) &! U + + 1.d0*dble(sze*m)*(N_st_diag*itermax) &! W + + 2.0d0*(N_st_diag*itermax)**2 &! h,y + + 2.d0*(N_st_diag*itermax) &! s2,lambda + + 1.d0*(N_st_diag) &! residual_norm + ! In H_S2_u_0_nstates_zmq + + 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on collector + + 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on slave + + 0.5d0*maxab &! idx0 in H_S2_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag,'Number of states in diagonalization') + call write_int(6,sze,'Number of basis functions') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + double precision, allocatable :: H_jj_tmp(:) + ASSERT (N_st > 0) + ASSERT (sze > 0) + allocate(H_jj_tmp(sze)) + + do i=1,sze + H_jj_tmp(i) = H_jj(i) + Dress_jj(i) + enddo + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag*itermax), & + ! Small + h(N_st_diag*itermax,N_st_diag*itermax), & + y(N_st_diag*itermax,N_st_diag*itermax), & + residual_norm(N_st_diag), & + lambda(N_st_diag*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag >= N_st) + ASSERT (sze > 0) + + ! Davidson iterations + ! =================== + + converged = .False. + + ! Initialize from N_st to N_st_diat with gaussian random numbers + ! to be sure to have overlap with any eigenvectors + do k=N_st+1,N_st_diag + u_in(k,k) = 10.d0 + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) + enddo + enddo + ! Normalize all states + do k=1,N_st_diag + call normalize(u_in(1,k),sze) + enddo + + ! Copy from the guess input "u_in" to the working vectors "U" + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 8) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag*(iter-1) + shift2 = N_st_diag*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + + ! Gram-Schmidt to orthogonalize all new guess with the previous vectors + call ortho_qr(U,size(U,1),sze,shift2) + call ortho_qr(U,size(U,1),sze,shift2) + ! it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag) + ! where sze is the size of the vector, N_st_diag is the number of states + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze) + ! Compute then the DIAGONAL PART OF THE DRESSING + ! += Dress_jj(i) * + call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze) + else + ! Already computed in update below + continue + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + + ! Diagonalize h y = lambda y + ! --------------- + + call lapack_diag(lambda,y,h,size(h,1),shift2) + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + do k=1,N_st + l = order(k) + if (k /= l) then + lambda(k) = overlap(l,1) + endif + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj_tmp(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + + do k=1,N_st + if (residual_norm(k) > 1.e8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) + do j=1,N_st_diag + k=1 + do while ((k 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- - ! Gram-Schmidt to orthogonalize all new guess with the previous vectors call ortho_qr(U,size(U,1),sze,shift2) call ortho_qr(U,size(U,1),sze,shift2) @@ -345,6 +344,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co enddo + ! Re-contract U and update W + ! -------------------------------- + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) do k=1,N_st_diag @@ -360,8 +362,8 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co U(i,k) = u_in(i,k) enddo enddo - call ortho_qr(U,size(U,1),sze,N_st_diag) - call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) do j=1,N_st_diag k=1 do while ((k Date: Fri, 5 Nov 2021 15:17:39 +0100 Subject: [PATCH 12/20] removed stupid assert --- .../dav_double_dress_ext_rout.irp.f | 2 -- src/mo_two_e_ints/map_integrals.irp.f | 30 +++++++++---------- 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/src/dav_general_mat/dav_double_dress_ext_rout.irp.f b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f index 4f4572fa..884fd672 100644 --- a/src/dav_general_mat/dav_double_dress_ext_rout.irp.f +++ b/src/dav_general_mat/dav_double_dress_ext_rout.irp.f @@ -207,8 +207,6 @@ subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies ASSERT (N_st > 0) ASSERT (N_st_diag_in >= N_st) ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) ! Davidson iterations ! =================== diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 9f73d518..272916e3 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -302,21 +302,21 @@ end integer(key_kind) :: idx double precision :: tmp - icount = 1 ! Avoid division by zero - do j=1,mo_num - do i=1,j-1 - call two_e_integrals_index(i,j,j,i,idx) - !DIR$ FORCEINLINE - call map_get(mo_integrals_map,idx,tmp) - banned_excitation(i,j) = dabs(tmp) < 1.d-14 - banned_excitation(j,i) = banned_excitation(i,j) - if (banned_excitation(i,j)) icount = icount+2 - enddo - enddo - use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1% - if (use_banned_excitation) then - print *, 'Using sparsity of exchange integrals' - endif +!icount = 1 ! Avoid division by zero +!do j=1,mo_num +! do i=1,j-1 +! call two_e_integrals_index(i,j,j,i,idx) +! !DIR$ FORCEINLINE +! call map_get(mo_integrals_map,idx,tmp) +! banned_excitation(i,j) = dabs(tmp) < 1.d-14 +! banned_excitation(j,i) = banned_excitation(i,j) +! if (banned_excitation(i,j)) icount = icount+2 +! enddo +!enddo +!use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1% +!if (use_banned_excitation) then +! print *, 'Using sparsity of exchange integrals' +!endif END_PROVIDER From af169a4e3cb3502df3ea0946f77f626e9d5a4aad Mon Sep 17 00:00:00 2001 From: FiletoRodriguez Date: Tue, 9 Nov 2021 15:48:49 +0100 Subject: [PATCH 13/20] Added grad_rho --- src/dft_one_e/mu_erf_dft.irp.f | 2 ++ src/dft_utils_func/mu_of_r_dft.irp.f | 2 +- src/kohn_sham_rs/rs_ks_scf.irp.f | 2 +- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index f9bff24d..bf590c84 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -22,6 +22,8 @@ BEGIN_PROVIDER [double precision, mu_of_r_dft, (n_points_final_grid)] mu_of_r_dft(i) = mu_of_r_hf(i) else if(mu_dft_type == "rsc")then mu_of_r_dft(i) = mu_rsc_of_r(i) + else if(mu_dft_type == "grad_rho")then + mu_of_r_dft(i) = mu_grad_rho(i) else print*,'mu_dft_type is not of good type = ',mu_dft_type print*,'it must be of type Read, cst, hf, rsc' diff --git a/src/dft_utils_func/mu_of_r_dft.irp.f b/src/dft_utils_func/mu_of_r_dft.irp.f index 7cba0a60..0e9a0f1b 100644 --- a/src/dft_utils_func/mu_of_r_dft.irp.f +++ b/src/dft_utils_func/mu_of_r_dft.irp.f @@ -17,7 +17,7 @@ double precision function mu_grad_rho_func(r) integer :: m double precision :: rho, dm_a, dm_b, grad_dm_a(3), grad_dm_b(3) double precision :: eta, grad_rho(3), grad_sqr - eta = 0.135d0 + eta = mu_erf call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b) rho = dm_a + dm_b do m = 1,3 diff --git a/src/kohn_sham_rs/rs_ks_scf.irp.f b/src/kohn_sham_rs/rs_ks_scf.irp.f index 5d23544e..84b85136 100644 --- a/src/kohn_sham_rs/rs_ks_scf.irp.f +++ b/src/kohn_sham_rs/rs_ks_scf.irp.f @@ -17,7 +17,7 @@ program rs_ks_scf print*, '**************************' print*, 'mu_erf_dft = ',mu_erf_dft print*, '**************************' - call check_coherence_functional +! call check_coherence_functional call create_guess call orthonormalize_mos call run From 8c521906488177db9dd3695cb49e08aadbd98e31 Mon Sep 17 00:00:00 2001 From: FiletoRodriguez Date: Tue, 30 Nov 2021 15:35:19 +0100 Subject: [PATCH 14/20] added src/tools/save_natorb_no_ref.irp.f --- src/determinants/density_matrix.irp.f | 38 +++++++++++++++++++++++++++ src/tools/save_natorb_no_ref.irp.f | 24 +++++++++++++++++ 2 files changed, 62 insertions(+) create mode 100644 src/tools/save_natorb_no_ref.irp.f diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index 6a474662..fa4b3328 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -268,6 +268,44 @@ subroutine set_natural_mos soft_touch mo_occ end + +subroutine save_natural_mos_canon_label + implicit none + BEGIN_DOC + ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in + ! the |MO| basis + END_DOC + call set_natural_mos_canon_label + call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) + call orthonormalize_mos + call save_mos +end + +subroutine set_natural_mos_canon_label + implicit none + BEGIN_DOC + ! Set natural orbitals, obtained by diagonalization of the one-body density matrix + ! in the |MO| basis + END_DOC + character*(64) :: label + double precision, allocatable :: tmp(:,:) + + label = "Canonical" + integer :: i,j,iorb,jorb + do i = 1, n_virt_orb + iorb = list_virt(i) + do j = 1, n_core_inact_act_orb + jorb = list_core_inact_act(j) + enddo + enddo + call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label) + soft_touch mo_occ +end + + + + + subroutine set_natorb_no_ov_rot implicit none BEGIN_DOC diff --git a/src/tools/save_natorb_no_ref.irp.f b/src/tools/save_natorb_no_ref.irp.f new file mode 100644 index 00000000..9d253fa0 --- /dev/null +++ b/src/tools/save_natorb_no_ref.irp.f @@ -0,0 +1,24 @@ +program save_natorb + implicit none + BEGIN_DOC +! Save natural |MOs| into the |EZFIO|. +! +! This program reads the wave function stored in the |EZFIO| directory, +! extracts the corresponding natural orbitals and setd them as the new +! |MOs|. +! +! If this is a multi-state calculation, the density matrix that produces +! the natural orbitals is obtained from an average of the density +! matrices of each state with the corresponding +! :option:`determinants state_average_weight` + END_DOC + read_wf = .True. + touch read_wf + call save_natural_mos_canon_label + call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('None') + call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None') +end + From 9d0e7c7034e0d0089fdd29ab68c61cb160cb733e Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 17 Dec 2021 18:15:47 +0100 Subject: [PATCH 15/20] moved a subroutine from determinant to bitmask --- src/bitmask/bitmasks_routines.irp.f | 15 ++++++++++++++ src/determinants/slater_rules.irp.f | 31 ++++++++++++++++++++++++++++- 2 files changed, 45 insertions(+), 1 deletion(-) diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index c34d54dc..9c6f4f0c 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -268,6 +268,21 @@ subroutine print_spindet(string,Nint) end +subroutine print_det_one_dimension(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant using the '+-' notation + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint) + character*(2048) :: output(1) + + call bitstring_to_str( output(1), string, Nint ) + print *, trim(output(1)) + +end + logical function is_integer_in_string(bite,string,Nint) use bitmasks implicit none diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 04cf861f..3a33a37d 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -438,7 +438,7 @@ subroutine bitstring_to_list_ab( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string ! For alpha/beta determinants. END_DOC integer, intent(in) :: Nint @@ -472,6 +472,35 @@ subroutine bitstring_to_list_ab( string, list, n_elements, Nint) end +!subroutine bitstring_to_list( string, list, n_elements, Nint) +! use bitmasks +! implicit none +! BEGIN_DOC +! ! Gives the indices(+1) of the bits set to 1 in the bit string +! END_DOC +! integer, intent(in) :: Nint +! integer(bit_kind), intent(in) :: string(Nint) +! integer, intent(out) :: list(Nint*bit_kind_size) +! integer, intent(out) :: n_elements +! +! integer :: i, j, ishift +! integer(bit_kind) :: l +! +! n_elements = 0 +! ishift = 1 +! do i=1,Nint +! l = string(i) +! do while (l /= 0_bit_kind) +! j = trailz(l) +! n_elements = n_elements + 1 +! l = ibclr(l,j) +! list(n_elements) = ishift+j +! enddo +! ishift = ishift + bit_kind_size +! enddo +! +!end + subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) use bitmasks From 7fba7a120907319fb346df605e07100763596501 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 23 Dec 2021 15:14:54 +0100 Subject: [PATCH 16/20] Removed qmkl --- config/ifort_2021_avx.cfg | 2 +- config/ifort_2021_avx_mpi.cfg | 2 +- config/ifort_2021_sse4.cfg | 2 +- config/ifort_2021_sse4_mpi.cfg | 2 +- config/ifort_2021_xHost.cfg | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/config/ifort_2021_avx.cfg b/config/ifort_2021_avx.cfg index 8fadda67..6f657052 100644 --- a/config/ifort_2021_avx.cfg +++ b/config/ifort_2021_avx.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_avx_mpi.cfg b/config/ifort_2021_avx_mpi.cfg index b6b74b73..c991a4a9 100644 --- a/config/ifort_2021_avx_mpi.cfg +++ b/config/ifort_2021_avx_mpi.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL diff --git a/config/ifort_2021_sse4.cfg b/config/ifort_2021_sse4.cfg index 269023da..a6299665 100644 --- a/config/ifort_2021_sse4.cfg +++ b/config/ifort_2021_sse4.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DINTEL diff --git a/config/ifort_2021_sse4_mpi.cfg b/config/ifort_2021_sse4_mpi.cfg index 41df87bc..6ae56d2a 100644 --- a/config/ifort_2021_sse4_mpi.cfg +++ b/config/ifort_2021_sse4_mpi.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL diff --git a/config/ifort_2021_xHost.cfg b/config/ifort_2021_xHost.cfg index 05c271f3..1e76a69d 100644 --- a/config/ifort_2021_xHost.cfg +++ b/config/ifort_2021_xHost.cfg @@ -7,7 +7,7 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -qmkl=parallel -lirc -lsvml -limf -lipps +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=64 -DINTEL From 2a98378fcab6454c060e3cd4ad543506c16e2103 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 30 Dec 2021 14:43:14 +0100 Subject: [PATCH 17/20] Fix basis bugs when working on trexio --- ocaml/qp_create_ezfio.ml | 60 +++++++++++++++++++++++----------------- src/pseudo/EZFIO.cfg | 6 ++-- 2 files changed, 38 insertions(+), 28 deletions(-) diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index 8fcaf5fc..be6c305b 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -91,7 +91,7 @@ let run ?o b au c d m p cart xyz_file = | Element e -> Element.to_string e | Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i) (Element.to_string e) in - Hashtbl.find basis_table key + Hashtbl.find basis_table key in let temp_filename = @@ -132,7 +132,7 @@ let run ?o b au c d m p cart xyz_file = Element.to_string elem.Atom.element in Hashtbl.add basis_table key new_channel - ) nuclei + ) nuclei end | Some (key, basis) -> (*Aux basis *) begin @@ -277,6 +277,16 @@ let run ?o b au c d m p cart xyz_file = ) nuclei in + let z_core = + List.map (fun x -> + Positive_int.to_int x.Pseudo.n_elec + |> float_of_int + ) pseudo + in + let nucl_num = (List.length z_core) in + Ezfio.set_pseudo_nucl_charge_remove (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:z_core); + let molecule = let n_elec_to_remove = List.fold_left (fun accu x -> @@ -293,13 +303,13 @@ let run ?o b au c d m p cart xyz_file = Molecule.nuclei = let charges = list_map (fun x -> Positive_int.to_int x.Pseudo.n_elec - |> Float.of_int) pseudo + |> Float.of_int) pseudo |> Array.of_list in List.mapi (fun i x -> { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) |> Charge.of_float } - ) molecule.Molecule.nuclei + ) molecule.Molecule.nuclei } in let nuclei = @@ -356,11 +366,11 @@ let run ?o b au c d m p cart xyz_file = in if (x > accu) then x else accu - ) 0 x.Pseudo.non_local + ) 0 x.Pseudo.non_local in if (x > accu) then x else accu - ) 0 pseudo + ) 0 pseudo in let kmax = @@ -368,10 +378,10 @@ let run ?o b au c d m p cart xyz_file = list_map (fun x -> List.filter (fun (y,_) -> (Positive_int.to_int y.Pseudo.GaussianPrimitive_non_local.proj) = i) - x.Pseudo.non_local - |> List.length ) pseudo + x.Pseudo.non_local + |> List.length ) pseudo |> List.fold_left (fun accu x -> - if accu > x then accu else x) 0 + if accu > x then accu else x) 0 ) |> Array.fold_left (fun accu i -> if i > accu then i else accu) 0 @@ -396,11 +406,11 @@ let run ?o b au c d m p cart xyz_file = in tmp_array_dz_k.(i).(j) <- y; tmp_array_n_k.(i).(j) <- z; - ) x.Pseudo.local + ) x.Pseudo.local ) pseudo ; let concat_2d tmp_array = let data = - Array.map Array.to_list tmp_array + Array.map Array.to_list tmp_array |> Array.to_list |> List.concat in @@ -438,14 +448,14 @@ let run ?o b au c d m p cart xyz_file = tmp_array_dz_kl.(k).(i).(j) <- y; tmp_array_n_kl.(k).(i).(j) <- z; last_idx.(k) <- i+1; - ) x.Pseudo.non_local + ) x.Pseudo.non_local ) pseudo ; let concat_3d tmp_array = let data = Array.map (fun x -> Array.map Array.to_list x |> Array.to_list - |> List.concat) tmp_array + |> List.concat) tmp_array |> Array.to_list |> List.concat in @@ -513,8 +523,8 @@ let run ?o b au c d m p cart xyz_file = Ezfio.set_ao_basis_ao_num ao_num; Ezfio.set_ao_basis_ao_basis b; Ezfio.set_basis_basis b; - let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis - and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis + let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis + and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis and ao_power= let l = list_map (fun (x,_,_) -> x) long_basis in (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.x)) l)@ @@ -526,7 +536,7 @@ let run ?o b au c d m p cart xyz_file = else s) 0 ao_prim_num in let gtos = - list_map (fun (_,x,_) -> x) long_basis + list_map (fun (_,x,_) -> x) long_basis in let create_expo_coef ec = @@ -534,10 +544,10 @@ let run ?o b au c d m p cart xyz_file = begin match ec with | `Coefs -> list_map (fun x-> list_map (fun (_,coef) -> - AO_coef.to_float coef) x.Gto.lc) gtos + AO_coef.to_float coef) x.Gto.lc) gtos | `Expos -> list_map (fun x-> list_map (fun (prim,_) -> AO_expo.to_float - prim.GaussianPrimitive.expo) x.Gto.lc) gtos + prim.GaussianPrimitive.expo) x.Gto.lc) gtos end in let rec get_n n accu = function @@ -567,7 +577,7 @@ let run ?o b au c d m p cart xyz_file = list_map ( fun (g,_) -> g.Gto.lc ) basis in let ang_mom = - list_map (fun (l : (GaussianPrimitive.t * Qptypes.AO_coef.t) list) -> + list_map (fun (l : (GaussianPrimitive.t * Qptypes.AO_coef.t) list) -> let x, _ = List.hd l in Angmom.to_l x.GaussianPrimitive.sym |> Qptypes.Positive_int.to_int ) lc @@ -577,7 +587,7 @@ let run ?o b au c d m p cart xyz_file = |> List.concat in let coef = - list_map (fun l -> + list_map (fun l -> list_map (fun (_,x) -> Qptypes.AO_coef.to_float x) l ) lc |> List.concat @@ -610,14 +620,14 @@ let run ?o b au c d m p cart xyz_file = ~rank:1 ~dim:[| prim_num |] ~data:shell_idx) ; Ezfio.set_basis_basis_nucleus_index (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| shell_num |] - ~data:( list_map (fun (_,n) -> Nucl_number.to_int n) basis) + ~data:( list_map (fun (_,n) -> Nucl_number.to_int n) basis) ) ; Ezfio.set_basis_nucleus_shell_num(Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| nucl_num |] ~data:( list_map (fun (_,n) -> Nucl_number.to_int n) basis - |> List.fold_left (fun accu i -> - match accu with + |> List.fold_left (fun accu i -> + match accu with | [] -> [(1,i)] | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((1,i)::(h,j)::rest) ) [] @@ -713,7 +723,7 @@ If a file with the same name as the basis set exists, this file will be read. O anonymous "FILE" Mandatory "Input file in xyz format or z-matrix."; ] - |> set_specs + |> set_specs end; @@ -737,7 +747,7 @@ If a file with the same name as the basis set exists, this file will be read. O | None -> 0 | Some x -> ( if x.[0] = 'm' then ~- (int_of_string (String.sub x 1 (String.length x - 1))) - else + else int_of_string x ) in diff --git a/src/pseudo/EZFIO.cfg b/src/pseudo/EZFIO.cfg index 952ebe04..2b76a5c3 100644 --- a/src/pseudo/EZFIO.cfg +++ b/src/pseudo/EZFIO.cfg @@ -10,7 +10,7 @@ type:integer interface: ezfio,provider [pseudo_n_k] -doc: Number of gaussians in the local component +doc: Powers of r - 2 in the local component type: integer interface: ezfio,provider size: (nuclei.nucl_num,pseudo.pseudo_klocmax) @@ -38,7 +38,7 @@ type:integer interface: ezfio,provider [pseudo_n_kl] -doc: Number of functions in the non-local component +doc: Powers of r - 2 in the non-local component type: integer interface: ezfio,provider size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax) @@ -69,7 +69,7 @@ default: 1000 [pseudo_grid_rmax] type: double precision -doc: R_max of the QMC grid +doc: R_max of the QMC grid interface: ezfio,provider,ocaml default: 10.0 From bc71a23ffb963234fe9b5652a248a832dc4c56bf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 4 Jan 2022 12:53:34 +0100 Subject: [PATCH 18/20] Bug in molden + bug in print --- src/cipsi/pt2_stoch_routines.irp.f | 2 +- src/tools/molden.irp.f | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 3fa2641a..14b1d060 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -131,7 +131,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_det_hii selection_weight pseudo_sym - PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max + PROVIDE list_act list_inact list_core list_virt list_del seniority_max PROVIDE excitation_beta_max excitation_alpha_max excitation_max if (h0_type == 'CFG') then diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f index 417b25ad..830a141e 100644 --- a/src/tools/molden.irp.f +++ b/src/tools/molden.irp.f @@ -52,8 +52,8 @@ program molden l += 1 if (l > ao_num) exit enddo - write(i_unit_output,*)'' enddo + write(i_unit_output,*)'' enddo From 1ee940018b8d778b3cbab3cbc4f849124d9f63cd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 7 Jan 2022 15:08:51 +0100 Subject: [PATCH 19/20] Added ifort_2019_mpi_rome.cfg --- config/ifort_2019_mpi_rome.cfg | 63 ++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 config/ifort_2019_mpi_rome.cfg diff --git a/config/ifort_2019_mpi_rome.cfg b/config/ifort_2019_mpi_rome.cfg new file mode 100644 index 00000000..171219e6 --- /dev/null +++ b/config/ifort_2019_mpi_rome.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : mpiifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback -shared-intel +FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + From cf3f510704b76cd23cc463aa9237327e1b82311d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 7 Jan 2022 15:11:05 +0100 Subject: [PATCH 20/20] Added ifort_2021_mpi_rome.cfg --- config/ifort_2021_mpi_rome.cfg | 63 ++++++++++++++++++++++++++++++++++ config/ifort_2021_rome.cfg | 63 ++++++++++++++++++++++++++++++++++ 2 files changed, 126 insertions(+) create mode 100644 config/ifort_2021_mpi_rome.cfg create mode 100644 config/ifort_2021_rome.cfg diff --git a/config/ifort_2021_mpi_rome.cfg b/config/ifort_2021_mpi_rome.cfg new file mode 100644 index 00000000..8413d23d --- /dev/null +++ b/config/ifort_2021_mpi_rome.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : mpiifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback -shared-intel +FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_2021_rome.cfg b/config/ifort_2021_rome.cfg new file mode 100644 index 00000000..b3023186 --- /dev/null +++ b/config/ifort_2021_rome.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : ifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback -shared-intel +FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp +