diff --git a/INSTALL.rst b/INSTALL.rst index c91b184f..825c0c4d 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -143,7 +143,7 @@ IRPF90 to Parameters (IRP) method. * Download the latest version of IRPF90 - here : ``_ and move + here : ``_ and move the downloaded archive in the :file:`${QP_ROOT}/external` directory * Extract the archive and go into the :file:`irpf90-*` directory to run diff --git a/configure b/configure index c343c243..c3ff8942 100755 --- a/configure +++ b/configure @@ -297,12 +297,13 @@ EOF ${QP_ROOT}/bin + EOF rm ${QP_ROOT}/external/opam_installer.sh source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true - ${QP_ROOT}/bin/opam init --verbose --yes + ${QP_ROOT}/bin/opam init --verbose --yes --comp=4.07.1 --disable-sandboxing eval $(${QP_ROOT}/bin/opam env) opam install -y ${OCAML_PACKAGES} || exit 1 @@ -310,13 +311,14 @@ EOF # Conventional commands execute << EOF chmod +x "\${QP_ROOT}"/external/opam_installer.sh + "\${QP_ROOT}"/external/opam_installer.sh --no-backup rm --force \${QP_ROOT}/bin/opam export OPAMROOT=\${OPAMROOT:-\${QP_ROOT}/external/opam} echo \${QP_ROOT}/bin \ | sh \${QP_ROOT}/external/opam_installer.sh rm \${QP_ROOT}/external/opam_installer.sh source \${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true - \${QP_ROOT}/bin/opam init --verbose --yes + \${QP_ROOT}/bin/opam init --verbose --yes --comp=4.07.1 --disable-sandboxing eval \$(\${QP_ROOT}/bin/opam env) opam install -y \${OCAML_PACKAGES} || exit 1 EOF diff --git a/src/bitmask/bitmasks.irp.f b/src/bitmask/bitmasks.irp.f index d425dda6..bbcff63c 100644 --- a/src/bitmask/bitmasks.irp.f +++ b/src/bitmask/bitmasks.irp.f @@ -11,7 +11,7 @@ BEGIN_PROVIDER [ integer, N_int ] if (N_int > N_int_max) then stop 'N_int > N_int_max' endif - + END_PROVIDER @@ -20,7 +20,7 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask, (N_int) ] BEGIN_DOC ! Bitmask to include all possible MOs END_DOC - + integer :: i,j,k k=0 do j=1,N_int @@ -37,34 +37,34 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask_4, (N_int,4) ] implicit none - integer :: i + integer :: i do i=1,N_int - full_ijkl_bitmask_4(i,1) = full_ijkl_bitmask(i) - full_ijkl_bitmask_4(i,2) = full_ijkl_bitmask(i) - full_ijkl_bitmask_4(i,3) = full_ijkl_bitmask(i) - full_ijkl_bitmask_4(i,4) = full_ijkl_bitmask(i) + full_ijkl_bitmask_4(i,1) = full_ijkl_bitmask(i) + full_ijkl_bitmask_4(i,2) = full_ijkl_bitmask(i) + full_ijkl_bitmask_4(i,3) = full_ijkl_bitmask(i) + full_ijkl_bitmask_4(i,4) = full_ijkl_bitmask(i) enddo END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), core_inact_act_bitmask_4, (N_int,4) ] implicit none - integer :: i + integer :: i do i=1,N_int - core_inact_act_bitmask_4(i,1) = reunion_of_core_inact_act_bitmask(i,1) - core_inact_act_bitmask_4(i,2) = reunion_of_core_inact_act_bitmask(i,1) - core_inact_act_bitmask_4(i,3) = reunion_of_core_inact_act_bitmask(i,1) - core_inact_act_bitmask_4(i,4) = reunion_of_core_inact_act_bitmask(i,1) + core_inact_act_bitmask_4(i,1) = reunion_of_core_inact_act_bitmask(i,1) + core_inact_act_bitmask_4(i,2) = reunion_of_core_inact_act_bitmask(i,1) + core_inact_act_bitmask_4(i,3) = reunion_of_core_inact_act_bitmask(i,1) + core_inact_act_bitmask_4(i,4) = reunion_of_core_inact_act_bitmask(i,1) enddo END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask_4, (N_int,4) ] implicit none - integer :: i + integer :: i do i=1,N_int - virt_bitmask_4(i,1) = virt_bitmask(i,1) - virt_bitmask_4(i,2) = virt_bitmask(i,1) - virt_bitmask_4(i,3) = virt_bitmask(i,1) - virt_bitmask_4(i,4) = virt_bitmask(i,1) + virt_bitmask_4(i,1) = virt_bitmask(i,1) + virt_bitmask_4(i,2) = virt_bitmask(i,1) + virt_bitmask_4(i,3) = virt_bitmask(i,1) + virt_bitmask_4(i,4) = virt_bitmask(i,1) enddo END_PROVIDER @@ -78,491 +78,480 @@ BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask, (N_int,2)] END_DOC integer :: i,j,n integer :: occ(elec_alpha_num) - + HF_bitmask = 0_bit_kind do i=1,elec_alpha_num - occ(i) = i + occ(i) = i enddo call list_to_bitstring( HF_bitmask(1,1), occ, elec_alpha_num, N_int) ! elec_alpha_num <= elec_beta_num, so occ is already OK. call list_to_bitstring( HF_bitmask(1,2), occ, elec_beta_num, N_int) - + END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), ref_bitmask, (N_int,2)] - implicit none - BEGIN_DOC -! Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask - END_DOC - ref_bitmask = HF_bitmask + implicit none + BEGIN_DOC + ! Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask + END_DOC + ref_bitmask = HF_bitmask END_PROVIDER BEGIN_PROVIDER [ integer, N_generators_bitmask ] - implicit none - BEGIN_DOC - ! Number of bitmasks for generators - END_DOC - logical :: exists - PROVIDE ezfio_filename N_int - - if (mpi_master) then - call ezfio_has_bitmasks_N_mask_gen(exists) - if (exists) then - call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask) - integer :: N_int_check - integer :: bit_kind_check - call ezfio_get_bitmasks_bit_kind(bit_kind_check) - if (bit_kind_check /= bit_kind) then - print *, bit_kind_check, bit_kind - print *, 'Error: bit_kind is not correct in EZFIO file' + implicit none + BEGIN_DOC + ! Number of bitmasks for generators + END_DOC + logical :: exists + PROVIDE ezfio_filename N_int + + if (mpi_master) then + call ezfio_has_bitmasks_N_mask_gen(exists) + if (exists) then + call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask) + integer :: N_int_check + integer :: bit_kind_check + call ezfio_get_bitmasks_bit_kind(bit_kind_check) + if (bit_kind_check /= bit_kind) then + print *, bit_kind_check, bit_kind + print *, 'Error: bit_kind is not correct in EZFIO file' + endif + call ezfio_get_bitmasks_N_int(N_int_check) + if (N_int_check /= N_int) then + print *, N_int_check, N_int + print *, 'Error: N_int is not correct in EZFIO file' + endif + else + N_generators_bitmask = 1 endif - call ezfio_get_bitmasks_N_int(N_int_check) - if (N_int_check /= N_int) then - print *, N_int_check, N_int - print *, 'Error: N_int is not correct in EZFIO file' - endif - else - N_generators_bitmask = 1 + ASSERT (N_generators_bitmask > 0) + call write_int(6,N_generators_bitmask,'N_generators_bitmask') endif - ASSERT (N_generators_bitmask > 0) - call write_int(6,N_generators_bitmask,'N_generators_bitmask') - endif IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) IRP_ENDIF IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( N_generators_bitmask, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read N_generators_bitmask with MPI' - endif + include 'mpif.h' + integer :: ierr + call MPI_BCAST( N_generators_bitmask, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read N_generators_bitmask with MPI' + endif IRP_ENDIF - - + + END_PROVIDER BEGIN_PROVIDER [ integer, N_generators_bitmask_restart ] - implicit none - BEGIN_DOC - ! Number of bitmasks for generators - END_DOC - logical :: exists - PROVIDE ezfio_filename N_int - - if (mpi_master) then - call ezfio_has_bitmasks_N_mask_gen(exists) - if (exists) then - call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask_restart) - integer :: N_int_check - integer :: bit_kind_check - call ezfio_get_bitmasks_bit_kind(bit_kind_check) - if (bit_kind_check /= bit_kind) then - print *, bit_kind_check, bit_kind - print *, 'Error: bit_kind is not correct in EZFIO file' + implicit none + BEGIN_DOC + ! Number of bitmasks for generators + END_DOC + logical :: exists + PROVIDE ezfio_filename N_int + + if (mpi_master) then + call ezfio_has_bitmasks_N_mask_gen(exists) + if (exists) then + call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask_restart) + integer :: N_int_check + integer :: bit_kind_check + call ezfio_get_bitmasks_bit_kind(bit_kind_check) + if (bit_kind_check /= bit_kind) then + print *, bit_kind_check, bit_kind + print *, 'Error: bit_kind is not correct in EZFIO file' + endif + call ezfio_get_bitmasks_N_int(N_int_check) + if (N_int_check /= N_int) then + print *, N_int_check, N_int + print *, 'Error: N_int is not correct in EZFIO file' + endif + else + N_generators_bitmask_restart = 1 endif - call ezfio_get_bitmasks_N_int(N_int_check) - if (N_int_check /= N_int) then - print *, N_int_check, N_int - print *, 'Error: N_int is not correct in EZFIO file' - endif - else - N_generators_bitmask_restart = 1 + ASSERT (N_generators_bitmask_restart > 0) + call write_int(6,N_generators_bitmask_restart,'N_generators_bitmask_restart') endif - ASSERT (N_generators_bitmask_restart > 0) - call write_int(6,N_generators_bitmask_restart,'N_generators_bitmask_restart') - endif IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) IRP_ENDIF - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( N_generators_bitmask_restart, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read N_generators_bitmask_restart with MPI' - endif - IRP_ENDIF - - + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( N_generators_bitmask_restart, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read N_generators_bitmask_restart with MPI' + endif + IRP_ENDIF + + END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask_restart) ] - implicit none - BEGIN_DOC - ! Bitmasks for generator determinants. - ! (N_int, alpha/beta, hole/particle, generator). - ! - ! 3rd index is : - ! - ! * 1 : hole for single exc - ! - ! * 2 : particle for single exc - ! - ! * 3 : hole for 1st exc of double - ! - ! * 4 : particle for 1st exc of double - ! - ! * 5 : hole for 2nd exc of double - ! - ! * 6 : particle for 2nd exc of double - ! - END_DOC - logical :: exists - PROVIDE ezfio_filename full_ijkl_bitmask N_generators_bitmask N_int - PROVIDE generators_bitmask_restart - - if (mpi_master) then - call ezfio_has_bitmasks_generators(exists) - if (exists) then - call ezfio_get_bitmasks_generators(generators_bitmask_restart) - else - integer :: k, ispin + implicit none + BEGIN_DOC + ! Bitmasks for generator determinants. + ! (N_int, alpha/beta, hole/particle, generator). + ! + ! 3rd index is : + ! + ! * 1 : hole for single exc + ! + ! * 2 : particle for single exc + ! + ! * 3 : hole for 1st exc of double + ! + ! * 4 : particle for 1st exc of double + ! + ! * 5 : hole for 2nd exc of double + ! + ! * 6 : particle for 2nd exc of double + ! + END_DOC + logical :: exists + PROVIDE ezfio_filename full_ijkl_bitmask N_generators_bitmask N_int + PROVIDE generators_bitmask_restart + + if (mpi_master) then + call ezfio_has_bitmasks_generators(exists) + if (exists) then + call ezfio_get_bitmasks_generators(generators_bitmask_restart) + else + integer :: k, ispin + do k=1,N_generators_bitmask + do ispin=1,2 + do i=1,N_int + generators_bitmask_restart(i,ispin,s_hole ,k) = full_ijkl_bitmask(i) + generators_bitmask_restart(i,ispin,s_part ,k) = full_ijkl_bitmask(i) + generators_bitmask_restart(i,ispin,d_hole1,k) = full_ijkl_bitmask(i) + generators_bitmask_restart(i,ispin,d_part1,k) = full_ijkl_bitmask(i) + generators_bitmask_restart(i,ispin,d_hole2,k) = full_ijkl_bitmask(i) + generators_bitmask_restart(i,ispin,d_part2,k) = full_ijkl_bitmask(i) + enddo + enddo + enddo + endif + + integer :: i do k=1,N_generators_bitmask do ispin=1,2 do i=1,N_int - generators_bitmask_restart(i,ispin,s_hole ,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,s_part ,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_hole1,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_part1,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_hole2,k) = full_ijkl_bitmask(i) - generators_bitmask_restart(i,ispin,d_part2,k) = full_ijkl_bitmask(i) + generators_bitmask_restart(i,ispin,s_hole ,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,s_hole,k) ) + generators_bitmask_restart(i,ispin,s_part ,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,s_part,k) ) + generators_bitmask_restart(i,ispin,d_hole1,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_hole1,k) ) + generators_bitmask_restart(i,ispin,d_part1,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_part1,k) ) + generators_bitmask_restart(i,ispin,d_hole2,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_hole2,k) ) + generators_bitmask_restart(i,ispin,d_part2,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_part2,k) ) enddo enddo enddo endif - - integer :: i - do k=1,N_generators_bitmask - do ispin=1,2 - do i=1,N_int - generators_bitmask_restart(i,ispin,s_hole ,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,s_hole,k) ) - generators_bitmask_restart(i,ispin,s_part ,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,s_part,k) ) - generators_bitmask_restart(i,ispin,d_hole1,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_hole1,k) ) - generators_bitmask_restart(i,ispin,d_part1,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_part1,k) ) - generators_bitmask_restart(i,ispin,d_hole2,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_hole2,k) ) - generators_bitmask_restart(i,ispin,d_part2,k) = iand(full_ijkl_bitmask(i),generators_bitmask_restart(i,ispin,d_part2,k) ) - enddo - enddo - enddo - endif IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) IRP_ENDIF IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( generators_bitmask_restart, N_int*2*6*N_generators_bitmask_restart, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read generators_bitmask_restart with MPI' - endif + include 'mpif.h' + integer :: ierr + call MPI_BCAST( generators_bitmask_restart, N_int*2*6*N_generators_bitmask_restart, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read generators_bitmask_restart with MPI' + endif IRP_ENDIF - + END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_bitmask) ] - implicit none - BEGIN_DOC - ! Bitmasks for generator determinants. - ! (N_int, alpha/beta, hole/particle, generator). - ! - ! 3rd index is : - ! - ! * 1 : hole for single exc - ! - ! * 2 : particle for single exc - ! - ! * 3 : hole for 1st exc of double - ! - ! * 4 : particle for 1st exc of double - ! - ! * 5 : hole for 2nd exc of double - ! - ! * 6 : particle for 2nd exc of double - ! - END_DOC - logical :: exists - PROVIDE ezfio_filename full_ijkl_bitmask N_generators_bitmask - -if (mpi_master) then - call ezfio_has_bitmasks_generators(exists) - if (exists) then - call ezfio_get_bitmasks_generators(generators_bitmask) - else - integer :: k, ispin, i - do k=1,N_generators_bitmask - do ispin=1,2 - do i=1,N_int - generators_bitmask(i,ispin,s_hole ,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,s_part ,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_hole1,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_part1,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_hole2,k) = full_ijkl_bitmask(i) - generators_bitmask(i,ispin,d_part2,k) = full_ijkl_bitmask(i) + implicit none + BEGIN_DOC + ! Bitmasks for generator determinants. + ! (N_int, alpha/beta, hole/particle, generator). + ! + ! 3rd index is : + ! + ! * 1 : hole for single exc + ! + ! * 2 : particle for single exc + ! + ! * 3 : hole for 1st exc of double + ! + ! * 4 : particle for 1st exc of double + ! + ! * 5 : hole for 2nd exc of double + ! + ! * 6 : particle for 2nd exc of double + ! + END_DOC + logical :: exists + PROVIDE ezfio_filename full_ijkl_bitmask N_generators_bitmask + + if (mpi_master) then + call ezfio_has_bitmasks_generators(exists) + if (exists) then + call ezfio_get_bitmasks_generators(generators_bitmask) + else + integer :: k, ispin, i + do k=1,N_generators_bitmask + do ispin=1,2 + do i=1,N_int + generators_bitmask(i,ispin,s_hole ,k) = full_ijkl_bitmask(i) + generators_bitmask(i,ispin,s_part ,k) = full_ijkl_bitmask(i) + generators_bitmask(i,ispin,d_hole1,k) = full_ijkl_bitmask(i) + generators_bitmask(i,ispin,d_part1,k) = full_ijkl_bitmask(i) + generators_bitmask(i,ispin,d_hole2,k) = full_ijkl_bitmask(i) + generators_bitmask(i,ispin,d_part2,k) = full_ijkl_bitmask(i) + enddo + enddo enddo - enddo - enddo - endif - - do k=1,N_generators_bitmask - do ispin=1,2 - do i=1,N_int - generators_bitmask(i,ispin,s_hole ,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,s_hole,k) ) - generators_bitmask(i,ispin,s_part ,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,s_part,k) ) - generators_bitmask(i,ispin,d_hole1,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_hole1,k) ) - generators_bitmask(i,ispin,d_part1,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_part1,k) ) - generators_bitmask(i,ispin,d_hole2,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_hole2,k) ) - generators_bitmask(i,ispin,d_part2,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_part2,k) ) - enddo - enddo - enddo - endif + endif + + do k=1,N_generators_bitmask + do ispin=1,2 + do i=1,N_int + generators_bitmask(i,ispin,s_hole ,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,s_hole,k) ) + generators_bitmask(i,ispin,s_part ,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,s_part,k) ) + generators_bitmask(i,ispin,d_hole1,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_hole1,k) ) + generators_bitmask(i,ispin,d_part1,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_part1,k) ) + generators_bitmask(i,ispin,d_hole2,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_hole2,k) ) + generators_bitmask(i,ispin,d_part2,k) = iand(full_ijkl_bitmask(i),generators_bitmask(i,ispin,d_part2,k) ) + enddo + enddo + enddo + endif IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) IRP_ENDIF IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( generators_bitmask, N_int*2*6*N_generators_bitmask, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read generators_bitmask with MPI' - endif + include 'mpif.h' + integer :: ierr + call MPI_BCAST( generators_bitmask, N_int*2*6*N_generators_bitmask, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read generators_bitmask with MPI' + endif IRP_ENDIF - + END_PROVIDER BEGIN_PROVIDER [ integer, N_cas_bitmask ] - implicit none - BEGIN_DOC - ! Number of bitmasks for CAS - END_DOC - logical :: exists - PROVIDE ezfio_filename - PROVIDE N_cas_bitmask N_int - if (mpi_master) then - call ezfio_has_bitmasks_N_mask_cas(exists) - if (exists) then - call ezfio_get_bitmasks_N_mask_cas(N_cas_bitmask) - integer :: N_int_check - integer :: bit_kind_check - call ezfio_get_bitmasks_bit_kind(bit_kind_check) - if (bit_kind_check /= bit_kind) then - print *, bit_kind_check, bit_kind - print *, 'Error: bit_kind is not correct in EZFIO file' + implicit none + BEGIN_DOC + ! Number of bitmasks for CAS + END_DOC + logical :: exists + PROVIDE ezfio_filename + PROVIDE N_cas_bitmask N_int + if (mpi_master) then + call ezfio_has_bitmasks_N_mask_cas(exists) + if (exists) then + call ezfio_get_bitmasks_N_mask_cas(N_cas_bitmask) + integer :: N_int_check + integer :: bit_kind_check + call ezfio_get_bitmasks_bit_kind(bit_kind_check) + if (bit_kind_check /= bit_kind) then + print *, bit_kind_check, bit_kind + print *, 'Error: bit_kind is not correct in EZFIO file' + endif + call ezfio_get_bitmasks_N_int(N_int_check) + if (N_int_check /= N_int) then + print *, N_int_check, N_int + print *, 'Error: N_int is not correct in EZFIO file' + endif + else + N_cas_bitmask = 1 endif - call ezfio_get_bitmasks_N_int(N_int_check) - if (N_int_check /= N_int) then - print *, N_int_check, N_int - print *, 'Error: N_int is not correct in EZFIO file' - endif - else - N_cas_bitmask = 1 + call write_int(6,N_cas_bitmask,'N_cas_bitmask') endif - call write_int(6,N_cas_bitmask,'N_cas_bitmask') - endif - ASSERT (N_cas_bitmask > 0) + ASSERT (N_cas_bitmask > 0) IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) IRP_ENDIF IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( N_cas_bitmask, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read N_cas_bitmask with MPI' - endif + include 'mpif.h' + integer :: ierr + call MPI_BCAST( N_cas_bitmask, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read N_cas_bitmask with MPI' + endif IRP_ENDIF - + END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] - implicit none - BEGIN_DOC - ! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference) - END_DOC - logical :: exists - integer :: i,i_part,i_gen,j,k - PROVIDE ezfio_filename generators_bitmask_restart full_ijkl_bitmask - PROVIDE n_generators_bitmask HF_bitmask - - if (mpi_master) then - call ezfio_has_bitmasks_cas(exists) - if (exists) then - call ezfio_get_bitmasks_cas(cas_bitmask) - else - if(N_generators_bitmask == 1)then - do j=1, N_cas_bitmask - do i=1, N_int - cas_bitmask(i,1,j) = iand(not(HF_bitmask(i,1)),full_ijkl_bitmask(i)) - cas_bitmask(i,2,j) = iand(not(HF_bitmask(i,2)),full_ijkl_bitmask(i)) - enddo - enddo + implicit none + BEGIN_DOC + ! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference) + END_DOC + logical :: exists + integer :: i,i_part,i_gen,j,k + PROVIDE ezfio_filename generators_bitmask_restart full_ijkl_bitmask + PROVIDE n_generators_bitmask HF_bitmask + + if (mpi_master) then + call ezfio_has_bitmasks_cas(exists) + if (exists) then + call ezfio_get_bitmasks_cas(cas_bitmask) else - i_part = 2 - i_gen = 1 - do j=1, N_cas_bitmask - do i=1, N_int - cas_bitmask(i,1,j) = generators_bitmask_restart(i,1,i_part,i_gen) - cas_bitmask(i,2,j) = generators_bitmask_restart(i,2,i_part,i_gen) - enddo - enddo + if(N_generators_bitmask == 1)then + do j=1, N_cas_bitmask + do i=1, N_int + cas_bitmask(i,1,j) = iand(not(HF_bitmask(i,1)),full_ijkl_bitmask(i)) + cas_bitmask(i,2,j) = iand(not(HF_bitmask(i,2)),full_ijkl_bitmask(i)) + enddo + enddo + else + i_part = 2 + i_gen = 1 + do j=1, N_cas_bitmask + do i=1, N_int + cas_bitmask(i,1,j) = generators_bitmask_restart(i,1,i_part,i_gen) + cas_bitmask(i,2,j) = generators_bitmask_restart(i,2,i_part,i_gen) + enddo + enddo + endif endif - endif - do i=1,N_cas_bitmask - do j = 1, N_cas_bitmask - do k=1,N_int - cas_bitmask(k,j,i) = iand(cas_bitmask(k,j,i),full_ijkl_bitmask(k)) + do i=1,N_cas_bitmask + do j = 1, N_cas_bitmask + do k=1,N_int + cas_bitmask(k,j,i) = iand(cas_bitmask(k,j,i),full_ijkl_bitmask(k)) + enddo enddo enddo - enddo - write(*,*) 'Read CAS bitmask' - endif + write(*,*) 'Read CAS bitmask' + endif IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) IRP_ENDIF IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( cas_bitmask, N_int*2*N_cas_bitmask, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read cas_bitmask with MPI' - endif + include 'mpif.h' + integer :: ierr + call MPI_BCAST( cas_bitmask, N_int*2*N_cas_bitmask, MPI_BIT_KIND, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read cas_bitmask with MPI' + endif IRP_ENDIF - - + + END_PROVIDER - BEGIN_PROVIDER [ integer, n_core_inact_orb ] - implicit none - integer :: i - n_core_inact_orb = 0 - do i = 1, N_int - n_core_inact_orb += popcnt(reunion_of_core_inact_bitmask(i,1)) - enddo - ENd_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)] - implicit none - BEGIN_DOC - ! Reunion of the core and inactive and virtual bitmasks - END_DOC - integer :: i - do i = 1, N_int - reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1)) - reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2)) - enddo - END_PROVIDER +BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the core and inactive and virtual bitmasks + END_DOC + integer :: i + do i = 1, N_int + reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1)) + reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2)) + enddo +END_PROVIDER - BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)] - implicit none - BEGIN_DOC - ! Reunion of the core, inactive and active bitmasks - END_DOC - integer :: i,j - - do i = 1, N_int - reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),act_bitmask(i,1)) - reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),act_bitmask(i,2)) - enddo - END_PROVIDER +BEGIN_PROVIDER [integer(bit_kind), reunion_of_inact_act_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the inactive and active bitmasks + END_DOC + integer :: i,j + + do i = 1, N_int + reunion_of_inact_act_bitmask(i,1) = ior(inact_bitmask(i,1),act_bitmask(i,1)) + reunion_of_inact_act_bitmask(i,2) = ior(inact_bitmask(i,2),act_bitmask(i,2)) + enddo +END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] - implicit none - BEGIN_DOC - ! Reunion of the inactive, active and virtual bitmasks - END_DOC - integer :: i,j - do i = 1, N_int - reunion_of_bitmask(i,1) = ior(ior(cas_bitmask(i,1,1),inact_bitmask(i,1)),virt_bitmask(i,1)) - reunion_of_bitmask(i,2) = ior(ior(cas_bitmask(i,2,1),inact_bitmask(i,2)),virt_bitmask(i,2)) - enddo - END_PROVIDER +BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the core, inactive and active bitmasks + END_DOC + integer :: i,j + + do i = 1, N_int + reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),act_bitmask(i,1)) + reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),act_bitmask(i,2)) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the inactive, active and virtual bitmasks + END_DOC + integer :: i,j + do i = 1, N_int + reunion_of_bitmask(i,1) = ior(ior(cas_bitmask(i,1,1),inact_bitmask(i,1)),virt_bitmask(i,1)) + reunion_of_bitmask(i,2) = ior(ior(cas_bitmask(i,2,1),inact_bitmask(i,2)),virt_bitmask(i,2)) + enddo +END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), inact_virt_bitmask, (N_int,2)] &BEGIN_PROVIDER [ integer(bit_kind), core_inact_virt_bitmask, (N_int,2)] - implicit none - BEGIN_DOC - ! Reunion of the inactive and virtual bitmasks - END_DOC - integer :: i,j - do i = 1, N_int - inact_virt_bitmask(i,1) = ior(inact_bitmask(i,1),virt_bitmask(i,1)) - inact_virt_bitmask(i,2) = ior(inact_bitmask(i,2),virt_bitmask(i,2)) - core_inact_virt_bitmask(i,1) = ior(core_bitmask(i,1),inact_virt_bitmask(i,1)) - core_inact_virt_bitmask(i,2) = ior(core_bitmask(i,2),inact_virt_bitmask(i,2)) - enddo - END_PROVIDER + implicit none + BEGIN_DOC + ! Reunion of the inactive and virtual bitmasks + END_DOC + integer :: i,j + do i = 1, N_int + inact_virt_bitmask(i,1) = ior(inact_bitmask(i,1),virt_bitmask(i,1)) + inact_virt_bitmask(i,2) = ior(inact_bitmask(i,2),virt_bitmask(i,2)) + core_inact_virt_bitmask(i,1) = ior(core_bitmask(i,1),inact_virt_bitmask(i,1)) + core_inact_virt_bitmask(i,2) = ior(core_bitmask(i,2),inact_virt_bitmask(i,2)) + enddo +END_PROVIDER BEGIN_PROVIDER [ integer, i_bitmask_gen ] - implicit none - BEGIN_DOC - ! Current bitmask for the generators - END_DOC - i_bitmask_gen = 1 + implicit none + BEGIN_DOC + ! Current bitmask for the generators + END_DOC + i_bitmask_gen = 1 END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), unpaired_alpha_electrons, (N_int)] - implicit none - BEGIN_DOC - ! Bitmask reprenting the unpaired alpha electrons in the HF_bitmask - END_DOC - integer :: i - unpaired_alpha_electrons = 0_bit_kind - do i = 1, N_int - unpaired_alpha_electrons(i) = xor(HF_bitmask(i,1),HF_bitmask(i,2)) - enddo - END_PROVIDER +BEGIN_PROVIDER [ integer(bit_kind), unpaired_alpha_electrons, (N_int)] + implicit none + BEGIN_DOC + ! Bitmask reprenting the unpaired alpha electrons in the HF_bitmask + END_DOC + integer :: i + unpaired_alpha_electrons = 0_bit_kind + do i = 1, N_int + unpaired_alpha_electrons(i) = xor(HF_bitmask(i,1),HF_bitmask(i,2)) + enddo +END_PROVIDER - BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)] - implicit none - integer :: i,j - do i = 1, N_int - closed_shell_ref_bitmask(i,1) = ior(ref_bitmask(i,1),cas_bitmask(i,1,1)) - closed_shell_ref_bitmask(i,2) = ior(ref_bitmask(i,2),cas_bitmask(i,2,1)) - enddo - END_PROVIDER +BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)] + implicit none + integer :: i,j + do i = 1, N_int + closed_shell_ref_bitmask(i,1) = ior(ref_bitmask(i,1),cas_bitmask(i,1,1)) + closed_shell_ref_bitmask(i,2) = ior(ref_bitmask(i,2),cas_bitmask(i,2,1)) + enddo +END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), reunion_of_cas_inact_bitmask, (N_int,2)] - implicit none - BEGIN_DOC - ! Reunion of the inactive, active and virtual bitmasks - END_DOC - integer :: i,j - do i = 1, N_int - reunion_of_cas_inact_bitmask(i,1) = ior(act_bitmask(i,1),inact_bitmask(i,1)) - reunion_of_cas_inact_bitmask(i,2) = ior(act_bitmask(i,2),inact_bitmask(i,2)) - enddo - END_PROVIDER - - - BEGIN_PROVIDER [integer, n_core_orb_allocate] - implicit none - n_core_orb_allocate = max(n_core_orb,1) - END_PROVIDER - - BEGIN_PROVIDER [integer, n_inact_orb_allocate] - implicit none - n_inact_orb_allocate = max(n_inact_orb,1) - END_PROVIDER - - BEGIN_PROVIDER [integer, n_virt_orb_allocate] - implicit none - n_virt_orb_allocate = max(n_virt_orb,1) - END_PROVIDER +BEGIN_PROVIDER [ integer(bit_kind), reunion_of_cas_inact_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the inactive, active and virtual bitmasks + END_DOC + integer :: i,j + do i = 1, N_int + reunion_of_cas_inact_bitmask(i,1) = ior(act_bitmask(i,1),inact_bitmask(i,1)) + reunion_of_cas_inact_bitmask(i,2) = ior(act_bitmask(i,2),inact_bitmask(i,2)) + enddo +END_PROVIDER diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index 177c3df5..ff7ee2de 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -1,250 +1,382 @@ use bitmasks +BEGIN_PROVIDER [ integer, n_core_orb] + implicit none + BEGIN_DOC + ! Number of core MOs + END_DOC + integer :: i + + n_core_orb = 0 + do i = 1, mo_num + if(mo_class(i) == 'Core')then + n_core_orb += 1 + endif + enddo + + call write_int(6,n_core_orb, 'Number of core MOs') + +END_PROVIDER - BEGIN_PROVIDER [ integer, n_core_orb] - &BEGIN_PROVIDER [ integer, n_inact_orb ] - &BEGIN_PROVIDER [ integer, n_act_orb] - &BEGIN_PROVIDER [ integer, n_virt_orb ] - &BEGIN_PROVIDER [ integer, n_del_orb ] - implicit none - BEGIN_DOC - ! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited - ! in post CAS methods - ! n_inact_orb : Number of inactive orbitals - ! virt_bitmask : Bitmaks of vritual orbitals which are supposed to be recieve electrons - ! in post CAS methods - ! n_virt_orb : Number of virtual orbitals - ! list_inact : List of the inactive orbitals which are supposed to be doubly excited - ! in post CAS methods - ! list_virt : List of vritual orbitals which are supposed to be recieve electrons - ! in post CAS methods - ! list_inact_reverse : reverse list of inactive orbitals - ! list_inact_reverse(i) = 0 ::> not an inactive - ! list_inact_reverse(i) = k ::> IS the kth inactive - ! list_virt_reverse : reverse list of virtual orbitals - ! list_virt_reverse(i) = 0 ::> not an virtual - ! list_virt_reverse(i) = k ::> IS the kth virtual - ! list_act(i) = index of the ith active orbital - ! - ! list_act_reverse : reverse list of active orbitals - ! list_act_reverse(i) = 0 ::> not an active - ! list_act_reverse(i) = k ::> IS the kth active orbital - END_DOC - logical :: exists - integer :: j,i +BEGIN_PROVIDER [ integer, n_inact_orb ] + implicit none + BEGIN_DOC + ! Number of inactive MOs + END_DOC + integer :: i + + n_inact_orb = 0 + do i = 1, mo_num + if (mo_class(i) == 'Inactive')then + n_inact_orb += 1 + endif + enddo + + call write_int(6,n_inact_orb,'Number of inactive MOs') + +END_PROVIDER - n_core_orb = 0 - n_inact_orb = 0 - n_act_orb = 0 - n_virt_orb = 0 - n_del_orb = 0 - do i = 1, mo_num - if(mo_class(i) == 'Core')then - n_core_orb += 1 - else if (mo_class(i) == 'Inactive')then - n_inact_orb += 1 - else if (mo_class(i) == 'Active')then - n_act_orb += 1 - else if (mo_class(i) == 'Virtual')then - n_virt_orb += 1 - else if (mo_class(i) == 'Deleted')then - n_del_orb += 1 - endif - enddo +BEGIN_PROVIDER [ integer, n_act_orb] + implicit none + BEGIN_DOC + ! Number of active MOs + END_DOC + integer :: i + + n_act_orb = 0 + do i = 1, mo_num + if (mo_class(i) == 'Active')then + n_act_orb += 1 + endif + enddo + + call write_int(6,n_act_orb, 'Number of active MOs') + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_virt_orb ] + implicit none + BEGIN_DOC + ! Number of virtual MOs + END_DOC + integer :: i + + n_virt_orb = 0 + do i = 1, mo_num + if (mo_class(i) == 'Virtual')then + n_virt_orb += 1 + endif + enddo + + call write_int(6,n_virt_orb, 'Number of virtual MOs') + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_del_orb ] + implicit none + BEGIN_DOC + ! Number of deleted MOs + END_DOC + integer :: i + + n_del_orb = 0 + do i = 1, mo_num + if (mo_class(i) == 'Deleted')then + n_del_orb += 1 + endif + enddo + + call write_int(6,n_del_orb, 'Number of deleted MOs') + +END_PROVIDER - call write_int(6,n_core_orb, 'Number of core MOs') - call write_int(6,n_inact_orb,'Number of inactive MOs') - call write_int(6,n_act_orb, 'Number of active MOs') - call write_int(6,n_virt_orb, 'Number of virtual MOs') - call write_int(6,n_del_orb, 'Number of deleted MOs') - - END_PROVIDER - - - BEGIN_PROVIDER [integer, dim_list_core_orb] -&BEGIN_PROVIDER [integer, dim_list_inact_orb] -&BEGIN_PROVIDER [integer, dim_list_virt_orb] -&BEGIN_PROVIDER [integer, dim_list_act_orb] -&BEGIN_PROVIDER [integer, dim_list_del_orb] - implicit none - BEGIN_DOC -! dimensions for the allocation of list_inact, list_virt, list_core and list_act -! it is at least 1 - END_DOC - dim_list_core_orb = max(n_core_orb,1) - dim_list_inact_orb = max(n_inact_orb,1) - dim_list_virt_orb = max(n_virt_orb,1) - dim_list_act_orb = max(n_act_orb,1) - dim_list_del_orb = max(n_del_orb,1) -END_PROVIDER - - BEGIN_PROVIDER [ integer, list_inact, (dim_list_inact_orb)] -&BEGIN_PROVIDER [ integer, list_virt, (dim_list_virt_orb)] -&BEGIN_PROVIDER [ integer, list_inact_reverse, (mo_num)] -&BEGIN_PROVIDER [ integer, list_virt_reverse, (mo_num)] -&BEGIN_PROVIDER [ integer, list_del_reverse, (mo_num)] -&BEGIN_PROVIDER [ integer, list_del, (mo_num)] -&BEGIN_PROVIDER [integer, list_core, (dim_list_core_orb)] -&BEGIN_PROVIDER [integer, list_core_reverse, (mo_num)] -&BEGIN_PROVIDER [integer, list_act, (dim_list_act_orb)] -&BEGIN_PROVIDER [integer, list_act_reverse, (mo_num)] -&BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)] -&BEGIN_PROVIDER [ integer(bit_kind), inact_bitmask, (N_int,2) ] -&BEGIN_PROVIDER [ integer(bit_kind), act_bitmask, (N_int,2) ] -&BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask, (N_int,2) ] -&BEGIN_PROVIDER [ integer(bit_kind), del_bitmask, (N_int,2) ] - implicit none - BEGIN_DOC - ! inact_bitmask : Bitmask of the inactive orbitals which are supposed to be doubly excited - ! in post CAS methods - ! n_inact_orb : Number of inactive orbitals - ! virt_bitmask : Bitmaks of vritual orbitals which are supposed to be recieve electrons - ! in post CAS methods - ! n_virt_orb : Number of virtual orbitals - ! list_inact : List of the inactive orbitals which are supposed to be doubly excited - ! in post CAS methods - ! list_virt : List of vritual orbitals which are supposed to be recieve electrons - ! in post CAS methods - ! list_inact_reverse : reverse list of inactive orbitals - ! list_inact_reverse(i) = 0 ::> not an inactive - ! list_inact_reverse(i) = k ::> IS the kth inactive - ! list_virt_reverse : reverse list of virtual orbitals - ! list_virt_reverse(i) = 0 ::> not an virtual - ! list_virt_reverse(i) = k ::> IS the kth virtual - ! list_act(i) = index of the ith active orbital - ! - ! list_act_reverse : reverse list of active orbitals - ! list_act_reverse(i) = 0 ::> not an active - ! list_act_reverse(i) = k ::> IS the kth active orbital - END_DOC - logical :: exists - integer :: j,i - integer :: n_core_orb_tmp, n_inact_orb_tmp, n_act_orb_tmp, n_virt_orb_tmp,n_del_orb_tmp - integer :: list_core_tmp(N_int*bit_kind_size) - integer :: list_inact_tmp(N_int*bit_kind_size) - integer :: list_act_tmp(N_int*bit_kind_size) - integer :: list_virt_tmp(N_int*bit_kind_size) - integer :: list_del_tmp(N_int*bit_kind_size) - list_core = 0 - list_inact = 0 - list_act = 0 - list_virt = 0 - list_del = 0 - list_core_reverse = 0 - list_inact_reverse = 0 - list_act_reverse = 0 - list_virt_reverse = 0 - list_del_reverse = 0 - n_core_orb_tmp = 0 - n_inact_orb_tmp = 0 - n_act_orb_tmp = 0 - n_virt_orb_tmp = 0 - n_del_orb_tmp = 0 - core_bitmask = 0_bit_kind - inact_bitmask = 0_bit_kind - act_bitmask = 0_bit_kind - virt_bitmask = 0_bit_kind - do i = 1, mo_num - if(mo_class(i) == 'Core')then - n_core_orb_tmp += 1 - list_core(n_core_orb_tmp) = i - list_core_tmp(n_core_orb_tmp) = i - list_core_reverse(i) = n_core_orb_tmp - else if (mo_class(i) == 'Inactive')then - n_inact_orb_tmp += 1 - list_inact(n_inact_orb_tmp) = i - list_inact_tmp(n_inact_orb_tmp) = i - list_inact_reverse(i) = n_inact_orb_tmp - else if (mo_class(i) == 'Active')then - n_act_orb_tmp += 1 - list_act(n_act_orb_tmp) = i - list_act_tmp(n_act_orb_tmp) = i - list_act_reverse(i) = n_act_orb_tmp - else if (mo_class(i) == 'Virtual')then - n_virt_orb_tmp += 1 - list_virt(n_virt_orb_tmp) = i - list_virt_tmp(n_virt_orb_tmp) = i - list_virt_reverse(i) = n_virt_orb_tmp - else if (mo_class(i) == 'Deleted')then - n_del_orb_tmp += 1 - list_del(n_del_orb_tmp) = i - list_del_tmp(n_del_orb_tmp) = i - list_del_reverse(i) = n_del_orb_tmp - endif - enddo - - if(n_core_orb.ne.0)then - call list_to_bitstring( core_bitmask(1,1), list_core, n_core_orb, N_int) - call list_to_bitstring( core_bitmask(1,2), list_core, n_core_orb, N_int) - endif - if(n_inact_orb.ne.0)then - call list_to_bitstring( inact_bitmask(1,1), list_inact, n_inact_orb, N_int) - call list_to_bitstring( inact_bitmask(1,2), list_inact, n_inact_orb, N_int) - endif - if(n_act_orb.ne.0)then - call list_to_bitstring( act_bitmask(1,1), list_act, n_act_orb, N_int) - call list_to_bitstring( act_bitmask(1,2), list_act, n_act_orb, N_int) - endif - if(n_virt_orb.ne.0)then - call list_to_bitstring( virt_bitmask(1,1), list_virt, n_virt_orb, N_int) - call list_to_bitstring( virt_bitmask(1,2), list_virt, n_virt_orb, N_int) - endif - if(n_del_orb.ne.0)then - call list_to_bitstring( del_bitmask(1,1), list_del, n_del_orb, N_int) - call list_to_bitstring( del_bitmask(1,2), list_del, n_del_orb, N_int) - endif - - -END_PROVIDER +BEGIN_PROVIDER [ integer, n_core_inact_orb ] + implicit none + BEGIN_DOC + ! n_core + n_inact + END_DOC + integer :: i + n_core_inact_orb = 0 + do i = 1, N_int + n_core_inact_orb += popcnt(reunion_of_core_inact_bitmask(i,1)) + enddo +END_PROVIDER BEGIN_PROVIDER [integer, n_inact_act_orb ] - implicit none - n_inact_act_orb = (n_inact_orb+n_act_orb) + implicit none + BEGIN_DOC + ! n_inact + n_act + END_DOC + n_inact_act_orb = (n_inact_orb+n_act_orb) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_core_orb] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_core. + ! it is at least 1 + END_DOC + dim_list_core_orb = max(n_core_orb,1) +END_PROVIDER -END_PROVIDER +BEGIN_PROVIDER [integer, dim_list_inact_orb] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_inact. + ! it is at least 1 + END_DOC + dim_list_inact_orb = max(n_inact_orb,1) +END_PROVIDER -BEGIN_PROVIDER [integer, list_inact_act, (n_inact_act_orb)] - integer :: i,itmp - itmp = 0 - do i = 1, n_inact_orb - itmp += 1 - list_inact_act(itmp) = list_inact(i) - enddo - do i = 1, n_act_orb - itmp += 1 - list_inact_act(itmp) = list_act(i) - enddo -END_PROVIDER +BEGIN_PROVIDER [integer, dim_list_act_orb] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_act. + ! it is at least 1 + END_DOC + dim_list_act_orb = max(n_act_orb,1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_virt_orb] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_virt. + ! it is at least 1 + END_DOC + dim_list_virt_orb = max(n_virt_orb,1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_del_orb] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_del. + ! it is at least 1 + END_DOC + dim_list_del_orb = max(n_del_orb,1) +END_PROVIDER BEGIN_PROVIDER [integer, n_core_inact_act_orb ] - implicit none - n_core_inact_act_orb = (n_core_orb + n_inact_orb + n_act_orb) + implicit none + BEGIN_DOC + ! Number of core inactive and active MOs + END_DOC + n_core_inact_act_orb = (n_core_orb + n_inact_orb + n_act_orb) +END_PROVIDER + -END_PROVIDER - BEGIN_PROVIDER [integer, list_core_inact_act, (n_core_inact_act_orb)] -&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (n_core_inact_act_orb)] - integer :: i,itmp - itmp = 0 - do i = 1, n_core_orb - itmp += 1 - list_core_inact_act(itmp) = list_core(i) - enddo - do i = 1, n_inact_orb - itmp += 1 - list_core_inact_act(itmp) = list_inact(i) - enddo - do i = 1, n_act_orb - itmp += 1 - list_core_inact_act(itmp) = list_act(i) - enddo + + BEGIN_PROVIDER [ integer(bit_kind), core_bitmask , (N_int,2) ] +&BEGIN_PROVIDER [ integer(bit_kind), inact_bitmask, (N_int,2) ] +&BEGIN_PROVIDER [ integer(bit_kind), act_bitmask , (N_int,2) ] +&BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask , (N_int,2) ] +&BEGIN_PROVIDER [ integer(bit_kind), del_bitmask , (N_int,2) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the core/inactive/active/virtual/deleted MOs + END_DOC - integer :: occ_inact(N_int*bit_kind_size) - occ_inact = 0 - call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), occ_inact(1), itest, N_int) - list_inact_reverse = 0 - do i = 1, n_core_inact_act_orb - list_core_inact_act_reverse(occ_inact(i)) = i - enddo -END_PROVIDER + core_bitmask = 0_bit_kind + inact_bitmask = 0_bit_kind + act_bitmask = 0_bit_kind + virt_bitmask = 0_bit_kind + del_bitmask = 0_bit_kind + + if(n_core_orb > 0)then + call list_to_bitstring( core_bitmask(1,1), list_core, n_core_orb, N_int) + call list_to_bitstring( core_bitmask(1,2), list_core, n_core_orb, N_int) + endif + if(n_inact_orb > 0)then + call list_to_bitstring( inact_bitmask(1,1), list_inact, n_inact_orb, N_int) + call list_to_bitstring( inact_bitmask(1,2), list_inact, n_inact_orb, N_int) + endif + if(n_act_orb > 0)then + call list_to_bitstring( act_bitmask(1,1), list_act, n_act_orb, N_int) + call list_to_bitstring( act_bitmask(1,2), list_act, n_act_orb, N_int) + endif + if(n_virt_orb > 0)then + call list_to_bitstring( virt_bitmask(1,1), list_virt, n_virt_orb, N_int) + call list_to_bitstring( virt_bitmask(1,2), list_virt, n_virt_orb, N_int) + endif + if(n_del_orb > 0)then + call list_to_bitstring( del_bitmask(1,1), list_del, n_del_orb, N_int) + call list_to_bitstring( del_bitmask(1,2), list_del, n_del_orb, N_int) + endif + +END_PROVIDER + + + + + + BEGIN_PROVIDER [ integer, list_core , (dim_list_core_orb) ] +&BEGIN_PROVIDER [ integer, list_core_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are in the core. + END_DOC + integer :: i, n + list_core = 0 + list_core_reverse = 0 + + n=0 + do i = 1, mo_num + if(mo_class(i) == 'Core')then + n += 1 + list_core(n) = i + list_core_reverse(i) = n + endif + enddo + print *, 'Core MOs:' + print *, list_core(1:n_core_orb) + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_inact , (dim_list_inact_orb) ] +&BEGIN_PROVIDER [ integer, list_inact_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are inactive. + END_DOC + integer :: i, n + list_inact = 0 + list_inact_reverse = 0 + + n=0 + do i = 1, mo_num + if (mo_class(i) == 'Inactive')then + n += 1 + list_inact(n) = i + list_inact_reverse(i) = n + endif + enddo + print *, 'Inactive MOs:' + print *, list_inact(1:n_inact_orb) + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_virt , (dim_list_virt_orb) ] +&BEGIN_PROVIDER [ integer, list_virt_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are virtual + END_DOC + integer :: i, n + list_virt = 0 + list_virt_reverse = 0 + + n=0 + do i = 1, mo_num + if (mo_class(i) == 'Virtual')then + n += 1 + list_virt(n) = i + list_virt_reverse(i) = n + endif + enddo + print *, 'Virtual MOs:' + print *, list_virt(1:n_virt_orb) + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_del , (dim_list_del_orb) ] +&BEGIN_PROVIDER [ integer, list_del_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are deleted. + END_DOC + integer :: i, n + list_del = 0 + list_del_reverse = 0 + + n=0 + do i = 1, mo_num + if (mo_class(i) == 'Deleted')then + n += 1 + list_del(n) = i + list_del_reverse(i) = n + endif + enddo + print *, 'Deleted MOs:' + print *, list_del(1:n_del_orb) + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_act , (dim_list_act_orb) ] +&BEGIN_PROVIDER [ integer, list_act_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are in the active. + END_DOC + integer :: i, n + list_act = 0 + list_act_reverse = 0 + + n=0 + do i = 1, mo_num + if (mo_class(i) == 'Active')then + n += 1 + list_act(n) = i + list_act_reverse(i) = n + endif + enddo + print *, 'Active MOs:' + print *, list_act(1:n_act_orb) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer, list_core_inact , (n_core_inact_orb) ] +&BEGIN_PROVIDER [ integer, list_core_inact_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of indices of the core and inactive MOs + END_DOC + integer :: i,itmp + call bitstring_to_list(reunion_of_core_inact_bitmask(1,1), list_core_inact, itmp, N_int) + list_core_inact_reverse = 0 + ASSERT (itmp == n_core_inact_orb) + do i = 1, n_core_inact_orb + list_core_inact_reverse(list_core_inact(i)) = i + enddo + print *, 'Core and Inactive MOs:' + print *, list_core_inact(1:n_core_inact_orb) +END_PROVIDER + + + BEGIN_PROVIDER [ integer, list_core_inact_act , (n_core_inact_act_orb) ] +&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of indices of the core inactive and active MOs + END_DOC + integer :: i,itmp + call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), list_core_inact_act, itmp, N_int) + list_core_inact_act_reverse = 0 + ASSERT (itmp == n_core_inact_act_orb) + do i = 1, n_core_inact_act_orb + list_core_inact_act_reverse(list_core_inact_act(i)) = i + enddo + print *, 'Core, Inactive and Active MOs:' + print *, list_core_inact_act(1:n_core_inact_act_orb) +END_PROVIDER + + + BEGIN_PROVIDER [ integer, list_inact_act , (n_inact_act_orb) ] +&BEGIN_PROVIDER [ integer, list_inact_act_reverse, (mo_num) ] + implicit none + BEGIN_DOC + ! List of indices of the inactive and active MOs + END_DOC + integer :: i,itmp + call bitstring_to_list(reunion_of_inact_act_bitmask(1,1), list_inact_act, itmp, N_int) + list_inact_act_reverse = 0 + ASSERT (itmp == n_inact_act_orb) + do i = 1, n_inact_act_orb + list_inact_act_reverse(list_inact_act(i)) = i + enddo + print *, 'Inactive and Active MOs:' + print *, list_inact_act(1:n_inact_act_orb) +END_PROVIDER + diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f index 74351760..daf3f68b 100644 --- a/src/casscf/bielec.irp.f +++ b/src/casscf/bielec.irp.f @@ -1,65 +1,58 @@ - BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)] +BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] BEGIN_DOC ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active ! indices are unshifted orbital numbers END_DOC implicit none integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 - double precision, allocatable :: integrals_array(:,:) real*8 :: mo_two_e_integral - allocate(integrals_array(mo_num,mo_num)) + bielec_PQxx(:,:,:,:) = 0.d0 + PROVIDE mo_two_e_integrals_in_map - bielec_PQxx = 0.d0 - - do i=1,n_core_orb - ii=list_core(i) - do j=i,n_core_orb - jj=list_core(j) - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num - bielec_PQxx(p,q,i,j)=integrals_array(p,q) - bielec_PQxx(p,q,j,i)=integrals_array(p,q) - end do - end do + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,ii,j,jj,i3,j3) & + !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx, & + !$OMP n_act_orb,mo_integrals_map,list_act) + + !$OMP DO + do i=1,n_core_inact_orb + ii=list_core_inact(i) + do j=i,n_core_inact_orb + jj=list_core_inact(j) + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j),mo_integrals_map) + bielec_PQxx(:,:,j,i)=bielec_PQxx(:,:,i,j) end do do j=1,n_act_orb jj=list_act(j) - j3=j+n_core_orb - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num - bielec_PQxx(p,q,i,j3)=integrals_array(p,q) - bielec_PQxx(p,q,j3,i)=integrals_array(p,q) - end do - end do + j3=j+n_core_inact_orb + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j3),mo_integrals_map) + bielec_PQxx(:,:,j3,i)=bielec_PQxx(:,:,i,j3) end do end do + !$OMP END DO - ! (ij|pq) + !$OMP DO do i=1,n_act_orb ii=list_act(i) - i3=i+n_core_orb + i3=i+n_core_inact_orb do j=i,n_act_orb jj=list_act(j) - j3=j+n_core_orb - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num - bielec_PQxx(p,q,i3,j3)=integrals_array(p,q) - bielec_PQxx(p,q,j3,i3)=integrals_array(p,q) - end do - end do + j3=j+n_core_inact_orb + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i3,j3),mo_integrals_map) + bielec_PQxx(:,:,j3,i3)=bielec_PQxx(:,:,i3,j3) end do end do + !$OMP END DO + + !$OMP END PARALLEL END_PROVIDER -BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)] +BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_orb+n_act_orb,n_core_inact_orb+n_act_orb, mo_num)] BEGIN_DOC ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active ! indices are unshifted orbital numbers @@ -69,17 +62,24 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_a double precision, allocatable :: integrals_array(:,:) real*8 :: mo_two_e_integral + PROVIDE mo_two_e_integrals_in_map + bielec_PxxQ = 0.d0 + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,ii,j,jj,i3,j3,integrals_array) & + !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ, & + !$OMP n_act_orb,mo_integrals_map,list_act) + allocate(integrals_array(mo_num,mo_num)) - bielec_PxxQ = 0.d0 - - do i=1,n_core_orb - ii=list_core(i) - do j=i,n_core_orb - jj=list_core(j) - call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num + !$OMP DO + do i=1,n_core_inact_orb + ii=list_core_inact(i) + do j=i,n_core_inact_orb + jj=list_core_inact(j) + call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) + do q=1,mo_num + do p=1,mo_num bielec_PxxQ(p,i,j,q)=integrals_array(p,q) bielec_PxxQ(p,j,i,q)=integrals_array(q,p) end do @@ -87,34 +87,41 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_a end do do j=1,n_act_orb jj=list_act(j) - j3=j+n_core_orb - call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num + j3=j+n_core_inact_orb + call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) + do q=1,mo_num + do p=1,mo_num bielec_PxxQ(p,i,j3,q)=integrals_array(p,q) bielec_PxxQ(p,j3,i,q)=integrals_array(q,p) end do end do end do end do + !$OMP END DO ! (ip|qj) + !$OMP DO do i=1,n_act_orb ii=list_act(i) - i3=i+n_core_orb + i3=i+n_core_inact_orb do j=i,n_act_orb jj=list_act(j) - j3=j+n_core_orb - call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num + j3=j+n_core_inact_orb + call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) + do q=1,mo_num + do p=1,mo_num bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q) bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p) end do end do end do end do + !$OMP END DO + + deallocate(integrals_array) + !$OMP END PARALLEL + END_PROVIDER @@ -125,24 +132,25 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] END_DOC implicit none integer :: i,j,k,p,t,u,v - double precision, allocatable :: integrals_array(:) - real*8 :: mo_two_e_integral + double precision, external :: mo_two_e_integral + PROVIDE mo_two_e_integrals_in_map - allocate(integrals_array(mo_num)) - - do i=1,n_act_orb - t=list_act(i) + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP PRIVATE(i,j,k,p,t,u,v) & + !$OMP SHARED(mo_num,n_act_orb,list_act,bielecCI) + do p=1,mo_num do j=1,n_act_orb u=list_act(j) do k=1,n_act_orb v=list_act(k) - ! (tu|vp) - call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array,mo_integrals_map) - do p=1,mo_num - bielecCI(i,k,j,p)=integrals_array(p) + do i=1,n_act_orb + t=list_act(i) + bielecCI(i,k,j,p) = mo_two_e_integral(t,u,v,p) end do end do end do end do -END_PROVIDER + !$OMP END PARALLEL DO +END_PROVIDER + diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f index ca1c8e9d..cb09be3e 100644 --- a/src/casscf/bielec_natorb.irp.f +++ b/src/casscf/bielec_natorb.irp.f @@ -1,180 +1,264 @@ - BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)] + BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] BEGIN_DOC ! integral (pq|xx) in the basis of natural MOs ! indices are unshifted orbital numbers END_DOC implicit none integer :: i,j,k,l,t,u,p,q,pp - real*8 :: d(n_act_orb) + double precision, allocatable :: f(:,:,:), d(:,:,:) - bielec_PQxx_no(:,:,:,:) = bielec_PQxx(:,:,:,:) - do j=1,mo_num - do k=1,n_core_orb+n_act_orb - do l=1,n_core_orb+n_act_orb + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,p,pp,d,f) & + !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & + !$OMP bielec_PQxx_no,bielec_PQxx,list_act,natorbsCI) + + allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), & + d(n_act_orb,mo_num,n_core_inact_act_orb)) + + !$OMP DO + do l=1,n_core_inact_act_orb + bielec_PQxx_no(:,:,:,l) = bielec_PQxx(:,:,:,l) + + do k=1,n_core_inact_act_orb + do j=1,mo_num do p=1,n_act_orb - d(p)=0.D0 + f(p,j,k)=bielec_PQxx_no(list_act(p),j,k,l) end do + end do + end do + call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0, & + natorbsCI, size(natorbsCI,1), & + f, n_act_orb, & + 0.d0, & + d, n_act_orb) + do k=1,n_core_inact_act_orb + do j=1,mo_num do p=1,n_act_orb pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PQxx_no(list_act(q),j,k,l)*natorbsCI(q,p) - end do + bielec_PQxx_no(list_act(p),j,k,l)=d(pp,j,k) end do + end do + + do j=1,mo_num do p=1,n_act_orb - bielec_PQxx_no(list_act(p),j,k,l)=d(p) + f(p,j,k)=bielec_PQxx_no(j,list_act(p),k,l) + end do + end do + end do + call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0, & + natorbsCI, n_act_orb, & + f, n_act_orb, & + 0.d0, & + d, n_act_orb) + do k=1,n_core_inact_act_orb + do p=1,n_act_orb + pp=n_act_orb-p+1 + do j=1,mo_num + bielec_PQxx_no(j,list_act(p),k,l)=d(pp,j,k) end do end do end do end do - ! 2nd quarter - do j=1,mo_num - do k=1,n_core_orb+n_act_orb - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d(p)=0.D0 + !$OMP END DO NOWAIT + + deallocate (f,d) + + allocate (f(mo_num,mo_num,n_act_orb),d(mo_num,mo_num,n_act_orb)) + + !$OMP DO + do l=1,n_core_inact_act_orb + + do p=1,n_act_orb + do k=1,mo_num + do j=1,mo_num + f(j,k,p) = bielec_PQxx_no(j,k,n_core_inact_orb+p,l) end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PQxx_no(j,list_act(q),k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PQxx_no(j,list_act(p),k,l)=d(p) + end do + end do + call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0, & + f, mo_num*mo_num, & + natorbsCI, n_act_orb, & + 0.d0, & + d, mo_num*mo_num) + do p=1,n_act_orb + pp=n_act_orb-p+1 + do k=1,mo_num + do j=1,mo_num + bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,pp) end do end do end do end do - ! 3rd quarter - do j=1,mo_num - do k=1,mo_num - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d(p)=0.D0 + !$OMP END DO NOWAIT + + !$OMP BARRIER + + !$OMP DO + do l=1,n_core_inact_act_orb + do p=1,n_act_orb + do k=1,mo_num + do j=1,mo_num + f(j,k,p) = bielec_PQxx_no(j,k,l,n_core_inact_orb+p) end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PQxx_no(j,k,n_core_orb+q,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PQxx_no(j,k,n_core_orb+p,l)=d(p) - end do - end do - end do - end do - ! 4th quarter - do j=1,mo_num - do k=1,mo_num - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PQxx_no(j,k,l,n_core_orb+q)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PQxx_no(j,k,l,n_core_orb+p)=d(p) + end do + end do + call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0, & + f, mo_num*mo_num, & + natorbsCI, n_act_orb, & + 0.d0, & + d, mo_num*mo_num) + do p=1,n_act_orb + pp=n_act_orb-p+1 + do k=1,mo_num + do j=1,mo_num + bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,pp) end do end do end do end do + !$OMP END DO + + deallocate (f,d) + !$OMP END PARALLEL END_PROVIDER -BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)] +BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] BEGIN_DOC ! integral (px|xq) in the basis of natural MOs ! indices are unshifted orbital numbers END_DOC implicit none integer :: i,j,k,l,t,u,p,q,pp - real*8 :: d(n_act_orb) + double precision, allocatable :: f(:,:,:), d(:,:,:) - bielec_PxxQ_no(:,:,:,:) = bielec_PxxQ(:,:,:,:) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,p,pp,d,f) & + !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & + !$OMP bielec_PxxQ_no,bielec_PxxQ,list_act,natorbsCI) + + allocate (f(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb), & + d(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb)) + + !$OMP DO do j=1,mo_num - do k=1,n_core_orb+n_act_orb - do l=1,n_core_orb+n_act_orb + bielec_PxxQ_no(:,:,:,j) = bielec_PxxQ(:,:,:,j) + do l=1,n_core_inact_act_orb + do k=1,n_core_inact_act_orb do p=1,n_act_orb - d(p)=0.D0 + f(p,k,l) = bielec_PxxQ_no(list_act(p),k,l,j) end do + end do + end do + call dgemm('T','N',n_act_orb,n_core_inact_act_orb**2,n_act_orb,1.d0, & + natorbsCI, size(natorbsCI,1), & + f, n_act_orb, & + 0.d0, & + d, n_act_orb) + do l=1,n_core_inact_act_orb + do k=1,n_core_inact_act_orb do p=1,n_act_orb pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PxxQ_no(list_act(q),k,l,j)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PxxQ_no(list_act(p),k,l,j)=d(p) + bielec_PxxQ_no(list_act(p),k,l,j)=d(pp,k,l) end do end do end do end do - ! 2nd quarter - do j=1,mo_num - do k=1,n_core_orb+n_act_orb - do l=1,n_core_orb+n_act_orb + !$OMP END DO NOWAIT + + deallocate (f,d) + + allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), & + d(n_act_orb,mo_num,n_core_inact_act_orb)) + + !$OMP DO + do k=1,mo_num + do l=1,n_core_inact_act_orb + do j=1,mo_num do p=1,n_act_orb - d(p)=0.D0 + f(p,j,l) = bielec_PxxQ_no(j,n_core_inact_orb+p,l,k) end do + end do + end do + call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0, & + natorbsCI, size(natorbsCI,1), & + f, n_act_orb, & + 0.d0, & + d, n_act_orb) + do l=1,n_core_inact_act_orb + do j=1,mo_num do p=1,n_act_orb pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PxxQ_no(j,k,l,list_act(p))=d(p) + bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(pp,j,l) end do end do end do end do - ! 3rd quarter - do j=1,mo_num - do k=1,mo_num - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d(p)=0.D0 + !$OMP END DO NOWAIT + + deallocate(f,d) + + allocate(f(mo_num,n_core_inact_act_orb,n_act_orb), & + d(mo_num,n_core_inact_act_orb,n_act_orb) ) + + !$OMP DO + do k=1,mo_num + do p=1,n_act_orb + do l=1,n_core_inact_act_orb + do j=1,mo_num + f(j,l,p) = bielec_PxxQ_no(j,l,n_core_inact_orb+p,k) end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PxxQ_no(j,n_core_orb+q,l,k)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PxxQ_no(j,n_core_orb+p,l,k)=d(p) + end do + end do + call dgemm('N','N',mo_num*n_core_inact_act_orb,n_act_orb,n_act_orb,1.d0, & + f, mo_num*n_core_inact_act_orb, & + natorbsCI, size(natorbsCI,1), & + 0.d0, & + d, mo_num*n_core_inact_act_orb) + do p=1,n_act_orb + pp=n_act_orb-p+1 + do l=1,n_core_inact_act_orb + do j=1,mo_num + bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(j,l,pp) end do end do end do end do - ! 4th quarter - do j=1,mo_num - do k=1,mo_num - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d(p)=0.D0 + !$OMP END DO NOWAIT + + !$OMP BARRIER + + !$OMP DO + do l=1,n_core_inact_act_orb + do p=1,n_act_orb + do k=1,n_core_inact_act_orb + do j=1,mo_num + f(j,k,p) = bielec_PxxQ_no(j,k,l,n_core_inact_orb+p) end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielec_PxxQ_no(j,l,n_core_orb+q,k)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PxxQ_no(j,l,n_core_orb+p,k)=d(p) + end do + end do + call dgemm('N','N',mo_num*n_core_inact_act_orb,n_act_orb,n_act_orb,1.d0, & + f, mo_num*n_core_inact_act_orb, & + natorbsCI, size(natorbsCI,1), & + 0.d0, & + d, mo_num*n_core_inact_act_orb) + do p=1,n_act_orb + pp=n_act_orb-p+1 + do k=1,n_core_inact_act_orb + do j=1,mo_num + bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)=d(j,k,pp) end do end do end do end do + !$OMP END DO NOWAIT + deallocate(f,d) + !$OMP END PARALLEL END_PROVIDER @@ -186,85 +270,112 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)] END_DOC implicit none integer :: i,j,k,l,t,u,p,q,pp - real*8 :: d(n_act_orb) + double precision, allocatable :: f(:,:,:), d(:,:,:) - bielecCI_no(:,:,:,:) = bielecCI(:,:,:,:) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,p,pp,d,f) & + !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & + !$OMP bielecCI_no,bielecCI,list_act,natorbsCI) - do j=1,n_act_orb + allocate (f(n_act_orb,n_act_orb,mo_num), & + d(n_act_orb,n_act_orb,mo_num)) + + !$OMP DO + do l=1,mo_num + bielecCI_no(:,:,:,l) = bielecCI(:,:,:,l) do k=1,n_act_orb - do l=1,mo_num + do j=1,n_act_orb do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielecCI_no(p,j,k,l)=d(p) + f(p,j,k)=bielecCI_no(p,j,k,l) end do end do end do - end do - ! 2nd quarter - do j=1,n_act_orb + call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0, & + natorbsCI, size(natorbsCI,1), & + f, n_act_orb, & + 0.d0, & + d, n_act_orb) do k=1,n_act_orb - do l=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do + do j=1,n_act_orb do p=1,n_act_orb pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p) - end do + bielecCI_no(p,j,k,l)=d(pp,j,k) end do + end do + + do j=1,n_act_orb do p=1,n_act_orb - bielecCI_no(j,p,k,l)=d(p) + f(p,j,k)=bielecCI_no(j,p,k,l) end do end do end do - end do - ! 3rd quarter - do j=1,n_act_orb + call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0, & + natorbsCI, n_act_orb, & + f, n_act_orb, & + 0.d0, & + d, n_act_orb) do k=1,n_act_orb - do l=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 + do p=1,n_act_orb + pp=n_act_orb-p+1 + do j=1,n_act_orb + bielecCI_no(j,p,k,l)=d(pp,j,k) end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p) - end do + end do + end do + + do p=1,n_act_orb + do k=1,n_act_orb + do j=1,n_act_orb + f(j,k,p)=bielecCI_no(j,k,p,l) end do - do p=1,n_act_orb - bielecCI_no(j,k,p,l)=d(p) + end do + end do + call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0, & + f, n_act_orb*n_act_orb, & + natorbsCI, n_act_orb, & + 0.d0, & + d, n_act_orb*n_act_orb) + + do p=1,n_act_orb + pp=n_act_orb-p+1 + do k=1,n_act_orb + do j=1,n_act_orb + bielecCI_no(j,k,p,l)=d(j,k,pp) end do end do end do end do - ! 4th quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 + !$OMP END DO + + !$OMP DO + do l=1,n_act_orb + do p=1,n_act_orb + do k=1,n_act_orb + do j=1,n_act_orb + f(j,k,p)=bielecCI_no(j,k,l,list_act(p)) end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielecCI_no(j,k,l,list_act(p))=d(p) + end do + end do + call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0, & + f, n_act_orb*n_act_orb, & + natorbsCI, n_act_orb, & + 0.d0, & + d, n_act_orb*n_act_orb) + + do p=1,n_act_orb + pp=n_act_orb-p+1 + do k=1,n_act_orb + do j=1,n_act_orb + bielecCI_no(j,k,l,list_act(p))=d(j,k,pp) end do end do end do end do + !$OMP END DO + + deallocate(d,f) + !$OMP END PARALLEL + END_PROVIDER diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f index b3c8988f..f00bc7c8 100644 --- a/src/casscf/gradient.irp.f +++ b/src/casscf/gradient.irp.f @@ -5,7 +5,7 @@ BEGIN_PROVIDER [ integer, nMonoEx ] ! Number of single excitations END_DOC implicit none - nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb + nMonoEx=n_core_inact_orb*n_act_orb+n_core_inact_orb*n_virt_orb+n_act_orb*n_virt_orb END_PROVIDER BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] @@ -17,8 +17,8 @@ END_PROVIDER implicit none integer :: i,t,a,ii,tt,aa,indx indx=0 - do ii=1,n_core_orb - i=list_core(ii) + do ii=1,n_core_inact_orb + i=list_core_inact(ii) do tt=1,n_act_orb t=list_act(tt) indx+=1 @@ -28,8 +28,8 @@ END_PROVIDER end do end do - do ii=1,n_core_orb - i=list_core(ii) + do ii=1,n_core_inact_orb + i=list_core_inact(ii) do aa=1,n_virt_orb a=list_virt(aa) indx+=1 @@ -145,14 +145,14 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] real*8 :: norm_grad indx=0 - do i=1,n_core_orb + do i=1,n_core_inact_orb do t=1,n_act_orb indx+=1 gradvec2(indx)=gradvec_it(i,t) end do end do - do i=1,n_core_orb + do i=1,n_core_inact_orb do a=1,n_virt_orb indx+=1 gradvec2(indx)=gradvec_ia(i,a) @@ -181,7 +181,7 @@ END_PROVIDER real*8 function gradvec_it(i,t) BEGIN_DOC - ! the orbital gradient core -> active + ! the orbital gradient core/inactive -> active ! we assume natural orbitals END_DOC implicit none @@ -190,16 +190,16 @@ real*8 function gradvec_it(i,t) integer :: ii,tt,v,vv,x,y integer :: x3,y3 - ii=list_core(i) + ii=list_core_inact(i) tt=list_act(t) gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii)) gradvec_it-=occnum(tt)*Fipq(ii,tt) do v=1,n_act_orb vv=list_act(v) do x=1,n_act_orb - x3=x+n_core_orb + x3=x+n_core_inact_orb do y=1,n_act_orb - y3=y+n_core_orb + y3=y+n_core_inact_orb gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3) end do end do @@ -209,12 +209,12 @@ end function gradvec_it real*8 function gradvec_ia(i,a) BEGIN_DOC - ! the orbital gradient core -> virtual + ! the orbital gradient core/inactive -> virtual END_DOC implicit none integer :: i,a,ii,aa - ii=list_core(i) + ii=list_core_inact(i) aa=list_virt(a) gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii)) gradvec_ia*=2.D0 diff --git a/src/casscf/hessian.irp.f b/src/casscf/hessian.irp.f index e047c5fd..75a27410 100644 --- a/src/casscf/hessian.irp.f +++ b/src/casscf/hessian.irp.f @@ -204,10 +204,10 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] endif indx=1 - do i=1,n_core_orb + do i=1,n_core_inact_orb do t=1,n_act_orb jndx=indx - do j=i,n_core_orb + do j=i,n_core_inact_orb if (i.eq.j) then ustart=t else @@ -219,7 +219,7 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] jndx+=1 end do end do - do j=1,n_core_orb + do j=1,n_core_inact_orb do a=1,n_virt_orb hessmat2(indx,jndx)=hessmat_itja(i,t,j,a) hessmat2(jndx,indx)=hessmat2(indx,jndx) @@ -237,10 +237,10 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] end do end do - do i=1,n_core_orb + do i=1,n_core_inact_orb do a=1,n_virt_orb jndx=indx - do j=i,n_core_orb + do j=i,n_core_inact_orb if (i.eq.j) then bstart=a else @@ -286,7 +286,7 @@ END_PROVIDER real*8 function hessmat_itju(i,t,j,u) BEGIN_DOC - ! the orbital hessian for core->act,core->act + ! the orbital hessian for core/inactive -> active, core/inactive -> active ! i, t, j, u are list indices, the corresponding orbitals are ii,tt,jj,uu ! ! we assume natural orbitals @@ -295,7 +295,7 @@ real*8 function hessmat_itju(i,t,j,u) integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj real*8 :: term,t2 - ii=list_core(i) + ii=list_core_inact(i) tt=list_act(t) if (i.eq.j) then if (t.eq.u) then @@ -343,7 +343,7 @@ real*8 function hessmat_itju(i,t,j,u) end if else ! it/ju - jj=list_core(j) + jj=list_core_inact(j) uu=list_act(u) if (t.eq.u) then term=occnum(tt)*Fipq(ii,jj) @@ -374,16 +374,16 @@ end function hessmat_itju real*8 function hessmat_itja(i,t,j,a) BEGIN_DOC - ! the orbital hessian for core->act,core->virt + ! the orbital hessian for core/inactive -> active, core/inactive -> virtual END_DOC implicit none integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y real*8 :: term ! it/ja - ii=list_core(i) + ii=list_core_inact(i) tt=list_act(t) - jj=list_core(j) + jj=list_core_inact(j) aa=list_virt(a) term=2.D0*(4.D0*bielec_pxxq_no(aa,j,i,tt) & -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt)) @@ -407,17 +407,17 @@ end function hessmat_itja real*8 function hessmat_itua(i,t,u,a) BEGIN_DOC - ! the orbital hessian for core->act,act->virt + ! the orbital hessian for core/inactive -> active, active -> virtual END_DOC implicit none integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3 real*8 :: term - ii=list_core(i) + ii=list_core_inact(i) tt=list_act(t) - t3=t+n_core_orb + t3=t+n_core_inact_orb uu=list_act(u) - u3=u+n_core_orb + u3=u+n_core_inact_orb aa=list_virt(a) if (t.eq.u) then term=-occnum(tt)*Fipq(aa,ii) @@ -428,11 +428,11 @@ real*8 function hessmat_itua(i,t,u,a) +bielec_pxxq_no(aa,t3,u3,ii)) do v=1,n_act_orb vv=list_act(v) - v3=v+n_core_orb + v3=v+n_core_inact_orb do x=1,n_act_orb integer :: x3 xx=list_act(x) - x3=x+n_core_orb + x3=x+n_core_inact_orb term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,ii,v3,x3) & +(P0tuvx_no(t,v,u,x)+P0tuvx_no(t,v,x,u)) & *bielec_pqxx_no(aa,xx,v3,i)) @@ -448,13 +448,13 @@ end function hessmat_itua real*8 function hessmat_iajb(i,a,j,b) BEGIN_DOC - ! the orbital hessian for core->virt,core->virt + ! the orbital hessian for core/inactive -> virtual, core/inactive -> virtual END_DOC implicit none integer :: i,a,j,b,ii,aa,jj,bb real*8 :: term - ii=list_core(i) + ii=list_core_inact(i) aa=list_virt(a) if (i.eq.j) then if (a.eq.b) then @@ -469,7 +469,7 @@ real*8 function hessmat_iajb(i,a,j,b) end if else ! ia/jb - jj=list_core(j) + jj=list_core_inact(j) bb=list_virt(b) term=2.D0*(4.D0*bielec_pxxq_no(aa,i,j,bb)-bielec_pqxx_no(aa,bb,i,j) & -bielec_pxxq_no(aa,j,i,bb)) @@ -484,17 +484,17 @@ end function hessmat_iajb real*8 function hessmat_iatb(i,a,t,b) BEGIN_DOC - ! the orbital hessian for core->virt,act->virt + ! the orbital hessian for core/inactive -> virtual, active -> virtual END_DOC implicit none integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3 real*8 :: term - ii=list_core(i) + ii=list_core_inact(i) aa=list_virt(a) tt=list_act(t) bb=list_virt(b) - t3=t+n_core_orb + t3=t+n_core_inact_orb term=occnum(tt)*(4.D0*bielec_pxxq_no(aa,i,t3,bb)-bielec_pxxq_no(aa,t3,i,bb)& -bielec_pqxx_no(aa,bb,i,t3)) if (a.eq.b) then @@ -533,10 +533,10 @@ real*8 function hessmat_taub(t,a,u,b) t1-=occnum(tt)*Fipq(tt,tt) do v=1,n_act_orb vv=list_act(v) - v3=v+n_core_orb + v3=v+n_core_inact_orb do x=1,n_act_orb xx=list_act(x) - x3=x+n_core_orb + x3=x+n_core_inact_orb t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) & +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & bielec_pxxq_no(aa,x3,v3,aa)) @@ -552,10 +552,10 @@ real*8 function hessmat_taub(t,a,u,b) term=occnum(tt)*Fipq(aa,bb) do v=1,n_act_orb vv=list_act(v) - v3=v+n_core_orb + v3=v+n_core_inact_orb do x=1,n_act_orb xx=list_act(x) - x3=x+n_core_orb + x3=x+n_core_inact_orb term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & *bielec_pxxq_no(aa,x3,v3,bb)) @@ -569,10 +569,10 @@ real*8 function hessmat_taub(t,a,u,b) term=0.D0 do v=1,n_act_orb vv=list_act(v) - v3=v+n_core_orb + v3=v+n_core_inact_orb do x=1,n_act_orb xx=list_act(x) - x3=x+n_core_orb + x3=x+n_core_inact_orb term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) & *bielec_pxxq_no(aa,x3,v3,bb)) @@ -606,14 +606,14 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub indx=0 - do i=1,n_core_orb + do i=1,n_core_inact_orb do t=1,n_act_orb indx+=1 hessdiag(indx)=hessmat_itju(i,t,i,t) end do end do - do i=1,n_core_orb + do i=1,n_core_inact_orb do a=1,n_virt_orb indx+=1 hessdiag(indx)=hessmat_iajb(i,a,i,a) diff --git a/src/casscf/mcscf_fock.irp.f b/src/casscf/mcscf_fock.irp.f index 84b87248..e4568405 100644 --- a/src/casscf/mcscf_fock.irp.f +++ b/src/casscf/mcscf_fock.irp.f @@ -12,8 +12,8 @@ BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] end do ! the inactive Fock matrix - do k=1,n_core_orb - kk=list_core(k) + do k=1,n_core_inact_orb + kk=list_core_inact(k) do q=1,mo_num do p=1,mo_num Fipq(p,q)+=2.D0*bielec_pqxx_no(p,q,k,k) -bielec_pxxq_no(p,k,k,q) diff --git a/src/casscf/natorb.irp.f b/src/casscf/natorb.irp.f index 52cd3747..dcbcd413 100644 --- a/src/casscf/natorb.irp.f +++ b/src/casscf/natorb.irp.f @@ -6,8 +6,8 @@ integer :: i occnum=0.D0 - do i=1,n_core_orb - occnum(list_core(i))=2.D0 + do i=1,n_core_inact_orb + occnum(list_core_inact(i))=2.D0 end do do i=1,n_act_orb diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f index 16680452..0f9df016 100644 --- a/src/casscf/neworbs.irp.f +++ b/src/casscf/neworbs.irp.f @@ -125,8 +125,8 @@ BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ] ! the orbital rotation matrix T Tmat(:,:)=0.D0 indx=1 - do i=1,n_core_orb - ii=list_core(i) + do i=1,n_core_inact_orb + ii=list_core_inact(i) do t=1,n_act_orb tt=list_act(t) indx+=1 @@ -134,8 +134,8 @@ BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ] Tmat(tt,ii)=-SXvector(indx) end do end do - do i=1,n_core_orb - ii=list_core(i) + do i=1,n_core_inact_orb + ii=list_core_inact(i) do a=1,n_virt_orb aa=list_virt(a) indx+=1 diff --git a/src/casscf/tot_en.irp.f b/src/casscf/tot_en.irp.f index ce787232..1d70e087 100644 --- a/src/casscf/tot_en.irp.f +++ b/src/casscf/tot_en.irp.f @@ -10,19 +10,19 @@ real*8 :: e_one_all,e_two_all e_one_all=0.D0 e_two_all=0.D0 - do i=1,n_core_orb - ii=list_core(i) + do i=1,n_core_inact_orb + ii=list_core_inact(i) e_one_all+=2.D0*mo_one_e_integrals(ii,ii) - do j=1,n_core_orb - jj=list_core(j) + do j=1,n_core_inact_orb + jj=list_core_inact(j) e_two_all+=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i) end do do t=1,n_act_orb tt=list_act(t) - t3=t+n_core_orb + t3=t+n_core_inact_orb do u=1,n_act_orb uu=list_act(u) - u3=u+n_core_orb + u3=u+n_core_inact_orb e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) & -bielec_PQxx(tt,ii,i,u3)) end do @@ -34,9 +34,9 @@ uu=list_act(u) e_one_all+=D0tu(t,u)*mo_one_e_integrals(tt,uu) do v=1,n_act_orb - v3=v+n_core_orb + v3=v+n_core_inact_orb do x=1,n_act_orb - x3=x+n_core_orb + x3=x+n_core_inact_orb e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxx(tt,uu,v3,x3) end do end do @@ -44,12 +44,12 @@ end do ecore =nuclear_repulsion ecore_bis=nuclear_repulsion - do i=1,n_core_orb - ii=list_core(i) + do i=1,n_core_inact_orb + ii=list_core_inact(i) ecore +=2.D0*mo_one_e_integrals(ii,ii) ecore_bis+=2.D0*mo_one_e_integrals(ii,ii) - do j=1,n_core_orb - jj=list_core(j) + do j=1,n_core_inact_orb + jj=list_core_inact(j) ecore +=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i) ecore_bis+=2.D0*bielec_PxxQ(ii,i,j,jj)-bielec_PxxQ(ii,j,j,ii) end do @@ -61,14 +61,14 @@ etwo_ter=0.D0 do t=1,n_act_orb tt=list_act(t) - t3=t+n_core_orb + t3=t+n_core_inact_orb do u=1,n_act_orb uu=list_act(u) - u3=u+n_core_orb + u3=u+n_core_inact_orb eone +=D0tu(t,u)*mo_one_e_integrals(tt,uu) eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu) - do i=1,n_core_orb - ii=list_core(i) + do i=1,n_core_inact_orb + ii=list_core_inact(i) eone +=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) & -bielec_PQxx(tt,ii,i,u3)) eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQ(tt,u3,i,ii) & @@ -76,10 +76,10 @@ end do do v=1,n_act_orb vv=list_act(v) - v3=v+n_core_orb + v3=v+n_core_inact_orb do x=1,n_act_orb xx=list_act(x) - x3=x+n_core_orb + x3=x+n_core_inact_orb real*8 :: h1,h2,h3 h1=bielec_PQxx(tt,uu,v3,x3) h2=bielec_PxxQ(tt,u3,v3,xx) diff --git a/src/selectors_full/selectors.irp.f b/src/selectors_full/selectors.irp.f index 4e14d65a..0531f731 100644 --- a/src/selectors_full/selectors.irp.f +++ b/src/selectors_full/selectors.irp.f @@ -38,35 +38,18 @@ END_PROVIDER END_DOC integer :: i,k -! if (threshold_selectors == 1.d0) then -! -! do i=1,N_det_selectors -! do k=1,N_int -! psi_selectors(k,1,i) = psi_det(k,1,i) -! psi_selectors(k,2,i) = psi_det(k,2,i) -! enddo -! enddo -! do k=1,N_states -! do i=1,N_det_selectors -! psi_selectors_coef(i,k) = psi_coef(i,k) -! enddo -! enddo -! -! else - + do i=1,N_det_selectors + do k=1,N_int + psi_selectors(k,1,i) = psi_det_sorted(k,1,i) + psi_selectors(k,2,i) = psi_det_sorted(k,2,i) + enddo + enddo + do k=1,N_states do i=1,N_det_selectors - do k=1,N_int - psi_selectors(k,1,i) = psi_det_sorted(k,1,i) - psi_selectors(k,2,i) = psi_det_sorted(k,2,i) - enddo - enddo - do k=1,N_states - do i=1,N_det_selectors - psi_selectors_coef(i,k) = psi_coef_sorted(i,k) - enddo + psi_selectors_coef(i,k) = psi_coef_sorted(i,k) enddo + enddo -! endif END_PROVIDER