From fe3c2bb875f822eda66af4c411b239c72e94338e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 3 Apr 2015 14:26:14 +0200 Subject: [PATCH] Minor changes --- scripts/ezfio_with_default.py | 1 + src/CISD/cisd.irp.f | 1 + src/MRCC/H_apply.irp.f | 10 +- src/MRCC/mrcc_dress.irp.f | 197 ++++++++++++++-------------------- src/MRCC/mrcc_utils.irp.f | 132 +++++++++++++++++++++++ 5 files changed, 223 insertions(+), 118 deletions(-) diff --git a/scripts/ezfio_with_default.py b/scripts/ezfio_with_default.py index 1208f4b7..63666abd 100755 --- a/scripts/ezfio_with_default.py +++ b/scripts/ezfio_with_default.py @@ -94,6 +94,7 @@ END_PROVIDER file = open(filename,'r') lines = file.readlines() file.close() + k=-1 # Search directory for k,line in enumerate(lines): if line[0] != ' ': diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index 0b6bc9fd..6d310d95 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -13,6 +13,7 @@ program cisd print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy enddo + call save_wavefunction ! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) ! do i = 1, N_states ! print*,'eigvalues(i) = ',eigvalues(i) diff --git a/src/MRCC/H_apply.irp.f b/src/MRCC/H_apply.irp.f index 4f284112..a8519c25 100644 --- a/src/MRCC/H_apply.irp.f +++ b/src/MRCC/H_apply.irp.f @@ -2,13 +2,13 @@ use bitmasks BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * -s = H_apply("mrcc") +s = H_apply("mrcc_simple") s.data["parameters"] = ", delta_ij_sd_, Ndet_sd" s.data["declarations"] += """ integer, intent(in) :: Ndet_sd double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) """ -s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" +s.data["keys_work"] = "call mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" s.data["params_post"] += ", delta_ij_sd_, Ndet_sd" s.data["params_main"] += "delta_ij_sd_, Ndet_sd" s.data["decls_main"] += """ @@ -19,7 +19,13 @@ s.data["finalization"] = "" s.data["copy_buffer"] = "" s.data["generate_psi_guess"] = "" s.data["size_max"] = "256" +print s + + + +s.data["subroutine"] = "H_apply_mrcc" +s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" print s END_SHELL diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index e5c8d233..d9c8978c 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -1,3 +1,84 @@ +subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet_sd + double precision, intent(inout) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + integer :: new_size + logical :: is_in_wavefunction + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) + logical :: good + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq, c_ref + integer :: connected_to_ref + + N_tq = 0 + do i=1,N_selected + c_ref = connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, & + i_generator,N_det_generators) + + if (c_ref /= 0) then + cycle + endif + + ! Select determinants that are triple or quadruple excitations + ! from the CAS + good = .True. + call get_excitation_degree_vector(psi_cas,det_buffer(1,1,i),degree,Nint,N_det_cas,idx) + do k=1,idx(0) + if (degree(k) < 3) then + good = .False. + exit + endif + enddo + if (good) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + endif + enddo + + ! Compute / (E0 - Haa) + double precision :: hka, haa + double precision :: haj + double precision :: f(N_states) + + do i=1,N_tq + call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx) + call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa) + do m=1,N_states + f(m) = 1.d0/(ci_electronic_energy(m)-haa) + enddo + do k=1,idx(0) + call i_h_j(tq(1,1,i),psi_sd(1,1,idx(k)),Nint,hka) + do j=k,idx(0) + call i_h_j(tq(1,1,i),psi_sd(1,1,idx(j)),Nint,haj) + do m=1,N_states + delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m) + delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m) + enddo + enddo + enddo + enddo +end + + + + + + + + subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none @@ -72,119 +153,3 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin enddo end -BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ] - implicit none - BEGIN_DOC - ! Dressing matrix in SD basis - END_DOC - delta_ij_sd = 0.d0 - call H_apply_mrcc(delta_ij_sd,N_det_sd) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] - implicit none - BEGIN_DOC - ! Dressing matrix in N_det basis - END_DOC - integer :: i,j,m - delta_ij = 0.d0 - do m=1,N_states - do j=1,N_det_sd - do i=1,N_det_sd - delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m) - enddo - enddo - enddo -END_PROVIDER - -BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] - implicit none - BEGIN_DOC - ! Dressed H with Delta_ij - END_DOC - integer :: i, j - do j=1,N_det - do i=1,N_det - h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1) - if (i==j) then - print *, i, delta_ij(i,j,1), h_matrix_all_dets(i,j) - endif - enddo - enddo - -END_PROVIDER - - - BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! Eigenvectors/values of the CI matrix - END_DOC - integer :: i,j - - do j=1,N_states_diag - do i=1,N_det - CI_eigenvectors_dressed(i,j) = psi_coef(i,j) - enddo - enddo - - if (diag_algorithm == "Davidson") then - - stop 'use Lapack' -! call davidson_diag(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, & -! size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_Dets) - - else if (diag_algorithm == "Lapack") then - - double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) - allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) - allocate (eigenvalues(N_det)) - call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_dressed,size(H_matrix_dressed,1),N_det) - CI_electronic_energy_dressed(:) = 0.d0 - do i=1,N_det - CI_eigenvectors_dressed(i,1) = eigenvectors(i,1) - enddo - integer :: i_state - double precision :: s2 - i_state = 0 - do j=1,N_det - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - if(dabs(s2-expected_s2).le.0.3d0)then - i_state += 1 - do i=1,N_det - CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(i_state) = eigenvalues(j) - CI_eigenvectors_s2_dressed(i_state) = s2 - endif - if (i_state.ge.N_states_diag) then - exit - endif - enddo - deallocate(eigenvectors,eigenvalues) - endif - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! N_states lowest eigenvalues of the dressed CI matrix - END_DOC - - integer :: j - character*(8) :: st - call write_time(output_Dets) - do j=1,N_states_diag - CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion - write(st,'(I4)') j - call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st)) - call write_double(output_Dets,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) - enddo - -END_PROVIDER - diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 2c5dc811..01e566ae 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -148,3 +148,135 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] enddo END_PROVIDER + + +BEGIN_PROVIDER [ character*(32), dressing_type ] + implicit none + BEGIN_DOC + ! [ Simple | MRCC ] + END_DOC + dressing_type = "MRCC" +END_PROVIDER + +BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ] + implicit none + BEGIN_DOC + ! Dressing matrix in SD basis + END_DOC + delta_ij_sd = 0.d0 + if (dressing_type == "MRCC") then + call H_apply_mrcc(delta_ij_sd,N_det_sd) + else if (dressing_type == "Simple") then + call H_apply_mrcc_simple(delta_ij_sd,N_det_sd) + else + print *, irp_here + stop 'dressing' + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] + implicit none + BEGIN_DOC + ! Dressing matrix in N_det basis + END_DOC + integer :: i,j,m + delta_ij = 0.d0 + do m=1,N_states + do j=1,N_det_sd + do i=1,N_det_sd + delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m) + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] + implicit none + BEGIN_DOC + ! Dressed H with Delta_ij + END_DOC + integer :: i, j + do j=1,N_det + do i=1,N_det + h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1) + if (i==j) then + print *, i, delta_ij(i,j,1), h_matrix_all_dets(i,j) + endif + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] + implicit none + BEGIN_DOC + ! Eigenvectors/values of the CI matrix + END_DOC + integer :: i,j + + do j=1,N_states_diag + do i=1,N_det + CI_eigenvectors_dressed(i,j) = psi_coef(i,j) + enddo + enddo + + if (diag_algorithm == "Davidson") then + + stop 'use Lapack' +! call davidson_diag(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, & +! size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_Dets) + + else if (diag_algorithm == "Lapack") then + + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_dressed,size(H_matrix_dressed,1),N_det) + CI_electronic_energy_dressed(:) = 0.d0 + do i=1,N_det + CI_eigenvectors_dressed(i,1) = eigenvectors(i,1) + enddo + integer :: i_state + double precision :: s2 + i_state = 0 + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(s2-expected_s2).le.0.3d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(i_state) = eigenvalues(j) + CI_eigenvectors_s2_dressed(i_state) = s2 + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the dressed CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(output_Dets) + do j=1,N_states_diag + CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion + write(st,'(I4)') j + call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st)) + call write_double(output_Dets,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) + enddo + +END_PROVIDER +