From e85a5927a1910a9e7a9946306ce64d214308b8b2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 8 Dec 2015 13:24:43 +0100 Subject: [PATCH] Natural orbitals using SVD --- config/ifort.cfg | 4 +-- install/Downloads/.gitignore | 0 ocaml/.gitignore | 1 + plugins/CAS_SD/cas_s_selected.irp.f | 2 +- scripts/compilation/qp_create_ninja.py | 2 +- src/Determinants/density_matrix.irp.f | 14 +++----- src/MO_Basis/mo_overlap.irp.f | 2 +- src/MO_Basis/utils.irp.f | 47 ++++++++++++++++++++++++++ src/Utils/LinearAlgebra.irp.f | 7 ++-- 9 files changed, 61 insertions(+), 18 deletions(-) delete mode 100644 install/Downloads/.gitignore diff --git a/config/ifort.cfg b/config/ifort.cfg index c1d7e968..bfa41c03 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -31,14 +31,14 @@ OPENMP : 1 ; Append OpenMP flags # -ftz : Flushes denormal results to zero # [OPT] -FCFLAGS : -xSSE4.2 -O2 -ip -opt-prefetch -ftz -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g # Profiling flags ################# # [PROFILE] FC : -p -g -FCFLAGS : -xSSE4.2 -O2 -ip -opt-prefetch -ftz +FCFLAGS : -xSSE4.2 -O2 -ip -ftz # Debugging flags ################# diff --git a/install/Downloads/.gitignore b/install/Downloads/.gitignore deleted file mode 100644 index e69de29b..00000000 diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 51fdb52b..761a0130 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -1,6 +1,7 @@ _build ezfio.ml .gitignore +Git.ml Input_auto_generated.ml Input_determinants.ml Input_hartree_fock.ml diff --git a/plugins/CAS_SD/cas_s_selected.irp.f b/plugins/CAS_SD/cas_s_selected.irp.f index 7a72a243..b1fd542a 100644 --- a/plugins/CAS_SD/cas_s_selected.irp.f +++ b/plugins/CAS_SD/cas_s_selected.irp.f @@ -30,7 +30,7 @@ program full_ci endif do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) - call H_apply_CAS_S_selected(pt2, norm_pert, H_pert_diag, N_st) + call H_apply_CAS_S_selected_no_skip(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef PROVIDE psi_det diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index 7fd5d248..325d7860 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -637,7 +637,7 @@ def ninja_binaries_rule(): # c m d # # ~#~#~ # - l_cmd = ["cd $module_abs/IRPF90_temp", "ninja $out && [[ -x $out ]] && touch $out"] + l_cmd = ["cd $module_abs/IRPF90_temp", "ninja $out && for i in $out ; do [ -x $$i ] && touch $$i ; done"] # ~#~#~#~#~#~ # # s t r i n g # diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 9aeb658e..e5d243f4 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -10,7 +10,7 @@ double precision :: ck, cl, ckl double precision :: phase integer :: h1,h2,p1,p2,s1,s2, degree - integer :: exc(0:2,2,2),n_occ_alpha + integer :: exc(0:2,2,2),n_occ(2) double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) if(only_single_double_dm)then @@ -22,7 +22,7 @@ one_body_dm_mo_beta = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & - !$OMP tmp_a, tmp_b, n_occ_alpha)& + !$OMP tmp_a, tmp_b, n_occ)& !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,& !$OMP mo_tot_num) @@ -31,8 +31,7 @@ tmp_b = 0.d0 !$OMP DO SCHEDULE(dynamic) do k=1,N_det - call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int) - call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int) + call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int) do m=1,N_states ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) do l=1,elec_alpha_num @@ -182,13 +181,10 @@ subroutine set_natural_mos END_DOC character*(64) :: label double precision, allocatable :: tmp(:,:) - allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2))) - ! Negation to have the occupied MOs first after the diagonalization - tmp = one_body_dm_mo label = "Natural" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,-1) - deallocate(tmp) +! call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,label,-1) + call mo_as_svd_vectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,label) end subroutine save_natural_mos diff --git a/src/MO_Basis/mo_overlap.irp.f b/src/MO_Basis/mo_overlap.irp.f index 9b8b48f0..c7e146bc 100644 --- a/src/MO_Basis/mo_overlap.irp.f +++ b/src/MO_Basis/mo_overlap.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, mo_overlap,(mo_tot_num_align,mo_tot_num)] integer :: i,j,n,l double precision :: f integer :: lmax - lmax = iand(ao_num,-4) + lmax = (ao_num/4) * 4 !$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) & !$OMP PRIVATE(i,j,n,l) & !$OMP SHARED(mo_overlap,mo_coef,ao_overlap, & diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 0d8ef5fa..462addc0 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -100,6 +100,53 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) mo_label = label end +subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label) + implicit none + integer,intent(in) :: lda,m,n + character*(64), intent(in) :: label + double precision, intent(in) :: matrix(lda,n) + + integer :: i,j + double precision, allocatable :: mo_coef_new(:,:), U(:,:),D(:), A(:,:), Vt(:,:), work(:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, U, Vt, A + + call write_time(output_mo_basis) + if (m /= mo_tot_num) then + print *, irp_here, ': Error : m/= mo_tot_num' + stop 1 + endif + + allocate(A(lda,n),U(lda,n),mo_coef_new(ao_num_align,m),D(m),Vt(lda,n)) + + do j=1,n + do i=1,m + A(i,j) = matrix(i,j) + enddo + enddo + mo_coef_new = mo_coef + + call svd(A,lda,U,lda,D,Vt,lda,m,n) + + write (output_mo_basis,'(A)'), 'MOs are now **'//trim(label)//'**' + write (output_mo_basis,'(A)'), '' + write (output_mo_basis,'(A)'), 'Eigenvalues' + write (output_mo_basis,'(A)'), '-----------' + write (output_mo_basis,'(A)'), '' + write (output_mo_basis,'(A)'), '======== ================' + + do i=1,m + write (output_mo_basis,'(I8,X,F16.10)'), i,D(i) + enddo + write (output_mo_basis,'(A)'), '======== ================' + write (output_mo_basis,'(A)'), '' + + call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1)) + deallocate(A,mo_coef_new,U,Vt,D) + call write_time(output_mo_basis) + + mo_label = label +end + subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label) implicit none integer,intent(in) :: n,m diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index e3ef0bfe..82e8a4c1 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -70,9 +70,8 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) double precision :: Vt(lda,n) double precision :: D(n) double precision :: S_half(lda,n) - double precision,allocatable :: work(:) - !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D, work - integer :: info, lwork, i, j, k + !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D + integer :: info, i, j, k call svd(overlap,lda,U,ldc,D,Vt,lda,m,n) @@ -82,7 +81,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP DO do i=1,n - if ( D(i) < 1.d-12 ) then + if ( D(i) < 1.d-6 ) then D(i) = 0.d0 else D(i) = 1.d0/dsqrt(D(i))