From 90e9d4bf7e282946742106c3cebaadc4e0268daa Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 30 Jun 2021 18:31:03 +0200 Subject: [PATCH] casscf works with full integral transformation --- devel/casscf/EZFIO.cfg | 6 ++ devel/casscf/TODO | 1 + devel/casscf/casscf.irp.f | 6 +- devel/casscf/densities.irp.f | 24 +++-- devel/casscf/densities_peter.irp.f | 150 +++++++++++++++++++++++++++++ 5 files changed, 174 insertions(+), 13 deletions(-) create mode 100644 devel/casscf/TODO create mode 100644 devel/casscf/densities_peter.irp.f diff --git a/devel/casscf/EZFIO.cfg b/devel/casscf/EZFIO.cfg index 4e4d3d3..7ff49a5 100644 --- a/devel/casscf/EZFIO.cfg +++ b/devel/casscf/EZFIO.cfg @@ -29,3 +29,9 @@ doc: Energy shift on the virtual MOs to improve SCF convergence interface: ezfio,provider,ocaml default: 0.005 + +[fast_2rdm] +type: logical +doc: If true, the two-rdm are computed with a fast algo +interface: ezfio,provider,ocaml +default: True diff --git a/devel/casscf/TODO b/devel/casscf/TODO new file mode 100644 index 0000000..e698199 --- /dev/null +++ b/devel/casscf/TODO @@ -0,0 +1 @@ +it recomputes the gradients and hessian also with only one determinant, useless and confusing diff --git a/devel/casscf/casscf.irp.f b/devel/casscf/casscf.irp.f index 950cfd5..76a6fcb 100644 --- a/devel/casscf/casscf.irp.f +++ b/devel/casscf/casscf.irp.f @@ -4,10 +4,10 @@ program casscf ! TODO : Put the documentation of the program here END_DOC call reorder_orbitals_for_casscf - no_vvvv_integrals = .True. - touch no_vvvv_integrals +! no_vvvv_integrals = .True. +! touch no_vvvv_integrals pt2_max = 0.02 - SOFT_TOUCH no_vvvv_integrals pt2_max + SOFT_TOUCH pt2_max call run_stochastic_cipsi call run end diff --git a/devel/casscf/densities.irp.f b/devel/casscf/densities.irp.f index d181d73..eb8e682 100644 --- a/devel/casscf/densities.irp.f +++ b/devel/casscf/densities.irp.f @@ -47,17 +47,21 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] endif P0tuvx= 0.d0 - do istate=1,N_states - do x = 1, n_act_orb - do v = 1, n_act_orb - do u = 1, n_act_orb - do t = 1, n_act_orb - ! 1 1 2 2 1 2 1 2 - P0tuvx(t,u,v,x) = state_av_act_2_rdm_spin_trace_mo(t,v,u,x) - enddo - enddo + if(fast_2rdm)then + do istate=1,N_states + do x = 1, n_act_orb + do v = 1, n_act_orb + do u = 1, n_act_orb + do t = 1, n_act_orb + ! 1 1 2 2 1 2 1 2 + P0tuvx(t,u,v,x) = state_av_act_2_rdm_spin_trace_mo(t,v,u,x) + enddo + enddo + enddo enddo enddo - enddo + else + P0tuvx = P0tuvx_peter + endif END_PROVIDER diff --git a/devel/casscf/densities_peter.irp.f b/devel/casscf/densities_peter.irp.f new file mode 100644 index 0000000..ee7414d --- /dev/null +++ b/devel/casscf/densities_peter.irp.f @@ -0,0 +1,150 @@ +use bitmasks + +BEGIN_PROVIDER [real*8, P0tuvx_peter, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] + BEGIN_DOC + ! the second-order density matrix in the basis of the starting MOs + ! matrices are state averaged + ! + ! we use the spin-free generators of mono-excitations + ! E_pq destroys q and creates p + ! D_pq = <0|E_pq|0> = D_qp + ! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0> + ! + END_DOC + implicit none + integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart + integer :: ierr + real*8 :: phase1,phase11,phase12,phase2,phase21,phase22 + integer :: nu1,nu2,nu11,nu12,nu21,nu22 + integer :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22 + real*8 :: cI_mu(N_states),term + integer(bit_kind), dimension(N_int,2) :: det_mu, det_mu_ex + integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12 + integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22 + + if (bavard) then + write(6,*) ' providing density matrix P0' + endif + + P0tuvx_peter = 0.d0 + + ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu + do mu=1,n_det + call det_extract(det_mu,mu,N_int) + do istate=1,n_states + cI_mu(istate)=psi_coef(mu,istate) + end do + do t=1,n_act_orb + ipart=list_act(t) + do u=1,n_act_orb + ihole=list_act(u) + ! apply E_tu + call det_copy(det_mu,det_mu_ex1,N_int) + call det_copy(det_mu,det_mu_ex2,N_int) + call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & + ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) + ! det_mu_ex1 is in the list + if (nu1.ne.-1) then + do istate=1,n_states + term=cI_mu(istate)*psi_coef(nu1,istate)*phase1 + ! and we fill P0_tvvu + do v=1,n_act_orb + P0tuvx_peter(t,v,v,u)-=term + end do + end do + end if + ! det_mu_ex2 is in the list + if (nu2.ne.-1) then + do istate=1,n_states + term=cI_mu(istate)*psi_coef(nu2,istate)*phase2 + do v=1,n_act_orb + P0tuvx_peter(t,v,v,u)-=term + end do + end do + end if + end do + end do + end do + ! now we do the double excitation E_tu E_vx |0> + do mu=1,n_det + call det_extract(det_mu,mu,N_int) + do istate=1,n_states + cI_mu(istate)=psi_coef(mu,istate) + end do + do v=1,n_act_orb + ipart=list_act(v) + do x=1,n_act_orb + ihole=list_act(x) + ! apply E_vx + call det_copy(det_mu,det_mu_ex1,N_int) + call det_copy(det_mu,det_mu_ex2,N_int) + call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & + ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) + ! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0> + if (ierr1.eq.1) then + do t=1,n_act_orb + jpart=list_act(t) + do u=1,n_act_orb + jhole=list_act(u) + call det_copy(det_mu_ex1,det_mu_ex11,N_int) + call det_copy(det_mu_ex1,det_mu_ex12,N_int) + call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11& + ,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12) + if (nu11.ne.-1) then + do istate=1,n_states + P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate)& + *phase11*phase1 + end do + end if + if (nu12.ne.-1) then + do istate=1,n_states + P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate)& + *phase12*phase1 + end do + end if + end do + end do + end if + + ! we apply E_tu to the second resultant determinant + if (ierr2.eq.1) then + do t=1,n_act_orb + jpart=list_act(t) + do u=1,n_act_orb + jhole=list_act(u) + call det_copy(det_mu_ex2,det_mu_ex21,N_int) + call det_copy(det_mu_ex2,det_mu_ex22,N_int) + call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21& + ,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22) + if (nu21.ne.-1) then + do istate=1,n_states + P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate)& + *phase21*phase2 + end do + end if + if (nu22.ne.-1) then + do istate=1,n_states + P0tuvx_peter(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate)& + *phase22*phase2 + end do + end if + end do + end do + end if + + end do + end do + end do + + ! we average by just dividing by the number of states + do x=1,n_act_orb + do v=1,n_act_orb + do u=1,n_act_orb + do t=1,n_act_orb + P0tuvx_peter(t,u,v,x)*=0.5D0/dble(N_states) + end do + end do + end do + end do + +END_PROVIDER