From b7992a11a95809003cdd5222fb0997200791b3ba Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 23 Oct 2019 00:11:55 +0200 Subject: [PATCH] working ... --- src/casscf/EZFIO.cfg | 5 + src/casscf/MORALITY | 1 + src/casscf/casscf.irp.f | 12 +- src/casscf/get_energy.irp.f | 50 +++--- src/casscf/neworbs.irp.f | 28 +++- src/casscf/superci_dm.irp.f | 207 ++++++++++++++++++++++++ src/casscf/swap_orb.irp.f | 307 +++++++++++++----------------------- src/mo_basis/utils.irp.f | 35 +++- 8 files changed, 417 insertions(+), 228 deletions(-) create mode 100644 src/casscf/MORALITY create mode 100644 src/casscf/superci_dm.irp.f diff --git a/src/casscf/EZFIO.cfg b/src/casscf/EZFIO.cfg index 48e296a5..529354fb 100644 --- a/src/casscf/EZFIO.cfg +++ b/src/casscf/EZFIO.cfg @@ -23,4 +23,9 @@ interface: ezfio,provider,ocaml default: False +[level_shift_casscf] +type: Positive_float +doc: Energy shift on the virtual MOs to improve SCF convergence +interface: ezfio,provider,ocaml +default: 0.05 diff --git a/src/casscf/MORALITY b/src/casscf/MORALITY new file mode 100644 index 00000000..9701a647 --- /dev/null +++ b/src/casscf/MORALITY @@ -0,0 +1 @@ +the CASCF can be obtained if a proper guess is given to the WF part diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index 8fe77fcc..dcf57835 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -6,19 +6,23 @@ program casscf no_vvvv_integrals = .True. pt2_max = 0.02 SOFT_TOUCH no_vvvv_integrals pt2_max + call run_stochastic_cipsi call run end subroutine run implicit none double precision :: energy_old, energy - logical :: converged + logical :: converged,state_following_casscf_save integer :: iteration converged = .False. energy = 0.d0 mo_label = "MCSCF" iteration = 1 + state_following_casscf_save = state_following_casscf + state_following_casscf = .True. + touch state_following_casscf do while (.not.converged) call run_stochastic_cipsi energy_old = energy @@ -35,12 +39,16 @@ subroutine run mo_coef = NewOrbs call save_mos iteration += 1 - N_det = N_det/2 + N_det = N_det/2 psi_det = psi_det_sorted psi_coef = psi_coef_sorted read_wf = .True. call clear_mo_map SOFT_TOUCH mo_coef N_det pt2_max psi_det psi_coef + if(iteration .gt. 3)then + state_following_casscf = state_following_casscf_save + touch state_following_casscf + endif enddo diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f index 4ac2b22e..3ff42fbc 100644 --- a/src/casscf/get_energy.irp.f +++ b/src/casscf/get_energy.irp.f @@ -16,30 +16,34 @@ subroutine routine_bis implicit none integer :: i,j double precision :: accu_d,accu_od - accu_d = 0.d0 - accu_od = 0.d0 - print*,'' - print*,'' - print*,'' - do i = 1, mo_num +!accu_d = 0.d0 +!accu_od = 0.d0 +!print*,'' +!print*,'' +!print*,'' +!do i = 1, mo_num ! write(*,'(100(F8.5,X))')super_ci_dm(i,:) - accu_d += super_ci_dm(i,i) - do j = i+1, mo_num - accu_od += dabs(super_ci_dm(i,j) - super_ci_dm(j,i)) - enddo - enddo - print*,'' - print*,'' - print*,'accu_d = ',accu_d - print*,'n_elec = ',elec_num - print*,'accu_od= ',accu_od - print*,'' - accu_d = 0.d0 - do i = 1, N_det - accu_d += psi_coef(i,1)**2 - enddo - print*,'accu_d = ',accu_d - provide superci_natorb +! accu_d += super_ci_dm(i,i) +! do j = i+1, mo_num +! accu_od += dabs(super_ci_dm(i,j) - super_ci_dm(j,i)) +! enddo +!enddo +!print*,'' +!print*,'' +!print*,'accu_d = ',accu_d +!print*,'n_elec = ',elec_num +!print*,'accu_od= ',accu_od +!print*,'' +!accu_d = 0.d0 +!do i = 1, N_det +! accu_d += psi_coef(i,1)**2 +!enddo +!print*,'accu_d = ',accu_d +!provide superci_natorb + + provide switch_mo_coef + mo_coef = switch_mo_coef + call save_mos end subroutine routine diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f index 980b4551..2af167be 100644 --- a/src/casscf/neworbs.irp.f +++ b/src/casscf/neworbs.irp.f @@ -24,6 +24,9 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] end do end do + do i = 1, nMonoEx + SXmatrix(i+1,i+1) += level_shift_casscf + enddo if (bavard) then do i=2,nMonoEx write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i) @@ -86,7 +89,7 @@ END_PROVIDER c0=SXeigenvec(1,best_vector_ovrlp_casscf) if (bavard) then write(6,*) ' SXdiag : eigenvalue for best overlap with ' - write(6,*) ' previous orbitals = ',SXeigenval(best_vector_ovrlp_casscf) + write(6,*) ' previous orbitals = ',SXeigenval(best_vector_ovrlp_casscf) write(6,*) ' weight of the 1st element ',c0 endif END_PROVIDER @@ -99,8 +102,10 @@ END_PROVIDER integer :: i double precision :: c0 c0=SXeigenvec(1,best_vector_ovrlp_casscf) + print*,'c0 = ',c0 do i=1,nMonoEx+1 SXvector(i)=SXeigenvec(i,best_vector_ovrlp_casscf)/c0 + print*,'',i,SXvector(i) end do END_PROVIDER @@ -113,12 +118,31 @@ BEGIN_PROVIDER [double precision, NewOrbs, (ao_num,mo_num) ] integer :: i,j,ialph if(state_following_casscf)then + print*,'Using the state following casscf ' call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & NatOrbsFCI, size(NatOrbsFCI,1), & Umat, size(Umat,1), 0.d0, & NewOrbs, size(NewOrbs,1)) else - NewOrbs = superci_natorb + double precision :: damp + print*,'Taking the lowest root for the CASSCF' + if(best_vector_ovrlp_casscf.ne.1)then + provide n_orb_swap + !call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & + ! NatOrbsFCI, size(NatOrbsFCI,1), & + ! Umat, size(Umat,1), 0.d0, & + ! NewOrbs, size(NewOrbs,1)) + NewOrbs = switch_mo_coef + mo_coef = switch_mo_coef + soft_touch mo_coef + call save_mos_no_occ + stop + else + call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & + NatOrbsFCI, size(NatOrbsFCI,1), & + Umat, size(Umat,1), 0.d0, & + NewOrbs, size(NewOrbs,1)) + endif endif END_PROVIDER diff --git a/src/casscf/superci_dm.irp.f b/src/casscf/superci_dm.irp.f new file mode 100644 index 00000000..0aef222b --- /dev/null +++ b/src/casscf/superci_dm.irp.f @@ -0,0 +1,207 @@ + BEGIN_PROVIDER [double precision, super_ci_dm, (mo_num,mo_num)] + implicit none + BEGIN_DOC +! density matrix of the super CI matrix, in the basis of NATURAL ORBITALS OF THE CASCI WF +! +! This is obtained from annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 +! +! WARNING ::: in the equation B3.d there is a TYPO with a forgotten MINUS SIGN (see variable mat_tmp_dm_super_ci ) + END_DOC + super_ci_dm = 0.d0 + integer :: i,j,iorb,jorb + integer :: a,aorb,b,borb + integer :: t,torb,v,vorb,u,uorb,x,xorb + double precision :: c0,ci + c0 = SXeigenvec(1,1) + ! equation B3.a of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + ! loop over the core/inact + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(iorb,iorb) = 2.d0 ! first term of B3.a + ! loop over the core/inact + do j = 1, n_core_inact_orb + jorb = list_core_inact(j) + ! loop over the virtual + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(jorb,iorb) += -2.d0 * lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,jorb) ! second term in B3.a + enddo + do t = 1, n_act_orb + torb = list_act(t) + ! thrid term of the B3.a + super_ci_dm(jorb,iorb) += - lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(jorb,torb) * (2.d0 - occ_act(t)) + enddo + enddo + enddo + + ! equation B3.b of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do t = 1, n_act_orb + torb = list_act(t) + super_ci_dm(iorb,torb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) + super_ci_dm(torb,iorb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(iorb,torb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + super_ci_dm(torb,iorb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + enddo + enddo + enddo + + ! equation B3.c of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(aorb,iorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) + super_ci_dm(iorb,aorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) + enddo + enddo + + ! equation B3.d of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do t = 1, n_act_orb + torb = list_act(t) + super_ci_dm(torb,torb) = occ_act(t) ! first term of equation B3.d + do x = 1, n_act_orb + xorb = list_act(x) + super_ci_dm(torb,torb) += - occ_act(x) * occ_act(t)* mat_tmp_dm_super_ci(x,x) ! second term involving the ONE-rdm + enddo + do u = 1, n_act_orb + uorb = list_act(u) + + ! second term of equation B3.d + do x = 1, n_act_orb + xorb = list_act(x) + do v = 1, n_act_orb + vorb = list_act(v) + super_ci_dm(torb,uorb) += 2.d0 * P0tuvx_no(v,x,t,u) * mat_tmp_dm_super_ci(v,x) ! second term involving the TWO-rdm + enddo + enddo + + ! third term of equation B3.d + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(torb,uorb) += lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(iorb,uorb) * (2.d0 - occ_act(t) - occ_act(u)) + enddo + + enddo + enddo + + ! equation B3.e of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do t = 1, n_act_orb + torb = list_act(t) + do a = 1, n_virt_orb + aorb = list_virt(a) + super_ci_dm(aorb,torb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + super_ci_dm(torb,aorb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(aorb,torb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) + super_ci_dm(torb,aorb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) + enddo + enddo + enddo + + ! equation B3.f of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 + do a = 1, n_virt_orb + aorb = list_virt(a) + do b = 1, n_virt_orb + borb= list_virt(b) + + ! First term of equation B3.f + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + super_ci_dm(borb,aorb) += 2.d0 * lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,borb) + enddo + + ! Second term of equation B3.f + do t = 1, n_act_orb + torb = list_act(t) + super_ci_dm(borb,aorb) += lowest_super_ci_coef_mo(torb,aorb) * lowest_super_ci_coef_mo(torb,borb) * occ_act(t) + enddo + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, superci_natorb, (ao_num,mo_num) +&BEGIN_PROVIDER [double precision, superci_nat_occ, (mo_num) + implicit none + call general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(super_ci_dm,mo_num,mo_num,mo_num,NatOrbsFCI,superci_nat_occ,superci_natorb) + +END_PROVIDER + + BEGIN_PROVIDER [double precision, mat_tmp_dm_super_ci, (n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC + ! computation of the term in [ ] in the equation B3.d of Roos et. al. Chemical Physics 48 (1980) 157-173 + ! + ! !!!!! WARNING !!!!!! there is a TYPO: a MINUS SIGN SHOULD APPEAR in that term + END_DOC + integer :: a,aorb,i,iorb + integer :: x,xorb,v,vorb + mat_tmp_dm_super_ci = 0.d0 + do v = 1, n_act_orb + vorb = list_act(v) + do x = 1, n_act_orb + xorb = list_act(x) + do a = 1, n_virt_orb + aorb = list_virt(a) + mat_tmp_dm_super_ci(x,v) += lowest_super_ci_coef_mo(aorb,vorb) * lowest_super_ci_coef_mo(aorb,xorb) + enddo + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + ! MARK THE MINUS SIGN HERE !!!!!!!!!!! BECAUSE OF TYPO IN THE ORIGINAL PAPER + mat_tmp_dm_super_ci(x,v) -= lowest_super_ci_coef_mo(iorb,vorb) * lowest_super_ci_coef_mo(iorb,xorb) + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, lowest_super_ci_coef_mo, (mo_num,mo_num)] + implicit none + integer :: i,j,iorb,jorb + integer :: a, aorb,t, torb + double precision :: sqrt2 + + sqrt2 = 1.d0/dsqrt(2.d0) + do i = 1, nMonoEx + iorb = excit(1,i) + jorb = excit(2,i) + lowest_super_ci_coef_mo(iorb,jorb) = SXeigenvec(i+1,1) + lowest_super_ci_coef_mo(jorb,iorb) = SXeigenvec(i+1,1) + enddo + + ! a_{it} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do t = 1, n_act_orb + torb = list_act(t) + lowest_super_ci_coef_mo(torb,iorb) *= (2.d0 - occ_act(t))**(-0.5d0) + lowest_super_ci_coef_mo(iorb,torb) *= (2.d0 - occ_act(t))**(-0.5d0) + enddo + enddo + + ! a_{ia} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 + do i = 1, n_core_inact_orb + iorb = list_core_inact(i) + do a = 1, n_virt_orb + aorb = list_virt(a) + lowest_super_ci_coef_mo(aorb,iorb) *= sqrt2 + lowest_super_ci_coef_mo(iorb,aorb) *= sqrt2 + enddo + enddo + + ! a_{ta} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 + do a = 1, n_virt_orb + aorb = list_virt(a) + do t = 1, n_act_orb + torb = list_act(t) + lowest_super_ci_coef_mo(torb,aorb) *= occ_act(t)**(-0.5d0) + lowest_super_ci_coef_mo(aorb,torb) *= occ_act(t)**(-0.5d0) + enddo + enddo + + END_PROVIDER + diff --git a/src/casscf/swap_orb.irp.f b/src/casscf/swap_orb.irp.f index c19bf811..6bbe733f 100644 --- a/src/casscf/swap_orb.irp.f +++ b/src/casscf/swap_orb.irp.f @@ -1,207 +1,122 @@ - BEGIN_PROVIDER [double precision, super_ci_dm, (mo_num,mo_num)] - implicit none - BEGIN_DOC -! density matrix of the super CI matrix, in the basis of NATURAL ORBITALS OF THE CASCI WF -! -! This is obtained from annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 -! -! WARNING ::: in the equation B3.d there is a TYPO with a forgotten MINUS SIGN (see variable mat_tmp_dm_super_ci ) - END_DOC - super_ci_dm = 0.d0 - integer :: i,j,iorb,jorb - integer :: a,aorb,b,borb - integer :: t,torb,v,vorb,u,uorb,x,xorb - double precision :: c0,ci - c0 = SXeigenvec(1,1) - ! equation B3.a of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - ! loop over the core/inact - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(iorb,iorb) = 2.d0 ! first term of B3.a - ! loop over the core/inact - do j = 1, n_core_inact_orb - jorb = list_core_inact(j) - ! loop over the virtual - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(jorb,iorb) += -2.d0 * lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,jorb) ! second term in B3.a - enddo - do t = 1, n_act_orb - torb = list_act(t) - ! thrid term of the B3.a - super_ci_dm(jorb,iorb) += - lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(jorb,torb) * (2.d0 - occ_act(t)) - enddo + BEGIN_PROVIDER [double precision, SXvector_lowest, (nMonoEx)] + implicit none + integer :: i + do i=2,nMonoEx+1 + SXvector_lowest(i-1)=SXeigenvec(i,1) enddo - enddo - - ! equation B3.b of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do t = 1, n_act_orb - torb = list_act(t) - super_ci_dm(iorb,torb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) - super_ci_dm(torb,iorb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(iorb,torb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - super_ci_dm(torb,iorb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - enddo - enddo - enddo - - ! equation B3.c of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(aorb,iorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) - super_ci_dm(iorb,aorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) - enddo - enddo - - ! equation B3.d of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do t = 1, n_act_orb - torb = list_act(t) - super_ci_dm(torb,torb) = occ_act(t) ! first term of equation B3.d - do x = 1, n_act_orb - xorb = list_act(x) - super_ci_dm(torb,torb) += - occ_act(x) * occ_act(t)* mat_tmp_dm_super_ci(x,x) ! second term involving the ONE-rdm - enddo - do u = 1, n_act_orb - uorb = list_act(u) - - ! second term of equation B3.d - do x = 1, n_act_orb - xorb = list_act(x) - do v = 1, n_act_orb - vorb = list_act(v) - super_ci_dm(torb,uorb) += 2.d0 * P0tuvx_no(v,x,t,u) * mat_tmp_dm_super_ci(v,x) ! second term involving the TWO-rdm - enddo - enddo - - ! third term of equation B3.d - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(torb,uorb) += lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(iorb,uorb) * (2.d0 - occ_act(t) - occ_act(u)) - enddo - - enddo - enddo - - ! equation B3.e of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do t = 1, n_act_orb - torb = list_act(t) - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(aorb,torb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - super_ci_dm(torb,aorb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(aorb,torb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) - super_ci_dm(torb,aorb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) - enddo - enddo - enddo - - ! equation B3.f of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do a = 1, n_virt_orb - aorb = list_virt(a) - do b = 1, n_virt_orb - borb= list_virt(b) - - ! First term of equation B3.f - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(borb,aorb) += 2.d0 * lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,borb) - enddo - - ! Second term of equation B3.f - do t = 1, n_act_orb - torb = list_act(t) - super_ci_dm(borb,aorb) += lowest_super_ci_coef_mo(torb,aorb) * lowest_super_ci_coef_mo(torb,borb) * occ_act(t) - enddo - enddo - enddo - END_PROVIDER - BEGIN_PROVIDER [double precision, superci_natorb, (ao_num,mo_num) -&BEGIN_PROVIDER [double precision, superci_nat_occ, (mo_num) + BEGIN_PROVIDER [double precision, thresh_overlap_switch] implicit none - call general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(super_ci_dm,mo_num,mo_num,mo_num,superci_nat_occ,superci_natorb) - -END_PROVIDER - - BEGIN_PROVIDER [double precision, mat_tmp_dm_super_ci, (n_act_orb,n_act_orb)] - implicit none - BEGIN_DOC - ! computation of the term in [ ] in the equation B3.d of Roos et. al. Chemical Physics 48 (1980) 157-173 - ! - ! !!!!! WARNING !!!!!! there is a TYPO: a MINUS SIGN SHOULD APPEAR in that term - END_DOC - integer :: a,aorb,i,iorb - integer :: x,xorb,v,vorb - mat_tmp_dm_super_ci = 0.d0 - do v = 1, n_act_orb - vorb = list_act(v) - do x = 1, n_act_orb - xorb = list_act(x) - do a = 1, n_virt_orb - aorb = list_virt(a) - mat_tmp_dm_super_ci(x,v) += lowest_super_ci_coef_mo(aorb,vorb) * lowest_super_ci_coef_mo(aorb,xorb) - enddo - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - ! MARK THE MINUS SIGN HERE !!!!!!!!!!! BECAUSE OF TYPO IN THE ORIGINAL PAPER - mat_tmp_dm_super_ci(x,v) -= lowest_super_ci_coef_mo(iorb,vorb) * lowest_super_ci_coef_mo(iorb,xorb) - enddo - enddo - enddo + thresh_overlap_switch = 0.5d0 END_PROVIDER - - BEGIN_PROVIDER [double precision, lowest_super_ci_coef_mo, (mo_num,mo_num)] - implicit none - integer :: i,j,iorb,jorb - integer :: a, aorb,t, torb - double precision :: sqrt2 - sqrt2 = 1.d0/dsqrt(2.d0) + BEGIN_PROVIDER [integer, max_overlap, (nMonoEx)] +&BEGIN_PROVIDER [integer, n_max_overlap] + implicit none + double precision, allocatable :: vec_tmp(:) + integer, allocatable :: iorder(:) + allocate(vec_tmp(nMonoEx),iorder(nMonoEx)) + integer :: i do i = 1, nMonoEx - iorb = excit(1,i) - jorb = excit(2,i) - lowest_super_ci_coef_mo(iorb,jorb) = SXeigenvec(i+1,1) - lowest_super_ci_coef_mo(jorb,iorb) = SXeigenvec(i+1,1) + iorder(i) = i + vec_tmp(i) = -dabs(SXvector_lowest(i)) + print*,'vec_tmp(i) = ',i,vec_tmp(i) enddo - - ! a_{it} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do t = 1, n_act_orb - torb = list_act(t) - lowest_super_ci_coef_mo(torb,iorb) *= (2.d0 - occ_act(t))**(-0.5d0) - lowest_super_ci_coef_mo(iorb,torb) *= (2.d0 - occ_act(t))**(-0.5d0) - enddo + call dsort(vec_tmp,iorder,nMonoEx) + n_max_overlap = 0 + do i = 1, nMonoEx + if(dabs(vec_tmp(i)).gt.thresh_overlap_switch)then + print*,vec_tmp(i),iorder(i) + n_max_overlap += 1 + max_overlap(n_max_overlap) = iorder(i) + endif enddo - - ! a_{ia} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do a = 1, n_virt_orb - aorb = list_virt(a) - lowest_super_ci_coef_mo(aorb,iorb) *= sqrt2 - lowest_super_ci_coef_mo(iorb,aorb) *= sqrt2 - enddo - enddo - - ! a_{ta} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 - do a = 1, n_virt_orb - aorb = list_virt(a) - do t = 1, n_act_orb - torb = list_act(t) - lowest_super_ci_coef_mo(torb,aorb) *= occ_act(t)**(-0.5d0) - lowest_super_ci_coef_mo(aorb,torb) *= occ_act(t)**(-0.5d0) - enddo - enddo - END_PROVIDER + BEGIN_PROVIDER [integer, orb_swap, (2,n_max_overlap)] +&BEGIN_PROVIDER [integer, n_orb_swap ] + implicit none + use bitmasks ! you need to include the bitmasks_module.f90 features + integer :: i,imono,iorb,jorb,j + n_orb_swap = 0 + do i = 1, n_max_overlap + imono = max_overlap(i) + iorb = excit(1,imono) + jorb = excit(2,imono) + if (excit_class(imono) == "c-a")then ! core --> active rotation + n_orb_swap += 1 + orb_swap(1,n_orb_swap) = iorb ! core + orb_swap(2,n_orb_swap) = jorb ! active + else if (excit_class(imono) == "a-v")then ! active --> virtual rotation + n_orb_swap += 1 + orb_swap(1,n_orb_swap) = jorb ! virtual + orb_swap(2,n_orb_swap) = iorb ! active + endif + enddo + print*,'n_orb_swap = ',n_orb_swap + do i = 1, n_orb_swap + print*,orb_swap(1,i),'-->',orb_swap(2,i) + enddo + + orb_swap_tmp = orb_swap + integer :: orb_swap_tmp(2,n_max_overlap) + integer(bit_kind), allocatable :: det_i(:),det_j(:) + allocate(det_i(N_int),det_j(N_int)) + logical, allocatable :: good_orb_rot(:) + allocate(good_orb_rot(n_orb_swap)) + good_orb_rot = .True. + integer :: icount,k + do i = 1, n_orb_swap + if(.not.good_orb_rot(i))cycle + det_i = 0_bit_kind + call set_bit_to_integer(orb_swap(1,i),det_i,N_int) + call set_bit_to_integer(orb_swap(2,i),det_i,N_int) + do j = i+1, n_orb_swap + det_j = 0_bit_kind + call set_bit_to_integer(orb_swap(1,j),det_j,N_int) + call set_bit_to_integer(orb_swap(2,j),det_j,N_int) + icount = 0 + do k = 1, N_int + icount += popcnt(ior(det_i(k),det_j(k))) + enddo + if (icount.ne.4)then + good_orb_rot(i) = .False. + good_orb_rot(j) = .False. + exit + endif + enddo + enddo + icount = n_orb_swap + n_orb_swap = 0 + do i = 1, icount + if(good_orb_rot(i))then + n_orb_swap += 1 + orb_swap(1,n_orb_swap) = orb_swap_tmp(1,i) + orb_swap(2,n_orb_swap) = orb_swap_tmp(2,i) + endif + enddo + + print*,'Cleaning !!' + print*,'n_orb_swap = ',n_orb_swap + do i = 1, n_orb_swap + print*,orb_swap(1,i),'-->',orb_swap(2,i) + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, switch_mo_coef, (ao_num,mo_num)] + implicit none + integer :: i,j,iorb,jorb + switch_mo_coef = NatOrbsFCI + do i = 1, n_orb_swap + iorb = orb_swap(1,i) + jorb = orb_swap(2,i) + do j = 1, ao_num + switch_mo_coef(j,jorb) = NatOrbsFCI(j,iorb) + enddo + do j = 1, ao_num + switch_mo_coef(j,iorb) = NatOrbsFCI(j,jorb) + enddo + enddo + + END_PROVIDER diff --git a/src/mo_basis/utils.irp.f b/src/mo_basis/utils.irp.f index 2cefcbd4..66ef5c1e 100644 --- a/src/mo_basis/utils.irp.f +++ b/src/mo_basis/utils.irp.f @@ -4,7 +4,7 @@ subroutine save_mos integer :: i,j call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) - + print*,'Saving MOs' call ezfio_set_mo_basis_mo_num(mo_num) call ezfio_set_mo_basis_mo_label(mo_label) call ezfio_set_mo_basis_ao_md5(ao_md5) @@ -18,6 +18,31 @@ subroutine save_mos call ezfio_set_mo_basis_mo_coef(buffer) call ezfio_set_mo_basis_mo_occ(mo_occ) deallocate (buffer) + print*,'End Saving MOs' + +end + + +subroutine save_mos_no_occ + implicit none + double precision, allocatable :: buffer(:,:) + integer :: i,j + + call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) + print*,'Saving MOs' + !call ezfio_set_mo_basis_mo_num(mo_num) + !call ezfio_set_mo_basis_mo_label(mo_label) + !call ezfio_set_mo_basis_ao_md5(ao_md5) + allocate ( buffer(ao_num,mo_num) ) + buffer = 0.d0 + do j = 1, mo_num + do i = 1, ao_num + buffer(i,j) = mo_coef(i,j) + enddo + enddo + call ezfio_set_mo_basis_mo_coef(buffer) + deallocate (buffer) + print*,'End Saving MOs' end @@ -217,10 +242,10 @@ subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label) end -subroutine general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,mo_coef_new) +subroutine general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,mo_coef_before,eig,mo_coef_new) implicit none integer,intent(in) :: lda,m,n - double precision, intent(in) :: matrix(lda,n) + double precision, intent(in) :: matrix(lda,n),mo_coef_before(ao_num,m) double precision, intent(out) :: eig(m),mo_coef_new(ao_num,m) integer :: i,j @@ -241,7 +266,7 @@ subroutine general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,ei A(i,j) = matrix(i,j) enddo enddo - mo_coef_tmp = mo_coef + mo_coef_tmp = mo_coef_before call svd(A,lda,U,lda,D,Vt,lda,m,n) @@ -261,7 +286,7 @@ subroutine general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,ei write (6,'(A)') '======== ================ ================' write (6,'(A)') '' - call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_tmp,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef_new,size(mo_coef,1)) + call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_tmp,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef_new,size(mo_coef_new,1)) do i=1,m eig(i) = D(i)