diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 2c5e260b..99b0beef 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -10,7 +10,7 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ] else diag_algorithm = "Lapack" endif - + if (N_det < N_states) then diag_algorithm = "Lapack" endif @@ -18,152 +18,152 @@ END_PROVIDER BEGIN_PROVIDER [ integer, N_det ] - implicit none - BEGIN_DOC - ! Number of determinants in the wave function - END_DOC - logical :: exists - character*(64) :: label - PROVIDE read_wf mo_label ezfio_filename nproc - if (mpi_master) then - if (read_wf) then - call ezfio_has_determinants_n_det(exists) - if (exists) then - call ezfio_has_determinants_mo_label(exists) - if (exists) then - call ezfio_get_determinants_mo_label(label) - exists = (label == mo_label) - endif - endif - if (exists) then - call ezfio_get_determinants_n_det(N_det) - else - N_det = 1 - endif - else - N_det = 1 - endif - call write_int(output_determinants,N_det,'Number of determinants') - endif + implicit none + BEGIN_DOC + ! Number of determinants in the wave function + END_DOC + logical :: exists + character*(64) :: label + PROVIDE read_wf mo_label ezfio_filename nproc + if (mpi_master) then + if (read_wf) then + call ezfio_has_determinants_n_det(exists) + if (exists) then + call ezfio_has_determinants_mo_label(exists) + if (exists) then + call ezfio_get_determinants_mo_label(label) + exists = (label == mo_label) + endif + endif + if (exists) then + call ezfio_get_determinants_n_det(N_det) + else + N_det = 1 + endif + else + N_det = 1 + endif + call write_int(output_determinants,N_det,'Number of determinants') + endif IRP_IF MPI include 'mpif.h' - integer :: ierr + integer :: ierr call MPI_BCAST( N_det, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) if (ierr /= MPI_SUCCESS) then stop 'Unable to read N_det with MPI' endif IRP_ENDIF - - ASSERT (N_det > 0) + + ASSERT (N_det > 0) END_PROVIDER BEGIN_PROVIDER [integer, max_degree_exc] - implicit none - integer :: i,degree - max_degree_exc = 0 - BEGIN_DOC - ! Maximum degree of excitation in the wf - END_DOC - do i = 1, N_det - call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) - if(degree.gt.max_degree_exc)then - max_degree_exc= degree - endif - enddo + implicit none + integer :: i,degree + max_degree_exc = 0 + BEGIN_DOC + ! Maximum degree of excitation in the wf + END_DOC + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree.gt.max_degree_exc)then + max_degree_exc= degree + endif + enddo END_PROVIDER BEGIN_PROVIDER [ integer, psi_det_size ] - implicit none - BEGIN_DOC - ! Size of the psi_det/psi_coef arrays - END_DOC - PROVIDE ezfio_filename output_determinants - logical :: exists - if (mpi_master) then - call ezfio_has_determinants_n_det(exists) - if (exists) then - call ezfio_get_determinants_n_det(psi_det_size) - else - psi_det_size = 1 + implicit none + BEGIN_DOC + ! Size of the psi_det/psi_coef arrays + END_DOC + PROVIDE ezfio_filename output_determinants + logical :: exists + if (mpi_master) then + call ezfio_has_determinants_n_det(exists) + if (exists) then + call ezfio_get_determinants_n_det(psi_det_size) + else + psi_det_size = 1 + endif + psi_det_size = max(psi_det_size,100000) + call write_int(output_determinants,psi_det_size,'Dimension of the psi arrays') endif - psi_det_size = max(psi_det_size,100000) - call write_int(output_determinants,psi_det_size,'Dimension of the psi arrays') - endif IRP_IF MPI include 'mpif.h' - integer :: ierr + integer :: ierr call MPI_BCAST( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) if (ierr /= MPI_SUCCESS) then stop 'Unable to read psi_det_size with MPI' endif IRP_ENDIF - - + + END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ] - implicit none - BEGIN_DOC - ! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file - ! is empty - END_DOC + implicit none + BEGIN_DOC + ! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file + ! is empty + END_DOC integer :: i logical :: exists - character*(64) :: label + character*(64) :: label - PROVIDE read_wf N_det mo_label ezfio_filename HF_bitmask mo_coef + PROVIDE read_wf N_det mo_label ezfio_filename HF_bitmask mo_coef psi_det = 0_bit_kind if (mpi_master) then - if (read_wf) then - call ezfio_has_determinants_N_int(exists) - if (exists) then - call ezfio_has_determinants_bit_kind(exists) - if (exists) then - call ezfio_has_determinants_N_det(exists) + if (read_wf) then + call ezfio_has_determinants_N_int(exists) if (exists) then - call ezfio_has_determinants_N_states(exists) - if (exists) then - call ezfio_has_determinants_psi_det(exists) + call ezfio_has_determinants_bit_kind(exists) if (exists) then - call ezfio_has_determinants_mo_label(exists) + call ezfio_has_determinants_N_det(exists) if (exists) then - call ezfio_get_determinants_mo_label(label) - exists = (label == mo_label) + call ezfio_has_determinants_N_states(exists) + if (exists) then + call ezfio_has_determinants_psi_det(exists) + if (exists) then + call ezfio_has_determinants_mo_label(exists) + if (exists) then + call ezfio_get_determinants_mo_label(label) + exists = (label == mo_label) + endif + endif + endif endif endif - endif endif - endif - endif - - if (exists) then - call read_dets(psi_det,N_int,N_det) - print *, 'Read psi_det' + + if (exists) then + call read_dets(psi_det,N_int,N_det) + print *, 'Read psi_det' + else + psi_det = 0_bit_kind + do i=1,N_int + psi_det(i,1,1) = HF_bitmask(i,1) + psi_det(i,2,1) = HF_bitmask(i,2) + enddo + endif else - psi_det = 0_bit_kind - do i=1,N_int - psi_det(i,1,1) = HF_bitmask(i,1) - psi_det(i,2,1) = HF_bitmask(i,2) - enddo + psi_det = 0_bit_kind + do i=1,N_int + psi_det(i,1,1) = HF_bitmask(i,1) + psi_det(i,2,1) = HF_bitmask(i,2) + enddo endif - else - psi_det = 0_bit_kind - do i=1,N_int - psi_det(i,1,1) = HF_bitmask(i,1) - psi_det(i,2,1) = HF_bitmask(i,2) - enddo - endif endif IRP_IF MPI include 'mpif.h' - integer :: ierr + integer :: ierr call MPI_BCAST( psi_det, N_int*2*N_det, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) if (ierr /= MPI_SUCCESS) then stop 'Unable to read psi_det with MPI' endif IRP_ENDIF - - + + END_PROVIDER @@ -178,60 +178,60 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ] integer :: i,k, N_int2 logical :: exists character*(64) :: label - - PROVIDE read_wf N_det mo_label ezfio_filename + + PROVIDE read_wf N_det mo_label ezfio_filename psi_coef = 0.d0 do i=1,min(N_states,psi_det_size) psi_coef(i,i) = 1.d0 enddo - + if (mpi_master) then - if (read_wf) then - call ezfio_has_determinants_psi_coef(exists) - if (exists) then - call ezfio_has_determinants_mo_label(exists) + if (read_wf) then + call ezfio_has_determinants_psi_coef(exists) if (exists) then - call ezfio_get_determinants_mo_label(label) - exists = (label == mo_label) + call ezfio_has_determinants_mo_label(exists) + if (exists) then + call ezfio_get_determinants_mo_label(label) + exists = (label == mo_label) + endif endif + if (exists) then + call ezfio_get_determinants_psi_coef(psi_coef) + endif + endif - if (exists) then - call ezfio_get_determinants_psi_coef(psi_coef) - endif - - endif - print *, 'Read psi_coef' + print *, 'Read psi_coef' endif IRP_IF MPI include 'mpif.h' - integer :: ierr + integer :: ierr call MPI_BCAST( psi_coef, N_states*psi_det_size, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) if (ierr /= MPI_SUCCESS) then stop 'Unable to read psi_coef with MPI' endif IRP_ENDIF - - + + END_PROVIDER subroutine update_psi_average_norm_contrib(w) implicit none BEGIN_DOC -! Compute psi_average_norm_contrib for different state average weights w(:) + ! Compute psi_average_norm_contrib for different state average weights w(:) END_DOC - double precision, intent(in) :: w(N_states) - double precision :: w0(N_states), f + double precision, intent(in) :: w(N_states) + double precision :: w0(N_states), f w0(:) = w(:)/sum(w(:)) - - integer :: i,j,k + + integer :: i,j,k do i=1,N_det psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*w(1) enddo do k=2,N_states do i=1,N_det - psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + & - psi_coef(i,k)*psi_coef(i,k)*w(k) + psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + & + psi_coef(i,k)*psi_coef(i,k)*w(k) enddo enddo f = 1.d0/sum(psi_average_norm_contrib(1:N_det)) @@ -239,31 +239,31 @@ subroutine update_psi_average_norm_contrib(w) psi_average_norm_contrib(i) = psi_average_norm_contrib(i)*f enddo SOFT_TOUCH psi_average_norm_contrib - + end subroutine BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ] - implicit none - BEGIN_DOC - ! Contribution of determinants to the state-averaged density - END_DOC - integer :: i,j,k - double precision :: f - f = 1.d0/dble(N_states) - do i=1,N_det - psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*f - enddo - do k=2,N_states - do i=1,N_det - psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + & - psi_coef(i,k)*psi_coef(i,k)*f - enddo - enddo - f = 1.d0/sum(psi_average_norm_contrib(1:N_det)) - do i=1,N_det - psi_average_norm_contrib(i) = psi_average_norm_contrib(i)*f - enddo + implicit none + BEGIN_DOC + ! Contribution of determinants to the state-averaged density + END_DOC + integer :: i,j,k + double precision :: f + f = 1.d0/dble(N_states) + do i=1,N_det + psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*f + enddo + do k=2,N_states + do i=1,N_det + psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + & + psi_coef(i,k)*psi_coef(i,k)*f + enddo + enddo + f = 1.d0/sum(psi_average_norm_contrib(1:N_det)) + do i=1,N_det + psi_average_norm_contrib(i) = psi_average_norm_contrib(i)*f + enddo END_PROVIDER @@ -279,134 +279,134 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, psi_coef_sorted, (psi_det_size,N_states) ] &BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted, (psi_det_size) ] &BEGIN_PROVIDER [ integer, psi_det_sorted_order, (psi_det_size) ] - implicit none - BEGIN_DOC - ! Wave function sorted by determinants contribution to the norm (state-averaged) - ! - ! psi_det_sorted_order(i) -> k : index in psi_det - END_DOC - integer :: i,j,k - integer, allocatable :: iorder(:) - allocate ( iorder(N_det) ) - do i=1,N_det - psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib(i) - iorder(i) = i - enddo - call dsort(psi_average_norm_contrib_sorted,iorder,N_det) - do i=1,N_det - do j=1,N_int - psi_det_sorted(j,1,i) = psi_det(j,1,iorder(i)) - psi_det_sorted(j,2,i) = psi_det(j,2,iorder(i)) - enddo - do k=1,N_states - psi_coef_sorted(i,k) = psi_coef(iorder(i),k) - enddo - psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib_sorted(i) - enddo - do i=1,N_det - psi_det_sorted_order(iorder(i)) = i - enddo - - - deallocate(iorder) - + implicit none + BEGIN_DOC +! Wave function sorted by determinants contribution to the norm (state-averaged) + ! + ! psi_det_sorted_order(i) -> k : index in psi_det + END_DOC + integer :: i,j,k + integer, allocatable :: iorder(:) + allocate ( iorder(N_det) ) + do i=1,N_det + psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib(i) + iorder(i) = i + enddo + call dsort(psi_average_norm_contrib_sorted,iorder,N_det) + do i=1,N_det + do j=1,N_int + psi_det_sorted(j,1,i) = psi_det(j,1,iorder(i)) + psi_det_sorted(j,2,i) = psi_det(j,2,iorder(i)) + enddo + do k=1,N_states + psi_coef_sorted(i,k) = psi_coef(iorder(i),k) + enddo + psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib_sorted(i) + enddo + do i=1,N_det + psi_det_sorted_order(iorder(i)) = i + enddo + + + deallocate(iorder) + END_PROVIDER - + subroutine flip_generators() - integer :: i,j,k - integer(bit_kind) :: detmp(N_int,2) - double precision :: tmp(N_states) - - do i=1,N_det_generators/2 - detmp(:,:) = psi_det_sorted(:,:,i) - tmp = psi_coef_sorted(i, :) - psi_det_sorted(:,:,i) = psi_det_sorted(:,:,N_det_generators+1-i) - psi_coef_sorted(i, :) = psi_coef_sorted(N_det_generators+1-i, :) - - psi_det_sorted(:,:,N_det_generators+1-i) = detmp(:,:) - psi_coef_sorted(N_det_generators+1-i, :) = tmp - end do - - TOUCH psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted + integer :: i,j,k + integer(bit_kind) :: detmp(N_int,2) + double precision :: tmp(N_states) + + do i=1,N_det_generators/2 + detmp(:,:) = psi_det_sorted(:,:,i) + tmp = psi_coef_sorted(i, :) + psi_det_sorted(:,:,i) = psi_det_sorted(:,:,N_det_generators+1-i) + psi_coef_sorted(i, :) = psi_coef_sorted(N_det_generators+1-i, :) + + psi_det_sorted(:,:,N_det_generators+1-i) = detmp(:,:) + psi_coef_sorted(N_det_generators+1-i, :) = tmp + end do + + TOUCH psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted end subroutine - + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_bit, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ double precision, psi_coef_sorted_bit, (psi_det_size,N_states) ] - implicit none - BEGIN_DOC - ! Determinants on which we apply for perturbation. - ! They are sorted by determinants interpreted as integers. Useful - ! to accelerate the search of a random determinant in the wave - ! function. - END_DOC - - call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, & - psi_det_sorted_bit, psi_coef_sorted_bit) + implicit none + BEGIN_DOC + ! Determinants on which we apply for perturbation. + ! They are sorted by determinants interpreted as integers. Useful + ! to accelerate the search of a random determinant in the wave + ! function. + END_DOC + + call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, & + psi_det_sorted_bit, psi_coef_sorted_bit) END_PROVIDER - + subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out) - use bitmasks - implicit none - integer, intent(in) :: Ndet - integer(bit_kind), intent(in) :: det_in (N_int,2,psi_det_size) - double precision , intent(in) :: coef_in(psi_det_size,N_states) - integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size) - double precision , intent(out) :: coef_out(psi_det_size,N_states) - BEGIN_DOC - ! Determinants are sorted are sorted according to their det_search_key. - ! Useful to accelerate the search of a random determinant in the wave - ! function. - END_DOC - integer :: i,j,k - integer, allocatable :: iorder(:) - integer*8, allocatable :: bit_tmp(:) - integer*8, external :: det_search_key - - allocate ( iorder(Ndet), bit_tmp(Ndet) ) - - do i=1,Ndet - iorder(i) = i - !$DIR FORCEINLINE - bit_tmp(i) = det_search_key(det_in(1,1,i),N_int) - enddo - call i8sort(bit_tmp,iorder,Ndet) - !DIR$ IVDEP - do i=1,Ndet - do j=1,N_int - det_out(j,1,i) = det_in(j,1,iorder(i)) - det_out(j,2,i) = det_in(j,2,iorder(i)) - enddo - do k=1,N_states - coef_out(i,k) = coef_in(iorder(i),k) - enddo - enddo - - deallocate(iorder, bit_tmp) - + use bitmasks + implicit none + integer, intent(in) :: Ndet + integer(bit_kind), intent(in) :: det_in (N_int,2,psi_det_size) + double precision , intent(in) :: coef_in(psi_det_size,N_states) + integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size) + double precision , intent(out) :: coef_out(psi_det_size,N_states) + BEGIN_DOC + ! Determinants are sorted are sorted according to their det_search_key. + ! Useful to accelerate the search of a random determinant in the wave + ! function. + END_DOC + integer :: i,j,k + integer, allocatable :: iorder(:) + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: det_search_key + + allocate ( iorder(Ndet), bit_tmp(Ndet) ) + + do i=1,Ndet + iorder(i) = i + !$DIR FORCEINLINE + bit_tmp(i) = det_search_key(det_in(1,1,i),N_int) + enddo + call i8sort(bit_tmp,iorder,Ndet) + !DIR$ IVDEP + do i=1,Ndet + do j=1,N_int + det_out(j,1,i) = det_in(j,1,iorder(i)) + det_out(j,2,i) = det_in(j,2,iorder(i)) + enddo + do k=1,N_states + coef_out(i,k) = coef_in(iorder(i),k) + enddo + enddo + + deallocate(iorder, bit_tmp) + end - - - - BEGIN_PROVIDER [ double precision, psi_coef_max, (N_states) ] -&BEGIN_PROVIDER [ double precision, psi_coef_min, (N_states) ] -&BEGIN_PROVIDER [ double precision, abs_psi_coef_max, (N_states) ] -&BEGIN_PROVIDER [ double precision, abs_psi_coef_min, (N_states) ] - implicit none - BEGIN_DOC - ! Max and min values of the coefficients - END_DOC - integer:: i - do i=1,N_states - psi_coef_min(i) = minval(psi_coef(:,i)) - psi_coef_max(i) = maxval(psi_coef(:,i)) - abs_psi_coef_min(i) = minval( dabs(psi_coef(:,i)) ) - abs_psi_coef_max(i) = maxval( dabs(psi_coef(:,i)) ) - call write_double(6,psi_coef_max(i), 'Max coef') - call write_double(6,psi_coef_min(i), 'Min coef') - call write_double(6,abs_psi_coef_max(i), 'Max abs coef') - call write_double(6,abs_psi_coef_min(i), 'Min abs coef') - enddo - + + + + BEGIN_PROVIDER [ double precision, psi_coef_max, (N_states) ] +&BEGIN_PROVIDER [ double precision, psi_coef_min, (N_states) ] +&BEGIN_PROVIDER [ double precision, abs_psi_coef_max, (N_states) ] +&BEGIN_PROVIDER [ double precision, abs_psi_coef_min, (N_states) ] + implicit none + BEGIN_DOC + ! Max and min values of the coefficients + END_DOC + integer :: i + do i=1,N_states + psi_coef_min(i) = minval(psi_coef(:,i)) + psi_coef_max(i) = maxval(psi_coef(:,i)) + abs_psi_coef_min(i) = minval( dabs(psi_coef(:,i)) ) + abs_psi_coef_max(i) = maxval( dabs(psi_coef(:,i)) ) + call write_double(6,psi_coef_max(i), 'Max coef') + call write_double(6,psi_coef_min(i), 'Min coef') + call write_double(6,abs_psi_coef_max(i), 'Max abs coef') + call write_double(6,abs_psi_coef_min(i), 'Min abs coef') + enddo + END_PROVIDER @@ -460,9 +460,9 @@ subroutine read_dets(det,Nint,Ndet) end subroutine save_ref_determinant - implicit none + implicit none use bitmasks - double precision :: buffer(1,N_states) + double precision :: buffer(1,N_states) buffer = 0.d0 buffer(1,1) = 1.d0 call save_wavefunction_general(1,N_states,ref_bitmask,1,buffer) @@ -475,7 +475,7 @@ subroutine save_wavefunction implicit none use bitmasks BEGIN_DOC -! Save the wave function into the EZFIO file + ! Save the wave function into the EZFIO file END_DOC if (mpi_master) then call save_wavefunction_general(N_det,min(N_states,N_det),psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) @@ -487,7 +487,7 @@ subroutine save_wavefunction_unsorted implicit none use bitmasks BEGIN_DOC -! Save the wave function into the EZFIO file + ! Save the wave function into the EZFIO file END_DOC if (mpi_master) then call save_wavefunction_general(N_det,min(N_states,N_det),psi_det,size(psi_coef,1),psi_coef) @@ -497,60 +497,60 @@ end subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) implicit none BEGIN_DOC -! Save the wave function into the EZFIO file + ! Save the wave function into the EZFIO file END_DOC use bitmasks include 'constants.include.F' - integer, intent(in) :: ndet,nstates,dim_psicoef - integer(bit_kind), intent(in) :: psidet(N_int,2,ndet) - double precision, intent(in) :: psicoef(dim_psicoef,nstates) + integer, intent(in) :: ndet,nstates,dim_psicoef + integer(bit_kind), intent(in) :: psidet(N_int,2,ndet) + double precision, intent(in) :: psicoef(dim_psicoef,nstates) integer*8, allocatable :: psi_det_save(:,:,:) double precision, allocatable :: psi_coef_save(:,:) - - integer :: i,j,k - + + integer :: i,j,k + if (mpi_master) then - call ezfio_set_determinants_N_int(N_int) - call ezfio_set_determinants_bit_kind(bit_kind) - call ezfio_set_determinants_N_det(ndet) - call ezfio_set_determinants_n_states(nstates) - call ezfio_set_determinants_mo_label(mo_label) - - allocate (psi_det_save(N_int,2,ndet)) - do i=1,ndet - do j=1,2 - do k=1,N_int - psi_det_save(k,j,i) = transfer(psidet(k,j,i),1_8) - enddo - enddo - enddo - call ezfio_set_determinants_psi_det(psi_det_save) - deallocate (psi_det_save) - - allocate (psi_coef_save(ndet,nstates)) - double precision :: accu_norm(nstates) - accu_norm = 0.d0 - do k=1,nstates + call ezfio_set_determinants_N_int(N_int) + call ezfio_set_determinants_bit_kind(bit_kind) + call ezfio_set_determinants_N_det(ndet) + call ezfio_set_determinants_n_states(nstates) + call ezfio_set_determinants_mo_label(mo_label) + + allocate (psi_det_save(N_int,2,ndet)) do i=1,ndet - accu_norm(k) = accu_norm(k) + psicoef(i,k) * psicoef(i,k) - psi_coef_save(i,k) = psicoef(i,k) + do j=1,2 + do k=1,N_int + psi_det_save(k,j,i) = transfer(psidet(k,j,i),1_8) + enddo + enddo enddo - if (accu_norm(k) == 0.d0) then - accu_norm(k) = 1.e-12 - endif - enddo - do k = 1, nstates - accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) - enddo - do k=1,nstates - do i=1,ndet - psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k) + call ezfio_set_determinants_psi_det(psi_det_save) + deallocate (psi_det_save) + + allocate (psi_coef_save(ndet,nstates)) + double precision :: accu_norm(nstates) + accu_norm = 0.d0 + do k=1,nstates + do i=1,ndet + accu_norm(k) = accu_norm(k) + psicoef(i,k) * psicoef(i,k) + psi_coef_save(i,k) = psicoef(i,k) + enddo + if (accu_norm(k) == 0.d0) then + accu_norm(k) = 1.e-12 + endif enddo - enddo - - call ezfio_set_determinants_psi_coef(psi_coef_save) - deallocate (psi_coef_save) - call write_int(output_determinants,ndet,'Saved determinants') + do k = 1, nstates + accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) + enddo + do k=1,nstates + do i=1,ndet + psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k) + enddo + enddo + + call ezfio_set_determinants_psi_coef(psi_coef_save) + deallocate (psi_coef_save) + call write_int(output_determinants,ndet,'Saved determinants') endif end @@ -559,29 +559,29 @@ end subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,index_det_save) implicit none BEGIN_DOC -! Save the wave function into the EZFIO file + ! Save the wave function into the EZFIO file END_DOC use bitmasks - integer, intent(in) :: ndet,nstates - integer(bit_kind), intent(in) :: psidet(N_int,2,ndet) - double precision, intent(in) :: psicoef(ndet,nstates) - integer, intent(in) :: index_det_save(ndet) - integer, intent(in) :: ndetsave + integer, intent(in) :: ndet,nstates + integer(bit_kind), intent(in) :: psidet(N_int,2,ndet) + double precision, intent(in) :: psicoef(ndet,nstates) + integer, intent(in) :: index_det_save(ndet) + integer, intent(in) :: ndetsave integer*8, allocatable :: psi_det_save(:,:,:) double precision, allocatable :: psi_coef_save(:,:) integer*8 :: det_8(100) integer(bit_kind) :: det_bk((100*8)/bit_kind) integer :: N_int2 equivalence (det_8, det_bk) - - integer :: i,k - + + integer :: i,k + call ezfio_set_determinants_N_int(N_int) call ezfio_set_determinants_bit_kind(bit_kind) call ezfio_set_determinants_N_det(ndetsave) call ezfio_set_determinants_n_states(nstates) call ezfio_set_determinants_mo_label(mo_label) - + N_int2 = (N_int*bit_kind)/8 allocate (psi_det_save(N_int2,2,ndetsave)) do i=1,ndetsave @@ -600,11 +600,11 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde enddo call ezfio_set_determinants_psi_det(psi_det_save) deallocate (psi_det_save) - + progress_bar(1) = 7 progress_value = dble(progress_bar(1)) allocate (psi_coef_save(ndetsave,nstates)) - double precision :: accu_norm(nstates) + double precision :: accu_norm(nstates) accu_norm = 0.d0 do k=1,nstates do i=1,ndetsave @@ -613,14 +613,14 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde enddo enddo do k = 1, nstates - accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) + accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) enddo do k=1,nstates do i=1,ndetsave psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k) enddo enddo - + call ezfio_set_determinants_psi_coef(psi_coef_save) call write_int(output_determinants,ndet,'Saved determinants') deallocate (psi_coef_save) @@ -628,43 +628,43 @@ end logical function detEq(a,b,Nint) - use bitmasks - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) - integer :: ni, i - - detEq = .false. - do i=1,2 - do ni=1,Nint - if(a(ni,i) /= b(ni,i)) return - end do - end do - detEq = .true. + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) + integer :: ni, i + + detEq = .false. + do i=1,2 + do ni=1,Nint + if(a(ni,i) /= b(ni,i)) return + end do + end do + detEq = .true. end function integer function detCmp(a,b,Nint) - use bitmasks - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) - integer :: ni, i - - detCmp = 0 - do i=1,2 - do ni=Nint,1,-1 - - if(a(ni,i) < b(ni,i)) then - detCmp = -1 - return - else if(a(ni,i) > b(ni,i)) then - detCmp = 1 - return - end if - - end do - end do + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) + integer :: ni, i + + detCmp = 0 + do i=1,2 + do ni=Nint,1,-1 + + if(a(ni,i) < b(ni,i)) then + detCmp = -1 + return + else if(a(ni,i) > b(ni,i)) then + detCmp = 1 + return + end if + + end do + end do end function @@ -672,20 +672,20 @@ subroutine apply_excitation(det, exc, res, ok, Nint) use bitmasks implicit none - integer, intent(in) :: Nint - integer, intent(in) :: exc(0:2,2,2) - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: h1,p1,h2,p2,s1,s2,degree - integer :: ii, pos + integer, intent(in) :: Nint + integer, intent(in) :: exc(0:2,2,2) + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: h1,p1,h2,p2,s1,s2,degree + integer :: ii, pos ok = .false. degree = exc(0,1,1) + exc(0,1,2) -! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) -! INLINE + ! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + ! INLINE select case(degree) case(2) if (exc(0,1,1) == 2) then @@ -733,32 +733,32 @@ subroutine apply_excitation(det, exc, res, ok, Nint) p2 = 0 s1 = 0 s2 = 0 - case default + case default print *, degree print *, "apply ex" STOP end select -! END INLINE - - res = det + ! END INLINE - ii = ishft(h1-1,-bit_kind_shift) + 1 + res = det + + ii = ishft(h1-1,-bit_kind_shift) + 1 pos = h1-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s1), ibset(0_bit_kind, pos)) == 0_8) return res(ii, s1) = ibclr(res(ii, s1), pos) - ii = ishft(p1-1,-bit_kind_shift) + 1 + ii = ishft(p1-1,-bit_kind_shift) + 1 pos = p1-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return res(ii, s1) = ibset(res(ii, s1), pos) if(degree == 2) then - ii = ishft(h2-1,-bit_kind_shift) + 1 + ii = ishft(h2-1,-bit_kind_shift) + 1 pos = h2-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return res(ii, s2) = ibclr(res(ii, s2), pos) - ii = ishft(p2-1,-bit_kind_shift) + 1 + ii = ishft(p2-1,-bit_kind_shift) + 1 pos = p2-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return res(ii, s2) = ibset(res(ii, s2), pos) @@ -770,28 +770,28 @@ end subroutine subroutine apply_particles(det, s1, p1, s2, p2, res, ok, Nint) use bitmasks implicit none - integer, intent(in) :: Nint - integer, intent(in) :: s1, p1, s2, p2 - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: ii, pos + integer, intent(in) :: Nint + integer, intent(in) :: s1, p1, s2, p2 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos ok = .false. - res = det + res = det if(p1 /= 0) then - ii = ishft(p1-1,-bit_kind_shift) + 1 - pos = p1-1-ishft(ii-1,bit_kind_shift) - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return - res(ii, s1) = ibset(res(ii, s1), pos) + ii = ishft(p1-1,-bit_kind_shift) + 1 + pos = p1-1-ishft(ii-1,bit_kind_shift) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s1) = ibset(res(ii, s1), pos) end if - - ii = ishft(p2-1,-bit_kind_shift) + 1 + + ii = ishft(p2-1,-bit_kind_shift) + 1 pos = p2-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return res(ii, s2) = ibset(res(ii, s2), pos) - + ok = .true. end subroutine @@ -799,49 +799,49 @@ end subroutine subroutine apply_holes(det, s1, h1, s2, h2, res, ok, Nint) use bitmasks implicit none - integer, intent(in) :: Nint - integer, intent(in) :: s1, h1, s2, h2 - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: ii, pos + integer, intent(in) :: Nint + integer, intent(in) :: s1, h1, s2, h2 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos ok = .false. - res = det + res = det if(h1 /= 0) then - ii = ishft(h1-1,-bit_kind_shift) + 1 - pos = h1-1-ishft(ii-1,bit_kind_shift) - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return - res(ii, s1) = ibclr(res(ii, s1), pos) + ii = ishft(h1-1,-bit_kind_shift) + 1 + pos = h1-1-ishft(ii-1,bit_kind_shift) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s1) = ibclr(res(ii, s1), pos) end if - - ii = ishft(h2-1,-bit_kind_shift) + 1 + + ii = ishft(h2-1,-bit_kind_shift) + 1 pos = h2-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return res(ii, s2) = ibclr(res(ii, s2), pos) - + ok = .true. end subroutine subroutine apply_particle(det, s1, p1, res, ok, Nint) use bitmasks implicit none - integer, intent(in) :: Nint - integer, intent(in) :: s1, p1 - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: ii, pos + integer, intent(in) :: Nint + integer, intent(in) :: s1, p1 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos ok = .false. - res = det + res = det - ii = ishft(p1-1,-bit_kind_shift) + 1 + ii = ishft(p1-1,-bit_kind_shift) + 1 pos = p1-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return res(ii, s1) = ibset(res(ii, s1), pos) - + ok = .true. end subroutine @@ -849,20 +849,20 @@ end subroutine subroutine apply_hole(det, s1, h1, res, ok, Nint) use bitmasks implicit none - integer, intent(in) :: Nint - integer, intent(in) :: s1, h1 - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: ii, pos + integer, intent(in) :: Nint + integer, intent(in) :: s1, h1 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos ok = .false. - res = det + res = det - ii = ishft(h1-1,-bit_kind_shift) + 1 + ii = ishft(h1-1,-bit_kind_shift) + 1 pos = h1-1-ishft(ii-1,bit_kind_shift) if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return res(ii, s1) = ibclr(res(ii, s1), pos) - + ok = .true. end subroutine diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index ef27cd01..5eb10b20 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -1,24 +1,24 @@ BEGIN_PROVIDER [ integer, mo_tot_num ] implicit none BEGIN_DOC -! Number of MOs + ! Number of MOs END_DOC - + logical :: has - PROVIDE ezfio_filename + PROVIDE ezfio_filename if (mpi_master) then call ezfio_has_mo_basis_mo_tot_num(has) endif IRP_IF MPI include 'mpif.h' - integer :: ierr + integer :: ierr call MPI_BCAST( has, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) if (ierr /= MPI_SUCCESS) then stop 'Unable to read mo_tot_num with MPI' endif IRP_ENDIF if (.not.has) then - mo_tot_num = ao_ortho_canonical_num + mo_tot_num = ao_ortho_canonical_num else if (mpi_master) then call ezfio_get_mo_basis_mo_tot_num(mo_tot_num) @@ -32,7 +32,7 @@ BEGIN_PROVIDER [ integer, mo_tot_num ] endif call write_int(6,mo_tot_num,'mo_tot_num') ASSERT (mo_tot_num > 0) - + END_PROVIDER diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index a382ef8d..307a9610 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -7,61 +7,61 @@ BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num,3) ] PROVIDE ezfio_filename nucl_label nucl_charge if (mpi_master) then - double precision, allocatable :: buffer(:,:) - nucl_coord = 0.d0 - allocate (buffer(nucl_num,3)) - buffer = 0.d0 - logical :: has - call ezfio_has_nuclei_nucl_coord(has) - if (.not.has) then - print *, irp_here - stop 1 - endif - call ezfio_get_nuclei_nucl_coord(buffer) - integer :: i,j - - do i=1,3 - do j=1,nucl_num - nucl_coord(j,i) = buffer(j,i) - enddo - enddo - deallocate(buffer) - - character*(64), parameter :: f = '(A16, 4(1X,F12.6))' - character*(64), parameter :: ft= '(A16, 4(1X,A12 ))' - double precision, parameter :: a0= 0.529177249d0 - - call write_time(output_Nuclei) - write(output_Nuclei,'(A)') '' - write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)' - write(output_Nuclei,'(A)') '===============================' - write(output_Nuclei,'(A)') '' - write(output_Nuclei,ft) & - '================','============','============','============','============' - write(output_Nuclei,*) & - ' Atom Charge X Y Z ' - write(output_Nuclei,ft) & - '================','============','============','============','============' - do i=1,nucl_num - write(output_Nuclei,f) nucl_label(i), nucl_charge(i), & - nucl_coord(i,1)*a0, & - nucl_coord(i,2)*a0, & - nucl_coord(i,3)*a0 - enddo - write(output_Nuclei,ft) & - '================','============','============','============','============' - write(output_Nuclei,'(A)') '' - + double precision, allocatable :: buffer(:,:) + nucl_coord = 0.d0 + allocate (buffer(nucl_num,3)) + buffer = 0.d0 + logical :: has + call ezfio_has_nuclei_nucl_coord(has) + if (.not.has) then + print *, irp_here + stop 1 + endif + call ezfio_get_nuclei_nucl_coord(buffer) + integer :: i,j + + do i=1,3 + do j=1,nucl_num + nucl_coord(j,i) = buffer(j,i) + enddo + enddo + deallocate(buffer) + + character*(64), parameter :: f = '(A16, 4(1X,F12.6))' + character*(64), parameter :: ft= '(A16, 4(1X,A12 ))' + double precision, parameter :: a0= 0.529177249d0 + + call write_time(output_Nuclei) + write(output_Nuclei,'(A)') '' + write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)' + write(output_Nuclei,'(A)') '===============================' + write(output_Nuclei,'(A)') '' + write(output_Nuclei,ft) & + '================','============','============','============','============' + write(output_Nuclei,*) & + ' Atom Charge X Y Z ' + write(output_Nuclei,ft) & + '================','============','============','============','============' + do i=1,nucl_num + write(output_Nuclei,f) nucl_label(i), nucl_charge(i), & + nucl_coord(i,1)*a0, & + nucl_coord(i,2)*a0, & + nucl_coord(i,3)*a0 + enddo + write(output_Nuclei,ft) & + '================','============','============','============','============' + write(output_Nuclei,'(A)') '' + endif - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( nucl_coord, 3*nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read nucl_coord with MPI' - endif - IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( nucl_coord, 3*nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read nucl_coord with MPI' + endif + IRP_ENDIF END_PROVIDER @@ -146,29 +146,30 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] END_DOC PROVIDE mpi_master nucl_coord nucl_charge nucl_num - IF (disk_access_nuclear_repulsion.EQ.'Read') THEN - LOGICAL :: has - if (mpi_master) then - call ezfio_has_nuclei_nuclear_repulsion(has) - if (has) then - call ezfio_get_nuclei_nuclear_repulsion(nuclear_repulsion) - else - print *, 'nuclei/nuclear_repulsion not found in EZFIO file' - stop 1 - endif - print*, 'Read nuclear_repulsion' - endif - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read nuclear_repulsion with MPI' - endif - IRP_ENDIF + if (disk_access_nuclear_repulsion.EQ.'Read') then + logical :: has - - ELSE + if (mpi_master) then + call ezfio_has_nuclei_nuclear_repulsion(has) + if (has) then + call ezfio_get_nuclei_nuclear_repulsion(nuclear_repulsion) + else + print *, 'nuclei/nuclear_repulsion not found in EZFIO file' + stop 1 + endif + print*, 'Read nuclear_repulsion' + endif + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( nuclear_repulsion, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read nuclear_repulsion with MPI' + endif + IRP_ENDIF + + + else integer :: k,l double precision :: Z12, r2, x(3) @@ -187,17 +188,17 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] enddo enddo nuclear_repulsion *= 0.5d0 - END IF + end if call write_time(output_Nuclei) call write_double(output_Nuclei,nuclear_repulsion, & 'Nuclear repulsion energy') - IF (disk_access_nuclear_repulsion.EQ.'Write') THEN + if (disk_access_nuclear_repulsion.EQ.'Write') then if (mpi_master) then call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) endif - END IF + endif END_PROVIDER BEGIN_PROVIDER [ character*(128), element_name, (78)]