From 08b3f247f040c793ef644a7bb22585fb3b4f5b4e Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 2 Jul 2021 17:46:24 +0200 Subject: [PATCH] added the possibility to diagonalize the SXmatrix with Davidson --- devel/casscf/EZFIO.cfg | 6 +++++ devel/casscf/NEED | 1 + devel/casscf/casscf.irp.f | 15 ++++++++---- devel/casscf/dav_sx_mat.irp.f | 44 +++++++++++++++++++++++++++++++++++ devel/casscf/gradient.irp.f | 8 +++++++ devel/casscf/hessian.irp.f | 42 +++++++++++++++++---------------- devel/casscf/neworbs.irp.f | 43 +++++++++++++++++++++++++++------- 7 files changed, 126 insertions(+), 33 deletions(-) create mode 100644 devel/casscf/dav_sx_mat.irp.f diff --git a/devel/casscf/EZFIO.cfg b/devel/casscf/EZFIO.cfg index 332c2d4..ff1f7b3 100644 --- a/devel/casscf/EZFIO.cfg +++ b/devel/casscf/EZFIO.cfg @@ -29,6 +29,12 @@ doc: If |true|, only the DIAGONAL part of the hessian is retained for the CASSCF interface: ezfio,provider,ocaml default: False +[hess_cv_cv] +type: logical +doc: If |true|, the core-virtual - core-virtual part of the hessian is computed +interface: ezfio,provider,ocaml +default: True + [level_shift_casscf] type: Positive_float diff --git a/devel/casscf/NEED b/devel/casscf/NEED index d9da718..dd91c7b 100644 --- a/devel/casscf/NEED +++ b/devel/casscf/NEED @@ -2,3 +2,4 @@ cipsi selectors_full generators_cas two_body_rdm +dav_general_mat diff --git a/devel/casscf/casscf.irp.f b/devel/casscf/casscf.irp.f index 355cd1d..27ffe74 100644 --- a/devel/casscf/casscf.irp.f +++ b/devel/casscf/casscf.irp.f @@ -6,8 +6,10 @@ program casscf call reorder_orbitals_for_casscf no_vvvv_integrals = .True. touch no_vvvv_integrals - pt2_max = 0.02 + pt2_max = 0.005 SOFT_TOUCH pt2_max + n_det_max_full = 500 + touch n_det_max_full call run_stochastic_cipsi call run end @@ -31,10 +33,13 @@ subroutine run energy = eone+etwo+ecore call write_time(6) - call write_int(6,iteration,'CAS-SCF iteration') - call write_double(6,energy,'CAS-SCF energy') - call write_double(6,norm_grad_vec2,'Norm of gradients') - call write_double(6,energy_improvement, 'Predicted energy improvement') + call write_int(6,iteration,'CAS-SCF iteration = ') + call write_double(6,energy,'CAS-SCF energy = ') + call write_double(6,norm_grad_vec2,'Norm of gradients = ') + call write_double(6,norm_grad_vec2_tab(1), ' Core-active gradients = ') + call write_double(6,norm_grad_vec2_tab(2), ' Core-virtual gradients = ') + call write_double(6,norm_grad_vec2_tab(3), ' Active-virtual gradients = ') + call write_double(6,energy_improvement, 'Predicted energy improvement = ') converged = dabs(energy_improvement) < thresh_scf pt2_max = dabs(energy_improvement / pt2_relative_error) diff --git a/devel/casscf/dav_sx_mat.irp.f b/devel/casscf/dav_sx_mat.irp.f new file mode 100644 index 0000000..1e83cdf --- /dev/null +++ b/devel/casscf/dav_sx_mat.irp.f @@ -0,0 +1,44 @@ + + +subroutine davidson_diag_sx_mat(N_st, u_in, energies) + implicit none + integer, intent(in) :: N_st + double precision, intent(out) :: u_in(nMonoEx+1,n_states_diag), energies(N_st) + integer :: i,j,N_st_tmp, dim_in, sze, N_st_diag_in + integer, allocatable :: list_guess(:) + double precision, allocatable :: H_jj(:) + logical :: converged + N_st_diag_in = n_states_diag + provide SXmatrix + sze = nMonoEx+1 + dim_in = sze + allocate(H_jj(sze), list_guess(sze)) + H_jj(1) = 0.d0 + N_st_tmp = 1 + list_guess(1) = 1 + do j = 2, nMonoEx+1 + H_jj(j) = SXmatrix(j,j) + if(H_jj(j).lt.0.d0)then + list_guess(N_st_tmp) = j + N_st_tmp += 1 + endif + enddo + if(N_st_tmp .ne. N_st)then + print*,'Pb in davidson_diag_sx_mat' + print*,'N_st_tmp .ne. N_st' + print*,N_st_tmp, N_st + stop + endif + print*,'Number of possibly interesting states = ',N_st + print*,'Corresponding diagonal elements of the SX matrix ' + u_in = 0.d0 + do i = 1, N_st + j = list_guess(i) + print*,'i,j',i,j + print*,'SX(i,i) = ',H_jj(j) + u_in(j,i) = 1.d0 + enddo + call davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,SXmatrix) + print*,'energies = ',energies + +end diff --git a/devel/casscf/gradient.irp.f b/devel/casscf/gradient.irp.f index 309371b..4c0caaf 100644 --- a/devel/casscf/gradient.irp.f +++ b/devel/casscf/gradient.irp.f @@ -95,6 +95,7 @@ END_PROVIDER BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] &BEGIN_PROVIDER [real*8, norm_grad_vec2] +&BEGIN_PROVIDER [real*8, norm_grad_vec2_tab, (3)] BEGIN_DOC ! calculate the orbital gradient from density ! matrices and integrals; Siegbahn et al, Phys Scr 1980 @@ -105,10 +106,12 @@ END_PROVIDER real*8 :: gradvec_it,gradvec_ia,gradvec_ta indx=0 + norm_grad_vec2_tab = 0.d0 do i=1,n_core_inact_orb do t=1,n_act_orb indx+=1 gradvec2(indx)=gradvec_it(i,t) + norm_grad_vec2_tab(1) += gradvec2(indx)*gradvec2(indx) end do end do @@ -116,6 +119,7 @@ END_PROVIDER do a=1,n_virt_orb indx+=1 gradvec2(indx)=gradvec_ia(i,a) + norm_grad_vec2_tab(2) += gradvec2(indx)*gradvec2(indx) end do end do @@ -123,6 +127,7 @@ END_PROVIDER do a=1,n_virt_orb indx+=1 gradvec2(indx)=gradvec_ta(t,a) + norm_grad_vec2_tab(3) += gradvec2(indx)*gradvec2(indx) end do end do @@ -130,6 +135,9 @@ END_PROVIDER do indx=1,nMonoEx norm_grad_vec2+=gradvec2(indx)*gradvec2(indx) end do + do i = 1, 3 + norm_grad_vec2_tab(i) = dsqrt(norm_grad_vec2_tab(i)) + enddo norm_grad_vec2=sqrt(norm_grad_vec2) if(bavard)then write(6,*) diff --git a/devel/casscf/hessian.irp.f b/devel/casscf/hessian.irp.f index 61864ef..458c6aa 100644 --- a/devel/casscf/hessian.irp.f +++ b/devel/casscf/hessian.irp.f @@ -453,34 +453,36 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] !$OMP END DO NOWAIT !$OMP END PARALLEL + + if(hess_cv_cv)then !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(hessmat,n_c_v_prov,list_idx_c_v,n_core_inact_orb,n_virt_orb,mat_idx_c_v) & !$OMP PRIVATE(indx_tmp,indx,i,a,j,b,bstart,jndx) - - !$OMP DO -!!!! < Core-VIRTUAL | H |Core-VIRTUAL > - ! Core-VIRTUAL excitations - do indx_tmp = 1, n_c_v_prov - indx = list_idx_c_v(1,indx_tmp) - i = list_idx_c_v(2,indx_tmp) - a = list_idx_c_v(3,indx_tmp) + !$OMP DO +!!!!! < Core-VIRTUAL | H |Core-VIRTUAL > ! Core-VIRTUAL excitations - do j = 1, n_core_inact_orb - if (i.eq.j) then - bstart=a - else - bstart=1 - end if - do b=bstart,n_virt_orb - jndx = mat_idx_c_v(j,b) - hessmat(jndx,indx) = hessmat_iajb(i,a,j,b) - hessmat(indx,jndx) = hessmat(jndx,indx) + do indx_tmp = 1, n_c_v_prov + indx = list_idx_c_v(1,indx_tmp) + i = list_idx_c_v(2,indx_tmp) + a = list_idx_c_v(3,indx_tmp) + ! Core-VIRTUAL excitations + do j = 1, n_core_inact_orb + if (i.eq.j) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + jndx = mat_idx_c_v(j,b) + hessmat(jndx,indx) = hessmat_iajb(i,a,j,b) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo enddo enddo - enddo - + !$OMP END DO NOWAIT !$OMP END PARALLEL + endif !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(hessmat,n_c_v_prov,n_a_v_prov,list_idx_c_v,list_idx_a_v) & diff --git a/devel/casscf/neworbs.irp.f b/devel/casscf/neworbs.irp.f index 3970395..d19e556 100644 --- a/devel/casscf/neworbs.irp.f +++ b/devel/casscf/neworbs.irp.f @@ -1,4 +1,5 @@ -BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] + BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] +&BEGIN_PROVIDER [integer, n_guess_sx_mat ] implicit none BEGIN_DOC ! Single-excitation matrix @@ -32,13 +33,18 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] do i = 1, nMonoEx SXmatrix(i+1,i+1) += level_shift_casscf enddo + n_guess_sx_mat = 1 + do i = 1, nMonoEx + if(SXmatrix(i+1,i+1).lt.0.d0 )then + n_guess_sx_mat += 1 + endif + enddo if (bavard) then do i=2,nMonoEx write(6,*) ' diagonal of the Hessian : ',i,hessmat(i,i) end do end if - END_PROVIDER BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)] @@ -47,19 +53,40 @@ END_PROVIDER BEGIN_DOC ! Eigenvectors/eigenvalues of the single-excitation matrix END_DOC - call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1) + if(nMonoEx+1.gt.n_det_max_full)then +! if(bavard)then + print*,'Using the Davidson algorithm to diagonalize the SXmatrix' +! endif + double precision, allocatable :: u_in(:,:),energies(:) + allocate(u_in(nMonoEx+1,n_states_diag),energies(n_guess_sx_mat)) + call davidson_diag_sx_mat(n_guess_sx_mat, u_in, energies) + integer :: i,j + SXeigenvec = 0.d0 + SXeigenval = 0.d0 + do i = 1, n_guess_sx_mat + SXeigenval(i) = energies(i) + do j = 1, nMonoEx+1 + SXeigenvec(j,i) = u_in(j,i) + enddo + enddo + else +! if(bavard)then + print*,'Diagonalize the SXmatrix with Jacobi' +! endif + call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1) + endif if (bavard) then - write(6,*) ' SXdiag : lowest 5 eigenvalues ' + write(6,*) ' SXdiag : lowest eigenvalues ' write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1) - if(nmonoex.gt.0)then + if(n_guess_sx_mat.gt.0)then write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2) write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3) write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4) write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5) endif + endif write(6,*) write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1) - endif END_PROVIDER BEGIN_PROVIDER [real*8, energy_improvement] @@ -82,8 +109,8 @@ END_PROVIDER best_vector_ovrlp_casscf = -1000 do i=1,nMonoEx+1 if (SXeigenval(i).lt.0.D0) then - if (abs(SXeigenvec(1,i)).gt.best_overlap_casscf) then - best_overlap_casscf=abs(SXeigenvec(1,i)) + if (dabs(SXeigenvec(1,i)).gt.best_overlap_casscf) then + best_overlap_casscf=dabs(SXeigenvec(1,i)) best_vector_ovrlp_casscf = i end if end if