diff --git a/plugins/CASSCF/EZFIO.cfg b/plugins/CASSCF/EZFIO.cfg new file mode 100644 index 00000000..e9e6e92e --- /dev/null +++ b/plugins/CASSCF/EZFIO.cfg @@ -0,0 +1,10 @@ +[energy] +type: double precision +doc: "Calculated CAS-SCF energy" +interface: ezfio + +[energy_pt2] +type: double precision +doc: "Calculated selected CAS-SCF energy with PT2 correction" +interface: ezfio + diff --git a/plugins/CASSCF/H_apply.irp.f b/plugins/CASSCF/H_apply.irp.f new file mode 100644 index 00000000..35c45fb6 --- /dev/null +++ b/plugins/CASSCF/H_apply.irp.f @@ -0,0 +1,39 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("CAS_SD") +print s + +s = H_apply("CAS_SD_selected_no_skip") +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +print s + +s = H_apply("CAS_SD_selected") +s.set_selection_pt2("epstein_nesbet_2x2") +print s + +s = H_apply("CAS_SD_PT2") +s.set_perturbation("epstein_nesbet_2x2") +print s + + +s = H_apply("CAS_S",do_double_exc=False) +print s + +s = H_apply("CAS_S_selected_no_skip",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +print s + +s = H_apply("CAS_S_selected",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +print s + +s = H_apply("CAS_S_PT2",do_double_exc=False) +s.set_perturbation("epstein_nesbet_2x2") +print s + +END_SHELL + diff --git a/plugins/CASSCF/NEEDED_CHILDREN_MODULES b/plugins/CASSCF/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..29e39f2f --- /dev/null +++ b/plugins/CASSCF/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Generators_CAS Perturbation Selectors_full diff --git a/plugins/CASSCF/README.rst b/plugins/CASSCF/README.rst new file mode 100644 index 00000000..ceeb7477 --- /dev/null +++ b/plugins/CASSCF/README.rst @@ -0,0 +1,20 @@ +====== +CASSCF +====== + +This module is not a "real" CAS-SCF. It is an orbital optimization step done by : + +1) Doing the CAS+SD +2) Taking one-electron density matrix +3) Cancelling all active-active rotations +4) Finding the order which matches with the input MOs + + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/CASSCF/casscf.irp.f b/plugins/CASSCF/casscf.irp.f new file mode 100644 index 00000000..4e7450dc --- /dev/null +++ b/plugins/CASSCF/casscf.irp.f @@ -0,0 +1,220 @@ +program casscf + implicit none + BEGIN_DOC +! Optimize MOs and CI coefficients of the CAS + END_DOC + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer(bit_kind), allocatable :: generators_bitmask_save(:,:,:,:) + + integer :: degree, N_generators_bitmask_save, N_det_ci + double precision :: E_old, E_CI + double precision :: selection_criterion_save, selection_criterion_min_save + + integer :: N_det_old + integer :: i, j, k, l + integer :: i_bit, j_bit, i_int, j_int + integer(bit_kind), allocatable :: bit_tmp(:), cas_bm(:) + character*(64) :: label + + allocate( pt2(N_states), norm_pert(N_states),H_pert_diag(N_states) ) + allocate( generators_bitmask_save(N_int,2,6,N_generators_bitmask) ) + allocate( bit_tmp(N_int), cas_bm(N_int) ) + + PROVIDE N_det_cas + N_det_old = 0 + pt2 = 1.d0 + E_CI = 1.d0 + E_old = 0.d0 + diag_algorithm = "Lapack" + selection_criterion_save = selection_criterion + selection_criterion_min_save = selection_criterion_min + + + cas_bm = 0_bit_kind + do i=1,N_cas_bitmask + do j=1,N_int + cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,1,i)) + cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,2,i)) + enddo + enddo + + ! Save CAS-SD bitmask + generators_bitmask_save = generators_bitmask + N_generators_bitmask_save = N_generators_bitmask + + ! Set the CAS bitmask + do i=1,6 + generators_bitmask(:,:,i,:) = cas_bitmask + enddo + N_generators_bitmask = N_cas_bitmask + SOFT_TOUCH generators_bitmask N_generators_bitmask + + + ! If the number of dets already in the file is larger than the requested + ! number of determinants, truncate the wf + if (N_det > N_det_max) then + call diagonalize_CI + call save_wavefunction + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = N_det_max + soft_touch N_det psi_det psi_coef + call diagonalize_CI + call save_wavefunction + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + endif + + ! Start MCSCF iteration + + ! CAS-CI + ! ------ + + E_old = E_CI + + ! Reset the selection criterion + selection_criterion = selection_criterion_save + selection_criterion_min = selection_criterion_min_save + SOFT_TOUCH selection_criterion_min selection_criterion selection_criterion_factor + + ! Set the CAS bitmask + do i=1,6 + generators_bitmask(:,:,i,:) = cas_bitmask + enddo + N_generators_bitmask = N_cas_bitmask + SOFT_TOUCH generators_bitmask N_generators_bitmask + + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_states))) > pt2_max) + N_det_old = N_det + call H_apply_CAS_SD_selected_no_skip(pt2, norm_pert, H_pert_diag, N_states) + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + if (N_det > N_det_max) then + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = N_det_max + soft_touch N_det psi_det psi_coef + endif + call diagonalize_CI + call save_wavefunction + print *, '======' + print *, 'CAS-CI' + print *, '======' + print *, '' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E(CAS) = ', CI_energy + print *, 'E(CAS)+PT2 = ', CI_energy+pt2 + print *, '-----' + print *, '' + E_CI = sum(CI_energy(1:N_states)+pt2(1:N_states))/dble(N_states) + + call ezfio_set_casscf_energy(CI_energy(1)) + if (abort_all) then + exit + endif + if (N_det == N_det_old) then + exit + endif + + enddo + + ! Super-CI + ! -------- + + selection_criterion_min = 1.d-12 + selection_criterion = 1.d-12 + + ! Set the CAS bitmask + generators_bitmask = generators_bitmask_save + N_generators_bitmask = N_generators_bitmask_save + SOFT_TOUCH generators_bitmask N_generators_bitmask selection_criterion selection_criterion_min selection_criterion_factor + + N_det_ci = N_det + + call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_states) + + + do i=1,mo_tot_num + i_int = ishft(i-1,-bit_kind_shift)+1 + i_bit = j-ishft(i_int-1,bit_kind_shift)-1 + bit_tmp(:) = 0_bit_kind + bit_tmp(i_int) = ibset(0_bit_kind,i_bit) + if (iand(bit_tmp(i_int), cas_bm(i_int)) == 0_bit_kind) then + ! Not a CAS MO + cycle + endif + + do j=1,mo_tot_num + if (j == i) then + cycle + endif + j_int = ishft(j-1,-bit_kind_shift)+1 + j_bit = j-ishft(j_int-1,bit_kind_shift)-1 + bit_tmp(:) = 0_bit_kind + bit_tmp(j_int) = ibset(0_bit_kind,j_bit) + if (iand(bit_tmp(j_int), cas_bm(j_int)) == 0_bit_kind) then + ! Not a CAS MO + cycle + endif + ! Now, both i and j are MOs of the CAS. De-couple them in the DM + one_body_dm_mo(i,j) = 0.d0 + enddo + + enddo + + SOFT_TOUCH one_body_dm_mo + + double precision :: mx, ov + double precision, allocatable :: mo_coef_old(:,:) + integer, allocatable :: iorder(:) + logical, allocatable :: selected(:) + allocate( mo_coef_old(size(mo_coef,1), size(mo_coef,2)), iorder(mo_tot_num), selected(mo_tot_num) ) + mo_coef_old = mo_coef + label = "Canonical" + call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,-1) + selected = .False. + do j=1,mo_tot_num + mx = -1.d0 + iorder(j) = j + do i=1,mo_tot_num + if (selected(i)) then + cycle + endif + ov = 0.d0 + do l=1,ao_num + do k=1,ao_num + ov = ov + mo_coef_old(k,j) * ao_overlap(k,l) * mo_coef(l,i) + enddo + enddo + ov= dabs(ov) + if (ov > mx) then + mx = ov + iorder(j) = i + endif + enddo + selected( iorder(j) ) = .True. + enddo + mo_coef_old = mo_coef + do i=1,mo_tot_num + mo_coef(:,i) = mo_coef_old(:,iorder(i)) + enddo + + call save_mos + + call write_double(6,E_CI,"Energy(CAS)") + + deallocate( mo_coef_old ) + deallocate( pt2, norm_pert,H_pert_diag ) + deallocate( generators_bitmask_save ) + deallocate( bit_tmp, cas_bm, iorder ) +end diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 044fa18b..29588369 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -262,13 +262,7 @@ END_PROVIDER logical :: exists integer :: j,i integer :: i_hole,i_part,i_gen - PROVIDE ezfio_filename -!do j = 1, N_int -! inact_bitmask(j,1) = xor(generators_bitmask(j,1,1,1),cas_bitmask(j,1,1)) -! inact_bitmask(j,2) = xor(generators_bitmask(j,2,1,1),cas_bitmask(j,2,1)) -! virt_bitmask(j,1) = xor(generators_bitmask(j,1,2,1),cas_bitmask(j,1,1)) -! virt_bitmask(j,2) = xor(generators_bitmask(j,2,2,1),cas_bitmask(j,2,1)) -!enddo + n_inact_orb = 0 n_virt_orb = 0 if(N_generators_bitmask == 1)then diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 3b737f5d..43c207d5 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -424,15 +424,6 @@ integer*8 function get_mo_map_size() get_mo_map_size = mo_integrals_map % n_elements end -subroutine clear_mo_map - implicit none - BEGIN_DOC - ! Frees the memory of the MO map - END_DOC - call map_deinit(mo_integrals_map) - FREE mo_integrals_map -end - BEGIN_TEMPLATE subroutine dump_$ao_integrals(filename) diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 1fa303b5..83f0ce05 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -503,3 +503,14 @@ BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_to enddo END_PROVIDER + + +subroutine clear_mo_map + implicit none + BEGIN_DOC + ! Frees the memory of the MO map + END_DOC + call map_deinit(mo_integrals_map) + FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti + FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map +end diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 3f4a260a..0d8ef5fa 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -98,7 +98,6 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) call write_time(output_mo_basis) mo_label = label - SOFT_TOUCH mo_coef mo_label end subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)