From ba3caadcc2631b1d07c66a5f3a93eb4be91cfd46 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Jun 2014 16:20:07 +0200 Subject: [PATCH 1/9] Added Cholesky MO routine --- src/MOs/cholesky_mo.irp.f | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 src/MOs/cholesky_mo.irp.f diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f new file mode 100644 index 00000000..d587d364 --- /dev/null +++ b/src/MOs/cholesky_mo.irp.f @@ -0,0 +1,39 @@ +subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank) + implicit none + integer, intent(in) :: n,m, LDC + double precision, intent(in) :: P(LDC,n) + double precision, intent(out) :: C(LDC,m) + double precision, intent(in) :: tol_in + integer, intent(out) :: rank + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:) + !DEC$ ATTRIBUTES ALIGN: 32 :: W + !DEC$ ATTRIBUTES ALIGN: 32 :: work + !DEC$ ATTRIBUTES ALIGN: 32 :: ipiv + + allocate(W(LDC,n),work(2*n)) + tol=tol_in + + info = 0 + do i=1,n + do k=1,i + W(i,k) = P(i,k) + enddo + do k=i+1,n + W(i,k) = 0. + enddo + enddo + call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info ) + do i=1,n + do k=1,min(m,rank) + C(ipiv(i),k) = W(i,k) + enddo + enddo + + deallocate(W,work) +end + From 40df4452cf15660c680351838cbb73848e8bbe12 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Jun 2014 16:21:20 +0200 Subject: [PATCH 2/9] Added missing NEEDED_MODULES --- src/MOGuess/NEEDED_MODULES | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 src/MOGuess/NEEDED_MODULES diff --git a/src/MOGuess/NEEDED_MODULES b/src/MOGuess/NEEDED_MODULES new file mode 100644 index 00000000..27c50e45 --- /dev/null +++ b/src/MOGuess/NEEDED_MODULES @@ -0,0 +1,2 @@ +AOs Ezfio_files MOs Nuclei Output Utils MonoInts + From 3beea8d2302ce4666d190116b71988976931e2af Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Jun 2014 16:34:47 +0200 Subject: [PATCH 3/9] Added mo_occ and mo_density_matrix providers --- src/MOs/NEEDED_MODULES | 2 +- src/MOs/README.rst | 15 +++++++++++-- src/MOs/cholesky_mo.irp.f | 44 +++++++++++++++++++++++++++++++++++++++ src/MOs/mos.ezfio_config | 1 + src/MOs/mos.irp.f | 26 +++++++++++++++++++++-- src/MOs/utils.irp.f | 1 + 6 files changed, 84 insertions(+), 5 deletions(-) diff --git a/src/MOs/NEEDED_MODULES b/src/MOs/NEEDED_MODULES index 6ec1892c..9bbf1b74 100644 --- a/src/MOs/NEEDED_MODULES +++ b/src/MOs/NEEDED_MODULES @@ -1 +1 @@ -AOs Ezfio_files Nuclei Output Utils +AOs Ezfio_files Nuclei Output Utils Electrons diff --git a/src/MOs/README.rst b/src/MOs/README.rst index 659db5d5..f27425bc 100644 --- a/src/MOs/README.rst +++ b/src/MOs/README.rst @@ -27,6 +27,7 @@ Needed Modules * `Nuclei `_ * `Output `_ * `Utils `_ +* `Electrons `_ Documentation ============= @@ -34,12 +35,19 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`cholesky_mo `_ + Cholesky decomposition of MO Density matrix to + generate MOs + +`mo_density_matrix `_ + Density matrix in MO basis + `mo_coef `_ Molecular orbital coefficients on AO basis set mo_coef(i,j) = coefficient of the ith ao on the jth mo mo_label : Label characterizing the MOS (local, canonical, natural, etc) -`mo_coef_transp `_ +`mo_coef_transp `_ Molecular orbital coefficients on AO basis set `mo_label `_ @@ -47,13 +55,16 @@ Documentation mo_coef(i,j) = coefficient of the ith ao on the jth mo mo_label : Label characterizing the MOS (local, canonical, natural, etc) +`mo_occ `_ + MO occupation numbers + `mo_tot_num `_ Total number of molecular orbitals and the size of the keys corresponding `mo_tot_num_align `_ Aligned variable for dimensioning of arrays -`mo_as_eigvectors_of_mo_matrix `_ +`mo_as_eigvectors_of_mo_matrix `_ Undocumented `save_mos `_ diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f index d587d364..d3e04dfb 100644 --- a/src/MOs/cholesky_mo.irp.f +++ b/src/MOs/cholesky_mo.irp.f @@ -1,5 +1,9 @@ subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank) implicit none + BEGIN_DOC +! Cholesky decomposition of MO Density matrix to +! generate MOs + END_DOC integer, intent(in) :: n,m, LDC double precision, intent(in) :: P(LDC,n) double precision, intent(out) :: C(LDC,m) @@ -37,3 +41,43 @@ subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank) deallocate(W,work) end +BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis + END_DOC + integer :: i,j,k + mo_density_matrix = 0.d0 + do k=1,mo_tot_num + if (mo_occ(k) == 0) then + cycle + endif + do j=1,ao_num + do i=1,ao_num + mo_density_matrix(i,j) = mo_density_matrix(i,j) + & + mo_occ(k) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis (virtual MOs) + END_DOC + integer :: i,j,k + mo_density_matrix = 0.d0 + do k=1,mo_tot_num + if (mo_occ(k) == 0) then + cycle + endif + do j=1,ao_num + do i=1,ao_num + mo_density_matrix(i,j) = mo_density_matrix(i,j) + & + (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/MOs/mos.ezfio_config b/src/MOs/mos.ezfio_config index a0eda491..b43f23b9 100644 --- a/src/MOs/mos.ezfio_config +++ b/src/MOs/mos.ezfio_config @@ -2,4 +2,5 @@ mo_basis mo_tot_num integer mo_coef double precision (ao_basis_ao_num,mo_basis_mo_tot_num) mo_label character*(64) + mo_occ double precision (mo_basis_mo_tot_num) diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f index b081d8ce..738d2c2a 100644 --- a/src/MOs/mos.irp.f +++ b/src/MOs/mos.irp.f @@ -32,7 +32,7 @@ END_PROVIDER !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer logical :: exists PROVIDE ezfio_filename - + !Label call ezfio_has_mo_basis_mo_label(exists) if (exists) then @@ -40,7 +40,7 @@ END_PROVIDER else mo_label = 'no_label' endif - + ! Coefs allocate(buffer(ao_num,mo_tot_num)) buffer = 0.d0 @@ -75,3 +75,25 @@ BEGIN_PROVIDER [ double precision, mo_coef_transp, (mo_tot_num_align,ao_num) ] END_PROVIDER +BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ] + implicit none + BEGIN_DOC + ! MO occupation numbers + END_DOC + PROVIDE ezfio_filename + logical :: exists + call ezfio_has_mo_basis_mo_occ(exists) + if (exists) then + call ezfio_get_mo_basis_mo_occ(mo_occ) + else + mo_occ = 0.d0 + integer :: i + do i=1,elec_beta_num + mo_occ(i) = 2.d0 + enddo + do i=elec_beta_num+1,elec_alpha_num + mo_occ(i) = 1.d0 + enddo + endif +END_PROVIDER + diff --git a/src/MOs/utils.irp.f b/src/MOs/utils.irp.f index 186f5efc..ee04ac4b 100644 --- a/src/MOs/utils.irp.f +++ b/src/MOs/utils.irp.f @@ -14,6 +14,7 @@ subroutine save_mos enddo enddo call ezfio_set_mo_basis_mo_coef(buffer) + call ezfio_set_mo_basis_mo_occ(mo_occ) deallocate (buffer) end From 42d8b4c4043d9b8b013303f403472ffc215bdf99 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Jun 2014 22:08:53 +0200 Subject: [PATCH 4/9] Improved Hartree-Fock --- src/AOs/aos.irp.f | 40 ---------------- src/Hartree_Fock/HF_density_matrix_ao.irp.f | 32 ++++++------- src/Hartree_Fock/mo_SCF_iterations.irp.f | 53 +++++++++++++++++---- src/MOs/cholesky_mo.irp.f | 15 +++--- 4 files changed, 66 insertions(+), 74 deletions(-) diff --git a/src/AOs/aos.irp.f b/src/AOs/aos.irp.f index 4b3a74c0..396d9aab 100644 --- a/src/AOs/aos.irp.f +++ b/src/AOs/aos.irp.f @@ -93,46 +93,6 @@ END_PROVIDER enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_overlap, (ao_num_align,ao_num) ] - implicit none - BEGIN_DOC -! matrix of the overlap for tha AOs -! S(i,j) = overlap between the ith and the jth atomic basis function - END_DOC - integer :: i,j,k,l,nz,num_i,num_j,powA(3),powB(3) - double precision :: accu,overlap_x,overlap_y,overlap_z,A_center(3),B_center(3),norm - nz=100 - do i = 1, ao_num - num_i = ao_nucl(i) - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) - A_center(1)=nucl_coord(num_i,1) - A_center(2)=nucl_coord(num_i,2) - A_center(3)=nucl_coord(num_i,3) - do j = i, ao_num - num_j = ao_nucl(j) - powB(1) = ao_power(j,1) - powB(2) = ao_power(j,2) - powB(3) = ao_power(j,3) - B_center(1)=nucl_coord(num_j,1) - B_center(2)=nucl_coord(num_j,2) - B_center(3)=nucl_coord(num_j,3) - accu = 0.d0 - do k = 1, ao_prim_num(i) - do l = 1, ao_prim_num(j) - call overlap_gaussian_xyz(A_center,B_center,ao_expo(i,k),ao_expo(j,l),powA,powB,overlap_x,overlap_y,overlap_z,norm,nz) - accu = accu + norm * ao_coef(i,k) * ao_coef(j,l) - enddo - enddo - ao_overlap(i,j) = accu - ao_overlap(j,i) = accu - enddo - enddo - -END_PROVIDER - - BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ] &BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ] implicit none diff --git a/src/Hartree_Fock/HF_density_matrix_ao.irp.f b/src/Hartree_Fock/HF_density_matrix_ao.irp.f index b0d58344..3183025b 100644 --- a/src/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/src/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -5,13 +5,13 @@ ! Alpha and Beta density matrix in the AO basis END_DOC integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) + integer, allocatable :: occ(:,:) - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) + allocate ( occ(elec_alpha_num,2) ) + call bitstring_to_list( HF_bitmask(1,1), occ(1,1), j, N_int) ASSERT ( j==elec_alpha_num ) - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) + call bitstring_to_list( HF_bitmask(1,2), occ(1,2), j, N_int) ASSERT ( j==elec_beta_num ) do j=1,ao_num @@ -21,8 +21,8 @@ HF_density_matrix_ao_beta (i,j) = 0.d0 enddo do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) + l1 = occ(k,1) + l2 = occ(k,2) !DIR$ VECTOR ALIGNED do i=1,ao_num HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& @@ -32,7 +32,7 @@ enddo enddo do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) + l1 = occ(k,1) !DIR$ VECTOR ALIGNED do i=1,ao_num HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& @@ -40,7 +40,7 @@ enddo enddo enddo - deallocate(mo_occ) + deallocate(occ) END_PROVIDER BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] @@ -49,13 +49,13 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] ! Density matrix in the AO basis END_DOC integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) + integer, allocatable :: occ(:,:) - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) + allocate ( occ(elec_alpha_num,2) ) + call bitstring_to_list( HF_bitmask(1,1), occ(1,1), j, N_int) ASSERT ( j==elec_alpha_num ) - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) + call bitstring_to_list( HF_bitmask(1,2), occ(1,2), j, N_int) ASSERT ( j==elec_beta_num ) do j=1,ao_num @@ -64,8 +64,8 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] HF_density_matrix_ao(i,j) = 0.d0 enddo do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) + l1 = occ(k,1) + l2 = occ(k,2) !DIR$ VECTOR ALIGNED do i=1,ao_num HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & @@ -74,7 +74,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] enddo enddo do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) + l1 = occ(k,1) !DIR$ VECTOR ALIGNED do i=1,ao_num HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & @@ -82,6 +82,6 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] enddo enddo enddo - deallocate(mo_occ) + deallocate(occ) END_PROVIDER diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 01cab413..cf38dbd1 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -3,21 +3,55 @@ program scf_iteration implicit none double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral double precision :: E0 - integer :: i_it - + integer :: i_it, i, j, k + integer, allocatable :: iorder(:) + double precision, allocatable :: DM_occ(:,:), E_new(:), R(:,:) + E0 = HF_energy - i_it = 0 - n_it_scf_max = 10 - SCF_energy_before = huge(1.d0) + i_it = 1 + n_it_scf_max = 100 + SCF_energy_before = 0.d0 SCF_energy_after = E0 print *, E0 mo_label = "Canonical" thresh_SCF = 1.d-10 - do while (i_it < 40 .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF) - SCF_energy_before = SCF_energy_after - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef mo_label + DM_occ = mo_density_matrix + allocate (DM_occ(size(mo_density_matrix,1),mo_tot_num), & + E_new(mo_tot_num), R(mo_tot_num,mo_tot_num), iorder(mo_tot_num)) + do while (i_it < n_it_scf_max .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF) + if (SCF_energy_after <= SCF_energy_before+1.d-4) then + mo_coef = eigenvectors_Fock_matrix_mo + TOUCH mo_coef mo_label + DM_occ = mo_density_matrix + else + DM_occ = mo_density_matrix + mo_coef = eigenvectors_Fock_matrix_mo + TOUCH mo_coef mo_label mo_integrals_map + DM_occ = DM_occ + 0.0d0*mo_density_matrix + integer :: rank + call cholesky_mo(ao_num,mo_tot_num,DM_occ,size(DM_occ,1),mo_coef,size(mo_coef,1),-1.d0,rank) + print *, rank + TOUCH mo_coef mo_label + call orthonormalize_mos + call find_rotation(eigenvectors_Fock_matrix_mo,mo_tot_num_align,mo_coef,ao_num,R, mo_tot_num) + do i=1,mo_tot_num + iorder(i) = i + E_new(i) = 0.d0 + do k=1,mo_tot_num + E_new(i) += R(k,i)*R(k,i)*diagonal_fock_matrix_mo(k) + enddo + enddo + call dsort(E_new(1),iorder(1),mo_tot_num) + eigenvectors_Fock_matrix_mo = mo_coef + do j=1,mo_tot_num + do i=1,ao_num + mo_coef(i,j) = eigenvectors_Fock_matrix_mo(i,iorder(j)) + enddo + enddo + TOUCH mo_coef mo_label mo_integrals_map + endif call clear_mo_map + SCF_energy_before = SCF_energy_after SCF_energy_after = HF_energy print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after) i_it +=1 @@ -31,6 +65,7 @@ program scf_iteration stop 'Failed' endif mo_label = "Canonical" + deallocate (DM_occ) TOUCH mo_label mo_coef call save_mos diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f index d3e04dfb..9afea3e6 100644 --- a/src/MOs/cholesky_mo.irp.f +++ b/src/MOs/cholesky_mo.irp.f @@ -1,11 +1,11 @@ -subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank) +subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank) implicit none BEGIN_DOC ! Cholesky decomposition of MO Density matrix to ! generate MOs END_DOC - integer, intent(in) :: n,m, LDC - double precision, intent(in) :: P(LDC,n) + integer, intent(in) :: n,m, LDC, LDP + double precision, intent(in) :: P(LDP,n) double precision, intent(out) :: C(LDC,m) double precision, intent(in) :: tol_in integer, intent(out) :: rank @@ -49,7 +49,7 @@ BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_ integer :: i,j,k mo_density_matrix = 0.d0 do k=1,mo_tot_num - if (mo_occ(k) == 0) then + if (mo_occ(k) == 0.d0) then cycle endif do j=1,ao_num @@ -67,14 +67,11 @@ BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, ! Density matrix in MO basis (virtual MOs) END_DOC integer :: i,j,k - mo_density_matrix = 0.d0 + mo_density_matrix_virtual = 0.d0 do k=1,mo_tot_num - if (mo_occ(k) == 0) then - cycle - endif do j=1,ao_num do i=1,ao_num - mo_density_matrix(i,j) = mo_density_matrix(i,j) + & + mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + & (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k) enddo enddo From 270dc053facd764aace2b68993526c1d4c9120d2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Jun 2014 23:02:41 +0200 Subject: [PATCH 5/9] Bug in generators --- src/Generators_full/generators.irp.f | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index 3f3fba86..9d8761d2 100644 --- a/src/Generators_full/generators.irp.f +++ b/src/Generators_full/generators.irp.f @@ -47,7 +47,13 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] ! For Single reference wave functions, the generator is the ! Hartree-Fock determinant END_DOC - psi_generators = psi_det_sorted + integer :: i, k + do i=1,N_det + do k=1,N_int + psi_generators(k,1,i) = psi_det_sorted(k,1,i) + psi_generators(k,2,i) = psi_det_sorted(k,2,i) + enddo + enddo END_PROVIDER From 89a7e3a6447bde2ef935e903b23ec481aa34452b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 19 Jun 2014 17:58:45 +0200 Subject: [PATCH 6/9] DIIS on the way... --- src/Dets/README.rst | 142 ++++++++++---------- src/Hartree_Fock/Fock_matrix.irp.f | 43 ++++++ src/Hartree_Fock/HF_density_matrix_ao.irp.f | 98 ++++---------- src/Hartree_Fock/README.rst | 49 +++---- src/Hartree_Fock/SCF.irp.f | 95 +++++++++++++ src/Hartree_Fock/diagonalize_fock.irp.f | 54 ++++++-- src/Hartree_Fock/mo_SCF_iterations.irp.f | 4 +- src/MOGuess/README.rst | 8 +- src/MOs/cholesky_mo.irp.f | 80 +++++++++++ src/MOs/mos.ezfio_config | 2 + src/MOs/mos.irp.f | 22 +++ src/Perturbation/README.rst | 42 ++++-- 12 files changed, 450 insertions(+), 189 deletions(-) create mode 100644 src/Hartree_Fock/SCF.irp.f create mode 100644 src/MOs/cholesky_mo.irp.f diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 35c61523..09352c44 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -50,24 +50,24 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`copy_h_apply_buffer_to_wf `_ +`copy_h_apply_buffer_to_wf `_ Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det after calling this function. -`fill_h_apply_buffer_no_selection `_ +`fill_h_apply_buffer_no_selection `_ Fill the H_apply buffer with determiants for CISD -`h_apply_buffer_allocated `_ +`h_apply_buffer_allocated `_ Buffer of determinants/coefficients/perturbative energy for H_apply. Uninitialized. Filled by H_apply subroutines. -`h_apply_threshold `_ +`h_apply_threshold `_ Theshold on | | -`resize_h_apply_buffer `_ +`resize_h_apply_buffer `_ Undocumented -`cisd_sc2 `_ +`cisd_sc2 `_ CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) .br dets_in : bitmasks corresponding to determinants @@ -83,29 +83,29 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`repeat_excitation `_ +`repeat_excitation `_ Undocumented -`connected_to_ref `_ +`connected_to_ref `_ Undocumented -`det_is_not_or_may_be_in_ref `_ +`det_is_not_or_may_be_in_ref `_ If true, det is not in ref If false, det may be in ref -`is_in_wavefunction `_ +`is_in_wavefunction `_ Undocumented -`key_pattern_not_in_ref `_ +`key_pattern_not_in_ref `_ Min and max values of the integers of the keys of the reference -`davidson_converged `_ +`davidson_converged `_ True if the Davidson algorithm is converged -`davidson_criterion `_ +`davidson_criterion `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`davidson_diag `_ +`davidson_diag `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -123,7 +123,7 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_diag_hjj `_ +`davidson_diag_hjj `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -143,114 +143,114 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_iter_max `_ +`davidson_iter_max `_ Max number of Davidson iterations -`davidson_sze_max `_ +`davidson_sze_max `_ Max number of Davidson sizes -`davidson_threshold `_ +`davidson_threshold `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`det_search_key `_ +`det_search_key `_ Return an integer*8 corresponding to a determinant index for searching -`n_det `_ +`n_det `_ Number of determinants in the wave function -`n_det_max_jacobi `_ +`n_det_max_jacobi `_ Maximum number of determinants diagonalized my jacobi -`n_states `_ +`n_states `_ Number of states to consider -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`psi_det `_ +`psi_det `_ The wave function determinants. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`double_exc_bitmask `_ +`double_exc_bitmask `_ double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1 double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1 double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2 double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2 for a given couple of hole/particle excitations i. -`n_double_exc_bitmasks `_ +`n_double_exc_bitmasks `_ Number of double excitation bitmasks -`n_single_exc_bitmasks `_ +`n_single_exc_bitmasks `_ Number of single excitation bitmasks -`single_exc_bitmask `_ +`single_exc_bitmask `_ single_exc_bitmask(:,1,i) is the bitmask for holes single_exc_bitmask(:,2,i) is the bitmask for particles for a given couple of hole/particle excitations i. -`ci_eigenvectors `_ +`ci_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_electronic_energy `_ +`ci_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_energy `_ +`ci_energy `_ N_states lowest eigenvalues of the CI matrix -`diag_algorithm `_ +`diag_algorithm `_ Diagonalization algorithm (Davidson or Lapack) -`diagonalize_ci `_ +`diagonalize_ci `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`ci_sc2_eigenvectors `_ +`ci_sc2_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_sc2_electronic_energy `_ +`ci_sc2_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_sc2_energy `_ +`ci_sc2_energy `_ N_states lowest eigenvalues of the CI matrix -`diagonalize_ci_sc2 `_ +`diagonalize_ci_sc2 `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`filter_connected `_ +`filter_connected `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -261,7 +261,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_davidson `_ +`filter_connected_davidson `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -272,7 +272,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0 `_ +`filter_connected_i_h_psi0 `_ returns the array idx which contains the index of the .br determinants in the array key1 that interact @@ -281,7 +281,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0_sc2 `_ +`filter_connected_i_h_psi0_sc2 `_ standard filter_connected_i_H_psi but returns in addition .br the array of the index of the non connected determinants to key1 @@ -292,69 +292,69 @@ Documentation .br to repeat the excitations -`get_s2 `_ +`get_s2 `_ Returns -`get_s2_u0 `_ +`get_s2_u0 `_ Undocumented -`s_z `_ +`s_z `_ Undocumented -`s_z2_sz `_ +`s_z2_sz `_ Undocumented -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem -`decode_exc `_ +`decode_exc `_ Decodes the exc arrays returned by get_excitation. h1,h2 : Holes p1,p2 : Particles s1,s2 : Spins (1:alpha, 2:beta) degree : Degree of excitation -`det_connections `_ +`det_connections `_ .br -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes -`get_double_excitation `_ +`get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase -`get_excitation `_ +`get_excitation `_ Returns the excitation operators between two determinants and the phase -`get_excitation_degree `_ +`get_excitation_degree `_ Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants -`get_mono_excitation `_ +`get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants .br H_jj : array of -`i_h_j `_ +`i_h_j `_ Returns where i and j are determinants -`i_h_psi `_ +`i_h_psi `_ for the various Nstates -`i_h_psi_sc2 `_ +`i_h_psi_sc2 `_ for the various Nstate .br returns in addition @@ -367,10 +367,10 @@ Documentation .br to repeat the excitations -`n_con_int `_ +`n_con_int `_ Number of integers to represent the connections between determinants -`h_matrix_all_dets `_ +`h_matrix_all_dets `_ H matrix on the basis of the slater deter;inants defined by psi_det diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 42bed645..69aa7e26 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -217,3 +217,46 @@ BEGIN_PROVIDER [ double precision, HF_energy ] END_PROVIDER +BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if (elec_alpha_num == elec_beta_num) then + integer :: i,j + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,ao_num_align + Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j) + enddo + enddo + else + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + Fock_matrix_mo, size(Fock_matrix_mo,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + Fock_matrix_ao, size(Fock_matrix_ao,1)) + + deallocate(T) + endif +END_PROVIDER + diff --git a/src/Hartree_Fock/HF_density_matrix_ao.irp.f b/src/Hartree_Fock/HF_density_matrix_ao.irp.f index b0d58344..e85e4a6c 100644 --- a/src/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/src/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -1,46 +1,27 @@ - BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ] -&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ] +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ] implicit none BEGIN_DOC - ! Alpha and Beta density matrix in the AO basis + ! Alpha density matrix in the AO basis END_DOC - integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) - ASSERT ( j==elec_alpha_num ) + call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ] + implicit none + BEGIN_DOC + ! Beta density matrix in the AO basis + END_DOC - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) - ASSERT ( j==elec_beta_num ) - - do j=1,ao_num - !DIR$ VECTOR ALIGNED - do i=1,ao_num_align - HF_density_matrix_ao_alpha(i,j) = 0.d0 - HF_density_matrix_ao_beta (i,j) = 0.d0 - enddo - do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& - mo_coef(i,l1) * mo_coef(j,l1) - HF_density_matrix_ao_beta (i,j) = HF_density_matrix_ao_beta (i,j) +& - mo_coef(i,l2) * mo_coef(j,l2) - enddo - enddo - do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& - mo_coef(i,l1) * mo_coef(j,l1) - enddo - enddo - enddo - deallocate(mo_occ) + call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1)) + END_PROVIDER BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] @@ -48,40 +29,13 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] BEGIN_DOC ! Density matrix in the AO basis END_DOC - integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1)) + if (elec_alpha_num== elec_beta_num) then + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_alpha + else + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_beta ,1)) + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_beta + endif - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) - ASSERT ( j==elec_alpha_num ) - - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) - ASSERT ( j==elec_beta_num ) - - do j=1,ao_num - !DIR$ VECTOR ALIGNED - do i=1,ao_num_align - HF_density_matrix_ao(i,j) = 0.d0 - enddo - do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & - mo_coef(i,l1) * mo_coef(j,l1) + & - mo_coef(i,l2) * mo_coef(j,l2) - enddo - enddo - do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & - mo_coef(i,l1) * mo_coef(j,l1) - enddo - enddo - enddo - deallocate(mo_occ) END_PROVIDER diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 074928ad..7173d418 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -26,19 +26,22 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fock_matrix_alpha_ao `_ +`fock_matrix_alpha_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_alpha_mo `_ +`fock_matrix_alpha_mo `_ Fock matrix on the MO basis -`fock_matrix_beta_ao `_ +`fock_matrix_ao `_ + Fock matrix in AO basis set + +`fock_matrix_beta_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_beta_mo `_ +`fock_matrix_beta_mo `_ Fock matrix on the MO basis -`fock_matrix_diag_mo `_ +`fock_matrix_diag_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -53,7 +56,7 @@ Documentation K = Fb - Fa .br -`fock_matrix_mo `_ +`fock_matrix_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -68,49 +71,49 @@ Documentation K = Fb - Fa .br -`hf_energy `_ +`hf_energy `_ Hartree-Fock energy -`hf_density_matrix_ao `_ +`hf_density_matrix_ao `_ Density matrix in the AO basis -`hf_density_matrix_ao_alpha `_ - Alpha and Beta density matrix in the AO basis +`hf_density_matrix_ao_alpha `_ + Alpha density matrix in the AO basis -`hf_density_matrix_ao_beta `_ - Alpha and Beta density matrix in the AO basis +`hf_density_matrix_ao_beta `_ + Beta density matrix in the AO basis -`diagonal_fock_matrix_mo `_ +`diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`eigenvectors_fock_matrix_mo `_ +`eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`scf_iteration `_ +`xcf_iteration `_ Undocumented -`do_diis `_ +`do_diis `_ If True, compute integrals on the fly -`n_it_scf_max `_ +`n_it_scf_max `_ Maximum number of SCF iterations -`thresh_scf `_ +`thresh_scf `_ Threshold on the convergence of the Hartree Fock energy -`bi_elec_ref_bitmask_energy `_ +`bi_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`kinetic_ref_bitmask_energy `_ +`kinetic_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`mono_elec_ref_bitmask_energy `_ +`mono_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`nucl_elec_ref_bitmask_energy `_ +`nucl_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`ref_bitmask_energy `_ +`ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f new file mode 100644 index 00000000..cb6e8cbf --- /dev/null +++ b/src/Hartree_Fock/SCF.irp.f @@ -0,0 +1,95 @@ +BEGIN_PROVIDER [ integer, it_scf ] + implicit none + BEGIN_DOC + ! Number of the current SCF iteration + END_DOC + it_scf = 0 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,0:n_it_scf_max) ] + implicit none + BEGIN_DOC + ! Density matrices at every SCF iteration + END_DOC + SCF_density_matrices = 0.d0 +END_PROVIDER + +subroutine insert_new_SCF_density_matrix + implicit none + integer :: i,j + do j=1,ao_num + do i=1,ao_num + SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,1,0) += HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) + SCF_density_matrices(i,j,2,0) += HF_density_matrix_ao_beta(i,j) + enddo + enddo +end + +subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) + implicit none + integer, intent(in) :: LDFMO ! size(FMO,1) + integer, intent(in) :: LDFAO ! size(FAO,1) + double precision, intent(in) :: FMO(LDFMO,*) + double precision, intent(out) :: FAO(LDFAO,*) + + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + FMO, size(FMO,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + FAO, size(FAO,1)) + deallocate(T,M) +end + +subroutine DIIS_step + implicit none + integer :: i,j + double precision :: c + c = 1.d0/dble(it_scf) + do j=1,ao_num + do i=1,ao_num + HF_density_matrix_ao_alpha(i,j) = SCF_density_matrices(i,j,1,0) * c + HF_density_matrix_ao_beta (i,j) = SCF_density_matrices(i,j,2,0) * c + enddo + enddo + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + +! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1), & +! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) +! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1), & +! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) +! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta +end + +subroutine scf_iteration + implicit none + integer :: i,j + do i=1,n_it_scf_max + it_scf += 1 + SOFT_TOUCH it_scf + mo_coef = eigenvectors_Fock_matrix_mo + TOUCH mo_coef + call insert_new_SCF_density_matrix + call DIIS_step + print *, HF_energy + enddo +end diff --git a/src/Hartree_Fock/diagonalize_fock.irp.f b/src/Hartree_Fock/diagonalize_fock.irp.f index c34a1415..c5faea9c 100644 --- a/src/Hartree_Fock/diagonalize_fock.irp.f +++ b/src/Hartree_Fock/diagonalize_fock.irp.f @@ -5,16 +5,52 @@ ! Diagonal Fock matrix in the MO basis END_DOC - double precision, allocatable :: R(:,:) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: R + integer :: i,j + integer :: liwork, lwork, n, info + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:), F(:,:), S(:,:) - allocate(R(mo_tot_num_align,mo_tot_num)) + allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) ) + do j=1,ao_num + do i=1,ao_num + S(i,j) = ao_overlap(i,j) + F(i,j) = Fock_matrix_ao(i,j) + enddo + enddo + + n = ao_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n - call lapack_diag(diagonal_Fock_matrix_mo,R,Fock_matrix_mo,size(Fock_matrix_mo,1),& - mo_tot_num) - - call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num,1.d0,mo_coef,size(mo_coef,1),& - R,size(R,1),0.d0,eigenvectors_Fock_matrix_mo,size(eigenvectors_Fock_matrix_mo,1)) - deallocate(R) + allocate(work(lwork), iwork(liwork) ) + + lwork = -1 + liwork = -1 + call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed' + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(work,iwork) + allocate(work(lwork), iwork(liwork) ) + + call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed' + stop 1 + endif + do j=1,ao_num + do i=1,ao_num + eigenvectors_Fock_matrix_mo(i,j) = F(i,j) + enddo + enddo + + deallocate(work, iwork, F, S) END_PROVIDER diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 01cab413..4854861f 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -1,4 +1,4 @@ -program scf_iteration +program xcf_iteration use bitmasks implicit none double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral @@ -32,6 +32,6 @@ program scf_iteration endif mo_label = "Canonical" TOUCH mo_label mo_coef - call save_mos +! call save_mos end diff --git a/src/MOGuess/README.rst b/src/MOGuess/README.rst index d8e72641..90785358 100644 --- a/src/MOGuess/README.rst +++ b/src/MOGuess/README.rst @@ -22,19 +22,19 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`h_core_guess `_ +`h_core_guess `_ Undocumented -`ao_ortho_lowdin_coef `_ +`ao_ortho_lowdin_coef `_ matrix of the coefficients of the mos generated by the orthonormalization by the S^{-1/2} canonical transformation of the aos ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital -`ao_ortho_lowdin_overlap `_ +`ao_ortho_lowdin_overlap `_ overlap matrix of the ao_ortho_lowdin supposed to be the Identity -`ao_ortho_lowdin_nucl_elec_integral `_ +`ao_ortho_lowdin_nucl_elec_integral `_ Undocumented diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f new file mode 100644 index 00000000..97b6abd2 --- /dev/null +++ b/src/MOs/cholesky_mo.irp.f @@ -0,0 +1,80 @@ +subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank) + implicit none + BEGIN_DOC +! Cholesky decomposition of AO Density matrix to +! generate MOs + END_DOC + integer, intent(in) :: n,m, LDC, LDP + double precision, intent(in) :: P(LDP,n) + double precision, intent(out) :: C(LDC,m) + double precision, intent(in) :: tol_in + integer, intent(out) :: rank + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:) + !DEC$ ATTRIBUTES ALIGN: 32 :: W + !DEC$ ATTRIBUTES ALIGN: 32 :: work + !DEC$ ATTRIBUTES ALIGN: 32 :: ipiv + + allocate(W(LDC,n),work(2*n)) + tol=tol_in + + info = 0 + do i=1,n + do k=1,i + W(i,k) = P(i,k) + enddo + do k=i+1,n + W(i,k) = 0. + enddo + enddo + call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info ) + do i=1,n + do k=1,min(m,rank) + C(ipiv(i),k) = W(i,k) + enddo + enddo + + deallocate(W,work) +end + +BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis + END_DOC + integer :: i,j,k + mo_density_matrix = 0.d0 + do k=1,mo_tot_num + if (mo_occ(k) == 0.d0) then + cycle + endif + do j=1,ao_num + do i=1,ao_num + mo_density_matrix(i,j) = mo_density_matrix(i,j) + & + mo_occ(k) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis (virtual MOs) + END_DOC + integer :: i,j,k + mo_density_matrix_virtual = 0.d0 + do k=1,mo_tot_num + do j=1,ao_num + do i=1,ao_num + mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + & + (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/MOs/mos.ezfio_config b/src/MOs/mos.ezfio_config index a0eda491..e34d1aab 100644 --- a/src/MOs/mos.ezfio_config +++ b/src/MOs/mos.ezfio_config @@ -1,5 +1,7 @@ mo_basis + ao_num integer mo_tot_num integer mo_coef double precision (ao_basis_ao_num,mo_basis_mo_tot_num) mo_label character*(64) + mo_occ double precision (mo_basis_mo_tot_num) diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f index b081d8ce..7895ac55 100644 --- a/src/MOs/mos.irp.f +++ b/src/MOs/mos.irp.f @@ -75,3 +75,25 @@ BEGIN_PROVIDER [ double precision, mo_coef_transp, (mo_tot_num_align,ao_num) ] END_PROVIDER +BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ] + implicit none + BEGIN_DOC + ! MO occupation numbers + END_DOC + PROVIDE ezfio_filename + logical :: exists + call ezfio_has_mo_basis_mo_occ(exists) + if (exists) then + call ezfio_get_mo_basis_mo_occ(mo_occ) + else + mo_occ = 0.d0 + integer :: i + do i=1,elec_beta_num + mo_occ(i) = 2.d0 + enddo + do i=elec_beta_num+1,elec_alpha_num + mo_occ(i) = 1.d0 + enddo + endif +END_PROVIDER + diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 806ad292..aad1d0c0 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -82,7 +82,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`pt2_moller_plesset `_ +`pt2_moller_plesset `_ compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution .br for the various n_st states. @@ -92,7 +92,7 @@ Documentation e_2_pert(i) = ^2/(difference of orbital energies) .br -`pt2_epstein_nesbet `_ +`pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states. @@ -102,7 +102,7 @@ Documentation e_2_pert(i) = ^2/( E(i) - ) .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various N_st states. @@ -112,20 +112,46 @@ Documentation c_pert(i) = e_2_pert(i)/ .br -`fill_h_apply_buffer_selection `_ +`pt2_epstein_nesbet_sc2_projected `_ + compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various N_st states, + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + .br + that could be repeated to this determinant. + .br + In addition, for the perturbative energetic contribution you have the standard second order + .br + e_2_pert = ^2/(Delta_E) + .br + and also the purely projected contribution + .br + H_pert_diag = c_pert + +`repeat_all_e_corr `_ + Undocumented + +`fill_h_apply_buffer_selection `_ Fill the H_apply buffer with determiants for the selection -`remove_small_contributions `_ +`remove_small_contributions `_ Remove determinants with small contributions. N_states is assumed to be provided. -`selection_criterion `_ +`selection_criterion `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_factor `_ +`selection_criterion_factor `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_min `_ +`selection_criterion_min `_ Threshold to select determinants. Set by selection routines. From 12c47364ca83852d01a1e6c550c57f385d8a17d7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 19 Jun 2014 22:34:56 +0200 Subject: [PATCH 7/9] Better Hartree-Fock --- src/Hartree_Fock/README.rst | 21 +++ src/Hartree_Fock/SCF.irp.f | 170 ++++++++++++----------- src/Hartree_Fock/mo_SCF_iterations.irp.f | 26 +--- 3 files changed, 114 insertions(+), 103 deletions(-) diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 7173d418..a9d70bd4 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -83,6 +83,27 @@ Documentation `hf_density_matrix_ao_beta `_ Beta density matrix in the AO basis +`fock_mo_to_ao `_ + Undocumented + +`insert_new_scf_density_matrix `_ + Undocumented + +`it_scf `_ + Number of the current SCF iteration + +`scf_density_matrices `_ + Density matrices at every SCF iteration + +`scf_energies `_ + Density matrices at every SCF iteration + +`scf_interpolation_step `_ + Undocumented + +`scf_iterations `_ + Undocumented + `diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f index cb6e8cbf..895ca546 100644 --- a/src/Hartree_Fock/SCF.irp.f +++ b/src/Hartree_Fock/SCF.irp.f @@ -1,95 +1,107 @@ BEGIN_PROVIDER [ integer, it_scf ] - implicit none - BEGIN_DOC - ! Number of the current SCF iteration - END_DOC - it_scf = 0 + implicit none + BEGIN_DOC + ! Number of the current SCF iteration + END_DOC + it_scf = 0 END_PROVIDER -BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,0:n_it_scf_max) ] - implicit none - BEGIN_DOC - ! Density matrices at every SCF iteration - END_DOC - SCF_density_matrices = 0.d0 + BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ] +&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ] + implicit none + BEGIN_DOC + ! Density matrices at every SCF iteration + END_DOC + SCF_density_matrices = 0.d0 + SCF_energies = 0.d0 END_PROVIDER subroutine insert_new_SCF_density_matrix - implicit none - integer :: i,j - do j=1,ao_num - do i=1,ao_num - SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) - SCF_density_matrices(i,j,1,0) += HF_density_matrix_ao_alpha(i,j) - SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) - SCF_density_matrices(i,j,2,0) += HF_density_matrix_ao_beta(i,j) + implicit none + integer :: i,j + do j=1,ao_num + do i=1,ao_num + SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) + enddo enddo - enddo + SCF_energies(it_scf) = HF_energy end subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) - implicit none - integer, intent(in) :: LDFMO ! size(FMO,1) - integer, intent(in) :: LDFAO ! size(FAO,1) - double precision, intent(in) :: FMO(LDFMO,*) - double precision, intent(out) :: FAO(LDFAO,*) - - double precision, allocatable :: T(:,:), M(:,:) - ! F_ao = S C F_mo C^t S - allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - ao_overlap, size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & - M, size(M,1), & - FMO, size(FMO,1), & - 0.d0, & - T, size(T,1)) - call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & - T, size(T,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - M, size(M,1), & - ao_overlap, size(ao_overlap,1), & - 0.d0, & - FAO, size(FAO,1)) - deallocate(T,M) + implicit none + integer, intent(in) :: LDFMO ! size(FMO,1) + integer, intent(in) :: LDFAO ! size(FAO,1) + double precision, intent(in) :: FMO(LDFMO,*) + double precision, intent(out) :: FAO(LDFAO,*) + + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + FMO, size(FMO,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + FAO, size(FAO,1)) + deallocate(T,M) end -subroutine DIIS_step - implicit none - integer :: i,j - double precision :: c - c = 1.d0/dble(it_scf) - do j=1,ao_num - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = SCF_density_matrices(i,j,1,0) * c - HF_density_matrix_ao_beta (i,j) = SCF_density_matrices(i,j,2,0) * c +subroutine SCF_interpolation_step + implicit none + integer :: i,j + double precision :: c + + if (it_scf == 1) then + return + endif + call random_number(c) + do j=1,ao_num + do i=1,ao_num + HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf)+SCF_density_matrices(i,j,1,it_scf-1) * (1.d0 - c) + HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf)+SCF_density_matrices(i,j,2,it_scf-1) * (1.d0 - c) + enddo enddo - enddo - TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta - -! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1), & -! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) -! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1), & -! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) -! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + + ! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1),& + ! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) + ! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1),& + ! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) + ! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta end -subroutine scf_iteration - implicit none - integer :: i,j - do i=1,n_it_scf_max - it_scf += 1 - SOFT_TOUCH it_scf - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef - call insert_new_SCF_density_matrix - call DIIS_step - print *, HF_energy - enddo +subroutine scf_iterations + implicit none + integer :: i,j + do i=1,n_it_scf_max + it_scf += 1 + SOFT_TOUCH it_scf + mo_coef = eigenvectors_Fock_matrix_mo + TOUCH mo_coef + call insert_new_SCF_density_matrix + print *, HF_energy + if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then + call SCF_interpolation_step + endif + if (it_scf>1 ) then + if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then + exit + endif + endif + enddo end diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 4854861f..bae0e6d4 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -6,32 +6,10 @@ program xcf_iteration integer :: i_it E0 = HF_energy - i_it = 0 - n_it_scf_max = 10 - SCF_energy_before = huge(1.d0) - SCF_energy_after = E0 - print *, E0 - mo_label = "Canonical" thresh_SCF = 1.d-10 - do while (i_it < 40 .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF) - SCF_energy_before = SCF_energy_after - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef mo_label - call clear_mo_map - SCF_energy_after = HF_energy - print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after) - i_it +=1 - if(i_it > n_it_scf_max)exit - enddo - - if (i_it >= n_it_scf_max) then - stop 'Failed' - endif - if (SCF_energy_after - E0 > thresh_SCF) then - stop 'Failed' - endif + call scf_iterations mo_label = "Canonical" TOUCH mo_label mo_coef -! call save_mos + call save_mos end From 0fb0b9e2ec49cd5247f4e9cec11f8c22d5980ee6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 20 Jun 2014 15:23:04 +0200 Subject: [PATCH 8/9] Put do_mono/diexcitations in generate_h_apply.py --- scripts/generate_h_apply.py | 5 +- src/CID/H_apply.irp.f | 2 +- src/CID/README.rst | 2 +- src/CID/do_mono_double.irp.f | 19 --- src/CID_SC2_selected/H_apply.irp.f | 2 +- src/CID_SC2_selected/README.rst | 2 +- src/CID_SC2_selected/do_mono_double.irp.f | 19 --- src/CID_selected/H_apply.irp.f | 2 +- src/CID_selected/README.rst | 4 +- src/CID_selected/do_mono_double.irp.f | 19 --- src/CISD/do_mono_double.irp.f | 19 --- src/CISD_SC2_selected/do_mono_double.irp.f | 19 --- src/CISD_selected/do_mono_double.irp.f | 19 --- src/Dets/H_apply_template.f | 8 +- src/Dets/README.rst | 145 +++++++++++---------- src/MonoInts/NEEDED_MODULES | 2 +- src/MonoInts/README.rst | 1 + 17 files changed, 91 insertions(+), 198 deletions(-) delete mode 100644 src/CID/do_mono_double.irp.f delete mode 100644 src/CID_SC2_selected/do_mono_double.irp.f delete mode 100644 src/CID_selected/do_mono_double.irp.f delete mode 100644 src/CISD/do_mono_double.irp.f delete mode 100644 src/CISD_SC2_selected/do_mono_double.irp.f delete mode 100644 src/CISD_selected/do_mono_double.irp.f diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index a9d1c966..a57f869f 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -27,7 +27,7 @@ filter_integrals class H_apply(object): - def __init__(self,sub,SingleRef=False): + def __init__(self,sub,SingleRef=False,do_mono_exc=True, do_double_exc=True): s = {} for k in keywords: s[k] = "" @@ -56,6 +56,9 @@ class H_apply(object): s["omp_do"] = "!$OMP DO SCHEDULE (static)" s["omp_enddo"] = "!$OMP ENDDO NOWAIT" + d = { True : '.True.', False : '.False.'} + s["do_mono_excitations"] = d[do_mono_exc] + s["do_double_excitations"] = d[do_double_exc] s["keys_work"] += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)" s["filter_integrals"] = "array_pairs = .True." diff --git a/src/CID/H_apply.irp.f b/src/CID/H_apply.irp.f index 0df1da38..6a1bb536 100644 --- a/src/CID/H_apply.irp.f +++ b/src/CID/H_apply.irp.f @@ -3,7 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import H_apply -H = H_apply("cisd") +H = H_apply("cisd",do_double_exc=True,do_mono_exc=False) print H END_SHELL diff --git a/src/CID/README.rst b/src/CID/README.rst index 003b51fa..5a0c5026 100644 --- a/src/CID/README.rst +++ b/src/CID/README.rst @@ -36,7 +36,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd `_ +`cisd `_ Undocumented diff --git a/src/CID/do_mono_double.irp.f b/src/CID/do_mono_double.irp.f deleted file mode 100644 index f211879e..00000000 --- a/src/CID/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .False. - -END_PROVIDER diff --git a/src/CID_SC2_selected/H_apply.irp.f b/src/CID_SC2_selected/H_apply.irp.f index 79668af7..e284ab43 100644 --- a/src/CID_SC2_selected/H_apply.irp.f +++ b/src/CID_SC2_selected/H_apply.irp.f @@ -3,7 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * from perturbation import perturbations -s = H_apply("PT2",SingleRef=True) +s = H_apply("PT2",SingleRef=True,do_mono_exc=False,do_double_exc=True) s.set_perturbation("epstein_nesbet_sc2_projected") print s END_SHELL diff --git a/src/CID_SC2_selected/README.rst b/src/CID_SC2_selected/README.rst index b6206850..8688fd41 100644 --- a/src/CID_SC2_selected/README.rst +++ b/src/CID_SC2_selected/README.rst @@ -8,7 +8,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd_sc2_selected `_ +`cisd_sc2_selected `_ Undocumented diff --git a/src/CID_SC2_selected/do_mono_double.irp.f b/src/CID_SC2_selected/do_mono_double.irp.f deleted file mode 100644 index f211879e..00000000 --- a/src/CID_SC2_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .False. - -END_PROVIDER diff --git a/src/CID_selected/H_apply.irp.f b/src/CID_selected/H_apply.irp.f index 91dfb9fc..e3afaa9d 100644 --- a/src/CID_selected/H_apply.irp.f +++ b/src/CID_selected/H_apply.irp.f @@ -4,7 +4,7 @@ from generate_h_apply import * from perturbation import perturbations for perturbation in perturbations: - s = H_apply("cisd_selection_"+perturbation) + s = H_apply("cisd_selection_"+perturbation,do_mono_exc=False) s.set_selection_pt2(perturbation) print s END_SHELL diff --git a/src/CID_selected/README.rst b/src/CID_selected/README.rst index 4fed26e5..8c64bc25 100644 --- a/src/CID_selected/README.rst +++ b/src/CID_selected/README.rst @@ -8,10 +8,10 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`h_apply_cisd_selection `_ +`h_apply_cisd_selection `_ Undocumented -`cisd `_ +`cisd `_ Undocumented diff --git a/src/CID_selected/do_mono_double.irp.f b/src/CID_selected/do_mono_double.irp.f deleted file mode 100644 index f211879e..00000000 --- a/src/CID_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .False. - -END_PROVIDER diff --git a/src/CISD/do_mono_double.irp.f b/src/CISD/do_mono_double.irp.f deleted file mode 100644 index 8a07c292..00000000 --- a/src/CISD/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .True. - -END_PROVIDER diff --git a/src/CISD_SC2_selected/do_mono_double.irp.f b/src/CISD_SC2_selected/do_mono_double.irp.f deleted file mode 100644 index 8a07c292..00000000 --- a/src/CISD_SC2_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .True. - -END_PROVIDER diff --git a/src/CISD_selected/do_mono_double.irp.f b/src/CISD_selected/do_mono_double.irp.f deleted file mode 100644 index 8a07c292..00000000 --- a/src/CISD_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .True. - -END_PROVIDER diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 773fbd92..b795223f 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -420,13 +420,13 @@ subroutine $subroutine($params_main) enddo enddo - if(do_double_excitations)then + if($do_double_excitations)then call $subroutine_diexc(psi_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator $params_post) endif - if(do_mono_excitations)then + if($do_mono_excitations)then call $subroutine_monoexc(psi_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator $params_post) @@ -475,13 +475,13 @@ subroutine $subroutine($params_main) not(psi_generators(k,ispin,i_generator)) ) enddo enddo - if(do_double_excitations)then + if($do_double_excitations)then call $subroutine_diexc(psi_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator $params_post) endif - if(do_mono_excitations)then + if($do_mono_excitations)then call $subroutine_monoexc(psi_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator $params_post) diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 09352c44..4a4e8ff2 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -50,24 +50,24 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`copy_h_apply_buffer_to_wf `_ +`copy_h_apply_buffer_to_wf `_ Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det after calling this function. -`fill_h_apply_buffer_no_selection `_ +`fill_h_apply_buffer_no_selection `_ Fill the H_apply buffer with determiants for CISD -`h_apply_buffer_allocated `_ +`h_apply_buffer_allocated `_ Buffer of determinants/coefficients/perturbative energy for H_apply. Uninitialized. Filled by H_apply subroutines. -`h_apply_threshold `_ +`h_apply_threshold `_ Theshold on | | -`resize_h_apply_buffer `_ +`resize_h_apply_buffer `_ Undocumented -`cisd_sc2 `_ +`cisd_sc2 `_ CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) .br dets_in : bitmasks corresponding to determinants @@ -83,29 +83,29 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`repeat_excitation `_ +`repeat_excitation `_ Undocumented -`connected_to_ref `_ +`connected_to_ref `_ Undocumented -`det_is_not_or_may_be_in_ref `_ +`det_is_not_or_may_be_in_ref `_ If true, det is not in ref If false, det may be in ref -`is_in_wavefunction `_ +`is_in_wavefunction `_ Undocumented -`key_pattern_not_in_ref `_ +`key_pattern_not_in_ref `_ Min and max values of the integers of the keys of the reference -`davidson_converged `_ +`davidson_converged `_ True if the Davidson algorithm is converged -`davidson_criterion `_ +`davidson_criterion `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`davidson_diag `_ +`davidson_diag `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -123,7 +123,7 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_diag_hjj `_ +`davidson_diag_hjj `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -143,114 +143,117 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_iter_max `_ +`davidson_iter_max `_ Max number of Davidson iterations -`davidson_sze_max `_ +`davidson_sze_max `_ Max number of Davidson sizes -`davidson_threshold `_ +`davidson_threshold `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`det_search_key `_ +`det_search_key `_ Return an integer*8 corresponding to a determinant index for searching -`n_det `_ +`n_det `_ Number of determinants in the wave function -`n_det_max_jacobi `_ +`n_det_max_jacobi `_ Maximum number of determinants diagonalized my jacobi -`n_states `_ +`n_states `_ Number of states to consider -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`psi_det `_ +`psi_det `_ The wave function determinants. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`double_exc_bitmask `_ +`double_exc_bitmask `_ double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1 double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1 double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2 double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2 for a given couple of hole/particle excitations i. -`n_double_exc_bitmasks `_ +`n_double_exc_bitmasks `_ Number of double excitation bitmasks -`n_single_exc_bitmasks `_ +`n_single_exc_bitmasks `_ Number of single excitation bitmasks -`single_exc_bitmask `_ +`single_exc_bitmask `_ single_exc_bitmask(:,1,i) is the bitmask for holes single_exc_bitmask(:,2,i) is the bitmask for particles for a given couple of hole/particle excitations i. -`ci_eigenvectors `_ +`ci_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_electronic_energy `_ +`ci_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_energy `_ +`ci_energy `_ N_states lowest eigenvalues of the CI matrix -`diag_algorithm `_ +`diag_algorithm `_ Diagonalization algorithm (Davidson or Lapack) -`diagonalize_ci `_ +`diagonalize_ci `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`ci_sc2_eigenvectors `_ +`ci_sc2_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_sc2_electronic_energy `_ +`ci_sc2_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_sc2_energy `_ +`ci_sc2_energy `_ N_states lowest eigenvalues of the CI matrix -`diagonalize_ci_sc2 `_ +`diagonalize_ci_sc2 `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`filter_connected `_ +`threshold_convergence_sc2 `_ + convergence of the correlation energy of SC2 iterations + +`filter_connected `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -261,7 +264,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_davidson `_ +`filter_connected_davidson `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -272,7 +275,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0 `_ +`filter_connected_i_h_psi0 `_ returns the array idx which contains the index of the .br determinants in the array key1 that interact @@ -281,7 +284,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0_sc2 `_ +`filter_connected_i_h_psi0_sc2 `_ standard filter_connected_i_H_psi but returns in addition .br the array of the index of the non connected determinants to key1 @@ -292,69 +295,69 @@ Documentation .br to repeat the excitations -`get_s2 `_ +`get_s2 `_ Returns -`get_s2_u0 `_ +`get_s2_u0 `_ Undocumented -`s_z `_ +`s_z `_ Undocumented -`s_z2_sz `_ +`s_z2_sz `_ Undocumented -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem -`decode_exc `_ +`decode_exc `_ Decodes the exc arrays returned by get_excitation. h1,h2 : Holes p1,p2 : Particles s1,s2 : Spins (1:alpha, 2:beta) degree : Degree of excitation -`det_connections `_ +`det_connections `_ .br -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes -`get_double_excitation `_ +`get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase -`get_excitation `_ +`get_excitation `_ Returns the excitation operators between two determinants and the phase -`get_excitation_degree `_ +`get_excitation_degree `_ Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants -`get_mono_excitation `_ +`get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants .br H_jj : array of -`i_h_j `_ +`i_h_j `_ Returns where i and j are determinants -`i_h_psi `_ +`i_h_psi `_ for the various Nstates -`i_h_psi_sc2 `_ +`i_h_psi_sc2 `_ for the various Nstate .br returns in addition @@ -367,10 +370,10 @@ Documentation .br to repeat the excitations -`n_con_int `_ +`n_con_int `_ Number of integers to represent the connections between determinants -`h_matrix_all_dets `_ +`h_matrix_all_dets `_ H matrix on the basis of the slater deter;inants defined by psi_det diff --git a/src/MonoInts/NEEDED_MODULES b/src/MonoInts/NEEDED_MODULES index 56587081..190f8c6e 100644 --- a/src/MonoInts/NEEDED_MODULES +++ b/src/MonoInts/NEEDED_MODULES @@ -1 +1 @@ -AOs Ezfio_files MOs Nuclei Output Utils +AOs Electrons Ezfio_files MOs Nuclei Output Utils diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index a9bd07bc..d10cee82 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -5,6 +5,7 @@ Needed Modules .. NEEDED_MODULES file. * `AOs `_ +* `Electrons `_ * `Ezfio_files `_ * `MOs `_ * `Nuclei `_ From f3fc0fdb8abd70cd2beeccd0db0ab159236e3769 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 20 Jun 2014 18:35:26 +0200 Subject: [PATCH 9/9] Natural orbitals implemented --- src/DensityMatrix/density_matrix.irp.f | 48 ---------- src/Dets/H_apply_template.f | 2 +- src/Dets/density_matrix.irp.f | 106 +++++++++++++++++++++++ src/Dets/save_natorb.irp.f | 4 + src/Full_CI/full_ci.irp.f | 18 ++-- src/Hartree_Fock/README.rst | 61 ++++++------- src/Hartree_Fock/SCF.irp.f | 5 +- src/Hartree_Fock/mo_SCF_iterations.irp.f | 11 ++- src/MOs/utils.irp.f | 16 +++- 9 files changed, 177 insertions(+), 94 deletions(-) create mode 100644 src/Dets/density_matrix.irp.f create mode 100644 src/Dets/save_natorb.irp.f diff --git a/src/DensityMatrix/density_matrix.irp.f b/src/DensityMatrix/density_matrix.irp.f index 5afcb523..3fc1d6f9 100644 --- a/src/DensityMatrix/density_matrix.irp.f +++ b/src/DensityMatrix/density_matrix.irp.f @@ -209,51 +209,3 @@ END_PROVIDER enddo enddo END_PROVIDER - - BEGIN_PROVIDER [ double precision, one_body_dm_a, (mo_tot_num_align,mo_tot_num) ] -&BEGIN_PROVIDER [ double precision, one_body_dm_b, (mo_tot_num_align,mo_tot_num) ] - implicit none - BEGIN_DOC - ! Alpha and beta one-body density matrix - END_DOC - - integer :: j,k,l - integer :: occ(N_int*bit_kind_size,2) - 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 - one_body_dm_a = 0.d0 - one_body_dm_b = 0.d0 - - do k=1,det_num - call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int) - call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int) - ck = det_coef_provider(k) - do l=1,elec_alpha_num - j = occ(l,1) - one_body_dm_a(j,j) += ck*ck - enddo - do l=1,elec_beta_num - j = occ(l,2) - one_body_dm_b(j,j) += ck*ck - enddo - do l=1,k-1 - call get_excitation_degree(det_provider(1,1,k),det_provider(1,1,l),degree,N_int) - if (degree /= 1) then - cycle - endif - call get_mono_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - ckl = ck * det_coef_provider(l) * phase - if (s1==1) then - one_body_dm_a(h1,p1) += ckl - one_body_dm_a(p1,h1) += ckl - else - one_body_dm_b(h1,p1) += ckl - one_body_dm_b(p1,h1) += ckl - endif - enddo - enddo -END_PROVIDER - diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index b795223f..1cb5d4e7 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -389,7 +389,7 @@ subroutine $subroutine($params_main) !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i_generator,wall_2,ispin,k,mask) allocate( mask(N_int,2,6) ) - !$OMP DO SCHEDULE(guided) + !$OMP DO SCHEDULE(dynamic,4) do i_generator=1,nmax if (abort_here) then cycle diff --git a/src/Dets/density_matrix.irp.f b/src/Dets/density_matrix.irp.f new file mode 100644 index 00000000..e5a907fb --- /dev/null +++ b/src/Dets/density_matrix.irp.f @@ -0,0 +1,106 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix for each state + END_DOC + + integer :: j,k,l,m + integer :: occ(N_int*bit_kind_size,2) + 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 + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + + one_body_dm_mo_alpha = 0.d0 + 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 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) + allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) + tmp_a = 0.d0 + 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) + do m=1,N_states + ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j) += ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j) += ck + enddo + enddo + do l=1,k-1 + call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m) + if (s1==1) then + tmp_a(h1,p1) += ckl + tmp_a(p1,h1) += ckl + else + tmp_b(h1,p1) += ckl + tmp_b(p1,h1) += ckl + endif + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a + !$OMP END CRITICAL + !$OMP CRITICAL + one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b + !$OMP END CRITICAL + deallocate(tmp_a,tmp_b) + !$OMP BARRIER + !$OMP END PARALLEL +END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! One-body density matrix + END_DOC + one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta +END_PROVIDER + +subroutine save_natural_mos + implicit none + BEGIN_DOC + ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + 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) + deallocate(tmp) + call save_mos + +end + +BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights in the state-average calculation of the density matrix + END_DOC + state_average_weight = 1.d0/dble(N_states) +END_PROVIDER + diff --git a/src/Dets/save_natorb.irp.f b/src/Dets/save_natorb.irp.f new file mode 100644 index 00000000..dd3d7684 --- /dev/null +++ b/src/Dets/save_natorb.irp.f @@ -0,0 +1,4 @@ +program save_natorb + call save_natural_mos +end + diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f index 532d3788..49a25b76 100644 --- a/src/Full_CI/full_ci.irp.f +++ b/src/Full_CI/full_ci.irp.f @@ -25,14 +25,14 @@ program cisd exit endif enddo - print *, 'Final step' - call remove_small_contributions - call diagonalize_CI - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 - print *, '-----' +! print *, 'Final step' +! call remove_small_contributions +! call diagonalize_CI +! print *, 'N_det = ', N_det +! print *, 'N_states = ', N_states +! print *, 'PT2 = ', pt2 +! print *, 'E = ', CI_energy +! print *, 'E+PT2 = ', CI_energy+pt2 +! print *, '-----' deallocate(pt2,norm_pert) end diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index a9d70bd4..341c6284 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -26,22 +26,22 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fock_matrix_alpha_ao `_ +`fock_matrix_alpha_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_alpha_mo `_ +`fock_matrix_alpha_mo `_ Fock matrix on the MO basis -`fock_matrix_ao `_ +`fock_matrix_ao `_ Fock matrix in AO basis set -`fock_matrix_beta_ao `_ +`fock_matrix_beta_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_beta_mo `_ +`fock_matrix_beta_mo `_ Fock matrix on the MO basis -`fock_matrix_diag_mo `_ +`fock_matrix_diag_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -56,7 +56,7 @@ Documentation K = Fb - Fa .br -`fock_matrix_mo `_ +`fock_matrix_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -71,70 +71,73 @@ Documentation K = Fb - Fa .br -`hf_energy `_ +`hf_energy `_ Hartree-Fock energy -`hf_density_matrix_ao `_ +`hf_density_matrix_ao `_ Density matrix in the AO basis -`hf_density_matrix_ao_alpha `_ +`hf_density_matrix_ao_alpha `_ Alpha density matrix in the AO basis -`hf_density_matrix_ao_beta `_ +`hf_density_matrix_ao_beta `_ Beta density matrix in the AO basis -`fock_mo_to_ao `_ +`fock_mo_to_ao `_ Undocumented -`insert_new_scf_density_matrix `_ +`insert_new_scf_density_matrix `_ Undocumented -`it_scf `_ +`it_scf `_ Number of the current SCF iteration -`scf_density_matrices `_ +`scf_density_matrices `_ Density matrices at every SCF iteration -`scf_energies `_ +`scf_energies `_ Density matrices at every SCF iteration -`scf_interpolation_step `_ +`scf_interpolation_step `_ Undocumented -`scf_iterations `_ +`scf_iterations `_ Undocumented -`diagonal_fock_matrix_mo `_ +`diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`eigenvectors_fock_matrix_mo `_ +`eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`xcf_iteration `_ +`run `_ Undocumented -`do_diis `_ +`scf `_ + Undocumented + +`do_diis `_ If True, compute integrals on the fly -`n_it_scf_max `_ +`n_it_scf_max `_ Maximum number of SCF iterations -`thresh_scf `_ +`thresh_scf `_ Threshold on the convergence of the Hartree Fock energy -`bi_elec_ref_bitmask_energy `_ +`bi_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`kinetic_ref_bitmask_energy `_ +`kinetic_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`mono_elec_ref_bitmask_energy `_ +`mono_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`nucl_elec_ref_bitmask_energy `_ +`nucl_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`ref_bitmask_energy `_ +`ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f index 895ca546..650aa365 100644 --- a/src/Hartree_Fock/SCF.irp.f +++ b/src/Hartree_Fock/SCF.irp.f @@ -70,10 +70,11 @@ subroutine SCF_interpolation_step return endif call random_number(c) + c = c*0.5d0 do j=1,ao_num do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf)+SCF_density_matrices(i,j,1,it_scf-1) * (1.d0 - c) - HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf)+SCF_density_matrices(i,j,2,it_scf-1) * (1.d0 - c) + HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf-1)+SCF_density_matrices(i,j,1,it_scf) * (1.d0 - c) + HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf-1)+SCF_density_matrices(i,j,2,it_scf) * (1.d0 - c) enddo enddo TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 2698e010..55c42178 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -1,4 +1,11 @@ -program xcf_iteration + +program scf + call orthonormalize_mos + call run +end + +subroutine run + use bitmasks implicit none double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral @@ -8,7 +15,7 @@ program xcf_iteration E0 = HF_energy thresh_SCF = 1.d-10 - call scf_iterations + call damping_SCF mo_label = "Canonical" TOUCH mo_label mo_coef call save_mos diff --git a/src/MOs/utils.irp.f b/src/MOs/utils.irp.f index ee04ac4b..5e7984ff 100644 --- a/src/MOs/utils.irp.f +++ b/src/MOs/utils.irp.f @@ -28,21 +28,31 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label) double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R + call write_time(output_mos) if (m /= mo_tot_num) then print *, irp_here, ': Error : m/= mo_tot_num' + stop 1 endif - allocate(R(n,m)) - allocate(mo_coef_new(ao_num_align,m),eigvalues(m)) + allocate(R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m)) mo_coef_new = mo_coef call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2)) integer :: i + write (output_mos,'(A)'), 'MOs are now **'//trim(label)//'**' + write (output_mos,'(A)'), '' + write (output_mos,'(A)'), 'Eigenvalues' + write (output_mos,'(A)'), '-----------' + write (output_mos,'(A)'), '' + write (output_mos,'(A)'), '======== ================' do i = 1, m - print*,'eigvalues(i) = ',eigvalues(i) + write (output_mos,'(I8,X,F16.10)'), i,eigvalues(i) enddo + write (output_mos,'(A)'), '======== ================' + write (output_mos,'(A)'), '' call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1)) deallocate(mo_coef_new,R,eigvalues) + call write_time(output_mos) mo_label = label SOFT_TOUCH mo_coef mo_label