From 71f6163c40d70b4f35bd65f221f4da7b370149df Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 18 Jun 2023 20:28:48 +0200 Subject: [PATCH 01/21] added some comments for normal ordering old --- src/tc_bi_ortho/normal_ordered_old.irp.f | 10 +++++- src/tc_bi_ortho/test_normal_order.irp.f | 43 ++++++++++++++++++++---- 2 files changed, 45 insertions(+), 8 deletions(-) diff --git a/src/tc_bi_ortho/normal_ordered_old.irp.f b/src/tc_bi_ortho/normal_ordered_old.irp.f index 417580dd..6ee21a14 100644 --- a/src/tc_bi_ortho/normal_ordered_old.irp.f +++ b/src/tc_bi_ortho/normal_ordered_old.irp.f @@ -120,6 +120,13 @@ END_PROVIDER subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) use bitmasks ! you need to include the bitmasks_module.f90 features + BEGIN_DOC +! give the contribution for a double excitation of opposite spin BUT averaged over spin +! +! it is the average of and +! +! because the orbitals h1,h2,p1,p2 are spatial orbitals and therefore can be of different spins + END_DOC implicit none integer, intent(in) :: Nint, h1, h2, p1, p2 @@ -158,7 +165,8 @@ subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) int_exc_12 = -1.d0 * integral - hthree += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12) + hthree += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12) ! spin average +! hthree += 1.d0 * int_direct - 1.0d0 * (int_exc_13 + int_exc_12) enddo return diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f index cb0c355c..ac84dbc6 100644 --- a/src/tc_bi_ortho/test_normal_order.irp.f +++ b/src/tc_bi_ortho/test_normal_order.irp.f @@ -20,7 +20,7 @@ subroutine test integer :: h1,h2,p1,p2,s1,s2,i_ok,degree,Ne(2) integer :: exc(0:2,2,2) integer(bit_kind), allocatable :: det_i(:,:) - double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal + double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal,hthree_tmp integer, allocatable :: occ(:,:) allocate( occ(N_int*bit_kind_size,2) ) call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int) @@ -32,15 +32,44 @@ subroutine test do p1 = elec_alpha_num+1, mo_num do h2 = 1, elec_beta_num do p2 = elec_beta_num+1, mo_num + hthree = 0.d0 + det_i = ref_bitmask + s1 = 1 + s2 = 2 call do_single_excitation(det_i,h1,p1,s1,i_ok) + if(i_ok.ne.1)cycle call do_single_excitation(det_i,h2,p2,s2,i_ok) - call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + if(i_ok.ne.1)cycle + call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,htilde_ij) call get_excitation_degree(ref_bitmask,det_i,degree,N_int) call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) - hthree *= phase -! !normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) - call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal) + hthree_tmp *= phase + hthree += 0.5d0 * hthree_tmp + det_i = ref_bitmask + s1 = 2 + s2 = 1 + call do_single_excitation(det_i,h1,p1,s1,i_ok) + if(i_ok.ne.1)cycle + call do_single_excitation(det_i,h2,p2,s2,i_ok) + if(i_ok.ne.1)cycle + call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,htilde_ij) + call get_excitation_degree(ref_bitmask,det_i,degree,N_int) + call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) + hthree_tmp *= phase + hthree += 0.5d0 * hthree_tmp + + +! normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) + call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, normal) + if(dabs(hthree).lt.1.d-10)cycle + if(dabs(hthree-normal).gt.1.d-10)then +! print*,pp2,pp1,hh2,hh1 + print*,p2,p1,h2,h1 + print*,hthree,normal,dabs(hthree-normal) + stop + endif +! call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal) ! normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1) accu += dabs(hthree-normal) enddo @@ -73,8 +102,8 @@ do h1 = 1, elec_alpha_num integer :: hh1, pp1, hh2, pp2, ss1, ss2 call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2) hthree *= phase -! normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) - normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1) + normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) +! normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1) if(dabs(hthree).lt.1.d-10)cycle if(dabs(hthree-normal).gt.1.d-10)then print*,pp2,pp1,hh2,hh1 From b2e44beb3e11cacc9594a58da5c7ed4295506092 Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 18 Jun 2023 21:42:40 +0200 Subject: [PATCH 02/21] added casscf_cipsi --- src/casscf_cipsi/50.casscf.bats | 49 +++ src/casscf_cipsi/EZFIO.cfg | 75 ++++ src/casscf_cipsi/NEED | 5 + src/casscf_cipsi/README.rst | 5 + src/casscf_cipsi/bavard.irp.f | 6 + src/casscf_cipsi/bielec.irp.f | 155 +++++++ src/casscf_cipsi/bielec_natorb.irp.f | 369 +++++++++++++++++ src/casscf_cipsi/casscf.irp.f | 110 +++++ src/casscf_cipsi/class.irp.f | 12 + src/casscf_cipsi/dav_sx_mat.irp.f | 45 +++ src/casscf_cipsi/densities.irp.f | 67 +++ src/casscf_cipsi/densities_peter.irp.f | 150 +++++++ src/casscf_cipsi/det_manip.irp.f | 125 ++++++ src/casscf_cipsi/driver_optorb.irp.f | 3 + src/casscf_cipsi/get_energy.irp.f | 51 +++ src/casscf_cipsi/grad_old.irp.f | 74 ++++ src/casscf_cipsi/gradient.irp.f | 215 ++++++++++ src/casscf_cipsi/hessian.irp.f | 539 +++++++++++++++++++++++++ src/casscf_cipsi/hessian_old.irp.f | 310 ++++++++++++++ src/casscf_cipsi/mcscf_fock.irp.f | 80 ++++ src/casscf_cipsi/natorb.irp.f | 231 +++++++++++ src/casscf_cipsi/neworbs.irp.f | 253 ++++++++++++ src/casscf_cipsi/reorder_orb.irp.f | 70 ++++ src/casscf_cipsi/save_energy.irp.f | 9 + src/casscf_cipsi/superci_dm.irp.f | 207 ++++++++++ src/casscf_cipsi/swap_orb.irp.f | 132 ++++++ src/casscf_cipsi/tot_en.irp.f | 101 +++++ 27 files changed, 3448 insertions(+) create mode 100644 src/casscf_cipsi/50.casscf.bats create mode 100644 src/casscf_cipsi/EZFIO.cfg create mode 100644 src/casscf_cipsi/NEED create mode 100644 src/casscf_cipsi/README.rst create mode 100644 src/casscf_cipsi/bavard.irp.f create mode 100644 src/casscf_cipsi/bielec.irp.f create mode 100644 src/casscf_cipsi/bielec_natorb.irp.f create mode 100644 src/casscf_cipsi/casscf.irp.f create mode 100644 src/casscf_cipsi/class.irp.f create mode 100644 src/casscf_cipsi/dav_sx_mat.irp.f create mode 100644 src/casscf_cipsi/densities.irp.f create mode 100644 src/casscf_cipsi/densities_peter.irp.f create mode 100644 src/casscf_cipsi/det_manip.irp.f create mode 100644 src/casscf_cipsi/driver_optorb.irp.f create mode 100644 src/casscf_cipsi/get_energy.irp.f create mode 100644 src/casscf_cipsi/grad_old.irp.f create mode 100644 src/casscf_cipsi/gradient.irp.f create mode 100644 src/casscf_cipsi/hessian.irp.f create mode 100644 src/casscf_cipsi/hessian_old.irp.f create mode 100644 src/casscf_cipsi/mcscf_fock.irp.f create mode 100644 src/casscf_cipsi/natorb.irp.f create mode 100644 src/casscf_cipsi/neworbs.irp.f create mode 100644 src/casscf_cipsi/reorder_orb.irp.f create mode 100644 src/casscf_cipsi/save_energy.irp.f create mode 100644 src/casscf_cipsi/superci_dm.irp.f create mode 100644 src/casscf_cipsi/swap_orb.irp.f create mode 100644 src/casscf_cipsi/tot_en.irp.f diff --git a/src/casscf_cipsi/50.casscf.bats b/src/casscf_cipsi/50.casscf.bats new file mode 100644 index 00000000..a0db725d --- /dev/null +++ b/src/casscf_cipsi/50.casscf.bats @@ -0,0 +1,49 @@ +#!/usr/bin/env bats + +source $QP_ROOT/tests/bats/common.bats.sh +source $QP_ROOT/quantum_package.rc + + +function run_stoch() { + thresh=$2 + test_exe casscf || skip + qp set perturbation do_pt2 True + qp set determinants n_det_max $3 + qp set davidson threshold_davidson 1.e-10 + qp set davidson n_states_diag 4 + qp run casscf | tee casscf.out + energy1="$(ezfio get casscf energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)" + eq $energy1 $1 $thresh +} + +@test "F2" { # 18.0198s + rm -rf f2_casscf + qp_create_ezfio -b aug-cc-pvdz ../input/f2.zmt -o f2_casscf + qp set_file f2_casscf + qp run scf + qp set_mo_class --core="[1-6,8-9]" --act="[7,10]" --virt="[11-46]" + run_stoch -198.773366970 1.e-4 100000 +} + +@test "N2" { # 18.0198s + rm -rf n2_casscf + qp_create_ezfio -b aug-cc-pvdz ../input/n2.xyz -o n2_casscf + qp set_file n2_casscf + qp run scf + qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]" + run_stoch -109.0961643162 1.e-4 100000 +} + +@test "N2_stretched" { + rm -rf n2_stretched_casscf + qp_create_ezfio -b aug-cc-pvdz -m 7 ../input/n2_stretched.xyz -o n2_stretched_casscf + qp set_file n2_stretched_casscf + qp run scf | tee scf.out + qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]" + qp set electrons elec_alpha_num 7 + qp set electrons elec_beta_num 7 + run_stoch -108.7860471300 1.e-4 100000 +# + +} + diff --git a/src/casscf_cipsi/EZFIO.cfg b/src/casscf_cipsi/EZFIO.cfg new file mode 100644 index 00000000..2a1f1926 --- /dev/null +++ b/src/casscf_cipsi/EZFIO.cfg @@ -0,0 +1,75 @@ +[energy] +type: double precision +doc: Calculated Selected |FCI| energy +interface: ezfio +size: (determinants.n_states) + +[energy_pt2] +type: double precision +doc: Calculated |FCI| energy + |PT2| +interface: ezfio +size: (determinants.n_states) + +[state_following_casscf] +type: logical +doc: If |true|, the CASSCF will try to follow the guess CI vector and orbitals +interface: ezfio,provider,ocaml +default: False + + +[diag_hess_cas] +type: logical +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 +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 + +[criterion_casscf] +type: character*(32) +doc: choice of the criterion for the convergence of the casscf: can be energy or gradients or e_pt2 +interface: ezfio, provider, ocaml +default: e_pt2 + +[thresh_casscf] +type: Threshold +doc: Threshold on the convergence of the CASCF energy. +interface: ezfio,provider,ocaml +default: 1.e-06 + + +[pt2_min_casscf] +type: Threshold +doc: Minimum value of the pt2_max parameter for the CIPSI in the CASSCF iterations. +interface: ezfio,provider,ocaml +default: 1.e-04 + +[n_big_act_orb] +type: integer +doc: Number of active orbitals from which the active space is considered as large, and therefore pt2_min_casscf is activated. +interface: ezfio,provider,ocaml +default: 16 + +[adaptive_pt2_max] +type: logical +doc: If |true|, the pt2_max value in the CIPSI iterations will automatically adapt, otherwise it is fixed at the value given in the EZFIO folder +interface: ezfio,provider,ocaml +default: True diff --git a/src/casscf_cipsi/NEED b/src/casscf_cipsi/NEED new file mode 100644 index 00000000..dd91c7bd --- /dev/null +++ b/src/casscf_cipsi/NEED @@ -0,0 +1,5 @@ +cipsi +selectors_full +generators_cas +two_body_rdm +dav_general_mat diff --git a/src/casscf_cipsi/README.rst b/src/casscf_cipsi/README.rst new file mode 100644 index 00000000..08bfd95b --- /dev/null +++ b/src/casscf_cipsi/README.rst @@ -0,0 +1,5 @@ +====== +casscf +====== + +|CASSCF| program with the CIPSI algorithm. diff --git a/src/casscf_cipsi/bavard.irp.f b/src/casscf_cipsi/bavard.irp.f new file mode 100644 index 00000000..463c3ea4 --- /dev/null +++ b/src/casscf_cipsi/bavard.irp.f @@ -0,0 +1,6 @@ +! -*- F90 -*- +BEGIN_PROVIDER [logical, bavard] +! bavard=.true. + bavard=.false. +END_PROVIDER + diff --git a/src/casscf_cipsi/bielec.irp.f b/src/casscf_cipsi/bielec.irp.f new file mode 100644 index 00000000..0a44f994 --- /dev/null +++ b/src/casscf_cipsi/bielec.irp.f @@ -0,0 +1,155 @@ +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 + real*8 :: mo_two_e_integral + + bielec_PQxx(:,:,:,:) = 0.d0 + PROVIDE mo_two_e_integrals_in_map + + !$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_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 + + + !$OMP DO + do i=1,n_act_orb + ii=list_act(i) + i3=i+n_core_inact_orb + do j=i,n_act_orb + jj=list_act(j) + 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_inact_act_orb,n_core_inact_act_orb, mo_num)] + BEGIN_DOC + ! bielec_PxxQ : integral (px|xq) 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 + + 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)) + + !$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 + end do + end do + do j=1,n_act_orb + jj=list_act(j) + 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_inact_orb + do j=i,n_act_orb + jj=list_act(j) + 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 + + +BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] + BEGIN_DOC + ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active + ! index p runs over the whole basis, t,u,v only over the active orbitals + END_DOC + implicit none + integer :: i,j,k,p,t,u,v + double precision, external :: mo_two_e_integral + PROVIDE mo_two_e_integrals_in_map + + !$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) + 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 + !$OMP END PARALLEL DO + +END_PROVIDER diff --git a/src/casscf_cipsi/bielec_natorb.irp.f b/src/casscf_cipsi/bielec_natorb.irp.f new file mode 100644 index 00000000..9968530c --- /dev/null +++ b/src/casscf_cipsi/bielec_natorb.irp.f @@ -0,0 +1,369 @@ + 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 + double precision, allocatable :: f(:,:,:), d(:,:,:) + + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,p,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 + 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 + bielec_PQxx_no(list_act(p),j,k,l)=d(p,j,k) + end do + end do + + do j=1,mo_num + do p=1,n_act_orb + 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 + do j=1,mo_num + bielec_PQxx_no(j,list_act(p),k,l)=d(p,j,k) + end do + end do + end do + end do + !$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 + 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 + do k=1,mo_num + do j=1,mo_num + bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,p) + end do + end do + end do + end do + !$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 + 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 + do k=1,mo_num + do j=1,mo_num + bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,p) + 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_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 + double precision, allocatable :: f(:,:,:), d(:,:,:) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,p,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 + 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 + 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 + bielec_PxxQ_no(list_act(p),k,l,j)=d(p,k,l) + end do + end do + end do + end do + !$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 + 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 + bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(p,j,l) + end do + end do + end do + end do + !$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 + 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 + 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,p) + end do + end do + end do + end do + !$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 + 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 + 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,p) + end do + end do + end do + end do + !$OMP END DO NOWAIT + deallocate(f,d) + !$OMP END PARALLEL + +END_PROVIDER + + +BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)] + BEGIN_DOC + ! integrals (tu|vp) in the basis of natural MOs + ! index p runs over the whole basis, t,u,v only over the active orbitals + END_DOC + implicit none + integer :: i,j,k,l,t,u,p,q + double precision, allocatable :: f(:,:,:), d(:,:,:) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,p,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) + + 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 j=1,n_act_orb + do p=1,n_act_orb + f(p,j,k)=bielecCI_no(p,j,k,l) + end do + end do + end do + 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 j=1,n_act_orb + do p=1,n_act_orb + bielecCI_no(p,j,k,l)=d(p,j,k) + end do + end do + + do j=1,n_act_orb + do p=1,n_act_orb + f(p,j,k)=bielecCI_no(j,p,k,l) + end do + end do + end do + 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 p=1,n_act_orb + do j=1,n_act_orb + bielecCI_no(j,p,k,l)=d(p,j,k) + 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 + 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 + do k=1,n_act_orb + do j=1,n_act_orb + bielecCI_no(j,k,p,l)=d(j,k,p) + end do + end do + end do + end do + !$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 + 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 + do k=1,n_act_orb + do j=1,n_act_orb + bielecCI_no(j,k,l,list_act(p))=d(j,k,p) + end do + end do + end do + end do + !$OMP END DO + + deallocate(d,f) + !$OMP END PARALLEL + + +END_PROVIDER + diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f new file mode 100644 index 00000000..a2f3c5a7 --- /dev/null +++ b/src/casscf_cipsi/casscf.irp.f @@ -0,0 +1,110 @@ +program casscf + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + call reorder_orbitals_for_casscf +! no_vvvv_integrals = .True. +! touch no_vvvv_integrals + n_det_max_full = 500 + touch n_det_max_full + pt2_relative_error = 0.04 + touch pt2_relative_error +! call run_stochastic_cipsi + call run +end + +subroutine run + implicit none + double precision :: energy_old, energy, pt2_max_before, ept2_before,delta_E + 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 + ept2_before = 0.d0 + if(adaptive_pt2_max)then + pt2_max = 0.005 + SOFT_TOUCH pt2_max + endif + do while (.not.converged) + print*,'pt2_max = ',pt2_max + call run_stochastic_cipsi + energy_old = energy + energy = eone+etwo+ecore + pt2_max_before = pt2_max + + call write_time(6) + call write_int(6,iteration,'CAS-SCF iteration = ') + call write_double(6,energy,'CAS-SCF energy = ') + if(n_states == 1)then + double precision :: E_PT2, PT2 + call ezfio_get_casscf_energy_pt2(E_PT2) + call ezfio_get_casscf_energy(PT2) + PT2 -= E_PT2 + call write_double(6,E_PT2,'E + PT2 energy = ') + call write_double(6,PT2,' PT2 = ') + call write_double(6,pt2_max,' PT2_MAX = ') + endif + + print*,'' + 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 = ') + print*,'' + call write_double(6,energy_improvement, 'Predicted energy improvement = ') + + if(criterion_casscf == "energy")then + converged = dabs(energy_improvement) < thresh_scf + else if (criterion_casscf == "gradients")then + converged = norm_grad_vec2 < thresh_scf + else if (criterion_casscf == "e_pt2")then + delta_E = dabs(E_PT2 - ept2_before) + converged = dabs(delta_E) < thresh_casscf + endif + ept2_before = E_PT2 + if(adaptive_pt2_max)then + pt2_max = dabs(energy_improvement / (pt2_relative_error)) + pt2_max = min(pt2_max, pt2_max_before) + if(n_act_orb.ge.n_big_act_orb)then + pt2_max = max(pt2_max,pt2_min_casscf) + endif + endif + print*,'' + call write_double(6,pt2_max, 'PT2_MAX for next iteration = ') + + mo_coef = NewOrbs + mo_occ = occnum + call save_mos + if(.not.converged)then + iteration += 1 + if(norm_grad_vec2.gt.0.01d0)then + N_det = N_states + else + N_det = max(N_det/8 ,N_states) + endif + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + read_wf = .True. + call clear_mo_map + SOFT_TOUCH mo_coef N_det psi_det psi_coef + if(adaptive_pt2_max)then + SOFT_TOUCH pt2_max + endif + if(iteration .gt. 3)then + state_following_casscf = state_following_casscf_save + soft_touch state_following_casscf + endif + endif + + enddo + +end + + diff --git a/src/casscf_cipsi/class.irp.f b/src/casscf_cipsi/class.irp.f new file mode 100644 index 00000000..7360a661 --- /dev/null +++ b/src/casscf_cipsi/class.irp.f @@ -0,0 +1,12 @@ + BEGIN_PROVIDER [ logical, do_only_1h1p ] +&BEGIN_PROVIDER [ logical, do_only_cas ] +&BEGIN_PROVIDER [ logical, do_ddci ] + implicit none + BEGIN_DOC + ! In the CAS case, all those are always false except do_only_cas + END_DOC + do_only_cas = .True. + do_only_1h1p = .False. + do_ddci = .False. +END_PROVIDER + diff --git a/src/casscf_cipsi/dav_sx_mat.irp.f b/src/casscf_cipsi/dav_sx_mat.irp.f new file mode 100644 index 00000000..1e24f0e2 --- /dev/null +++ b/src/casscf_cipsi/dav_sx_mat.irp.f @@ -0,0 +1,45 @@ + + +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, min(N_st, N_st_diag_in) +! 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/src/casscf_cipsi/densities.irp.f b/src/casscf_cipsi/densities.irp.f new file mode 100644 index 00000000..bebcf5d7 --- /dev/null +++ b/src/casscf_cipsi/densities.irp.f @@ -0,0 +1,67 @@ +use bitmasks + +BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ] + implicit none + BEGIN_DOC + ! the first-order density matrix in the basis of the starting MOs. + ! matrix is state averaged. + END_DOC + integer :: t,u + + do u=1,n_act_orb + do t=1,n_act_orb + D0tu(t,u) = one_e_dm_mo_alpha_average( list_act(t), list_act(u) ) + & + one_e_dm_mo_beta_average ( list_act(t), list_act(u) ) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [real*8, P0tuvx, (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 ONLY IN THE RANGE OF ACTIVE MOS + ! The values 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> + ! + ! P0tuvx(p,q,r,s) = chemist notation : 1/2 <0|E_pq E_rs - delta_qr E_ps|0> + END_DOC + implicit none + integer :: t,u,v,x + integer :: tt,uu,vv,xx + integer :: 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 the 2 body RDM on the active part' + endif + + P0tuvx= 0.d0 + 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) = 0.5d0 * state_av_act_2_rdm_spin_trace_mo(t,v,u,x) + enddo + enddo + enddo + enddo + enddo + else + P0tuvx = P0tuvx_peter + endif + +END_PROVIDER diff --git a/src/casscf_cipsi/densities_peter.irp.f b/src/casscf_cipsi/densities_peter.irp.f new file mode 100644 index 00000000..ee7414da --- /dev/null +++ b/src/casscf_cipsi/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 diff --git a/src/casscf_cipsi/det_manip.irp.f b/src/casscf_cipsi/det_manip.irp.f new file mode 100644 index 00000000..d8c309a4 --- /dev/null +++ b/src/casscf_cipsi/det_manip.irp.f @@ -0,0 +1,125 @@ +use bitmasks + +subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, & + ispin,phase,ierr) + BEGIN_DOC + ! we create the mono-excitation, and determine, if possible, + ! the phase and the number in the list of determinants + END_DOC + implicit none + integer(bit_kind) :: key1(N_int,2),key2(N_int,2) + integer(bit_kind), allocatable :: keytmp(:,:) + integer :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin + real*8 :: phase + logical :: found + allocate(keytmp(N_int,2)) + + nu=-1 + phase=1.D0 + ierr=0 + call det_copy(key1,key2,N_int) + ! write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin + ! call print_det(key2,N_int) + call do_single_excitation(key2,ihole,ipart,ispin,ierr) + ! write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin + ! call print_det(key2,N_int) + ! write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr + if (ierr.eq.1) then + ! excitation is possible + ! get the phase + call get_single_excitation(key1,key2,exc,phase,N_int) + ! get the number in the list + found=.false. + nu=0 + + !TODO BOTTLENECK + do while (.not.found) + nu+=1 + if (nu.gt.N_det) then + ! the determinant is possible, but not in the list + found=.true. + nu=-1 + else + call det_extract(keytmp,nu,N_int) + integer :: i,ii + found=.true. + do ii=1,2 + do i=1,N_int + if (keytmp(i,ii).ne.key2(i,ii)) then + found=.false. + end if + end do + end do + end if + end do + end if + ! + ! we found the new string, the phase, and possibly the number in the list + ! +end subroutine do_signed_mono_excitation + +subroutine det_extract(key,nu,Nint) + BEGIN_DOC + ! extract a determinant from the list of determinants + END_DOC + implicit none + integer :: ispin,i,nu,Nint + integer(bit_kind) :: key(Nint,2) + do ispin=1,2 + do i=1,Nint + key(i,ispin)=psi_det(i,ispin,nu) + end do + end do +end subroutine det_extract + +subroutine det_copy(key1,key2,Nint) + use bitmasks ! you need to include the bitmasks_module.f90 features + BEGIN_DOC + ! copy a determinant from key1 to key2 + END_DOC + implicit none + integer :: ispin,i,Nint + integer(bit_kind) :: key1(Nint,2),key2(Nint,2) + do ispin=1,2 + do i=1,Nint + key2(i,ispin)=key1(i,ispin) + end do + end do +end subroutine det_copy + +subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2 & + ,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr) + BEGIN_DOC + ! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q) + ! we may create two determinants as result + ! + END_DOC + implicit none + integer(bit_kind) :: key_in(N_int,2),key_out1(N_int,2) + integer(bit_kind) :: key_out2(N_int,2) + integer :: ihole,ipart,ierr,jerr,nu1,nu2 + integer :: ispin + real*8 :: phase1,phase2 + + ! write(6,*) ' applying E_',ipart,ihole,' on determinant ' + ! call print_det(key_in,N_int) + + ! spin alpha + ispin=1 + call do_signed_mono_excitation(key_in,key_out1,nu1,ihole & + ,ipart,ispin,phase1,ierr) + ! if (ierr.eq.1) then + ! write(6,*) ' 1 result is ',nu1,phase1 + ! call print_det(key_out1,N_int) + ! end if + ! spin beta + ispin=2 + call do_signed_mono_excitation(key_in,key_out2,nu2,ihole & + ,ipart,ispin,phase2,jerr) + ! if (jerr.eq.1) then + ! write(6,*) ' 2 result is ',nu2,phase2 + ! call print_det(key_out2,N_int) + ! end if + +end subroutine do_spinfree_mono_excitation + diff --git a/src/casscf_cipsi/driver_optorb.irp.f b/src/casscf_cipsi/driver_optorb.irp.f new file mode 100644 index 00000000..2e3e02dc --- /dev/null +++ b/src/casscf_cipsi/driver_optorb.irp.f @@ -0,0 +1,3 @@ +subroutine driver_optorb + implicit none +end diff --git a/src/casscf_cipsi/get_energy.irp.f b/src/casscf_cipsi/get_energy.irp.f new file mode 100644 index 00000000..cfb26b59 --- /dev/null +++ b/src/casscf_cipsi/get_energy.irp.f @@ -0,0 +1,51 @@ +program print_2rdm + implicit none + BEGIN_DOC + ! get the active part of the bielectronic energy on a given wave function. + ! + ! useful to test the active part of the spin trace 2 rdms + END_DOC +!no_vvvv_integrals = .True. + read_wf = .True. +!touch read_wf no_vvvv_integrals +!call routine +!call routine_bis + call print_grad +end + +subroutine print_grad + implicit none + integer :: i + do i = 1, nMonoEx + if(dabs(gradvec2(i)).gt.1.d-5)then + print*,'' + print*,i,gradvec2(i),excit(:,i) + endif + enddo +end + +subroutine routine + integer :: i,j,k,l + integer :: ii,jj,kk,ll + double precision :: accu(4),twodm,thr,act_twodm2,integral,get_two_e_integral + thr = 1.d-10 + + + accu = 0.d0 + do ll = 1, n_act_orb + l = list_act(ll) + do kk = 1, n_act_orb + k = list_act(kk) + do jj = 1, n_act_orb + j = list_act(jj) + do ii = 1, n_act_orb + i = list_act(ii) + integral = get_two_e_integral(i,j,k,l,mo_integrals_map) + accu(1) += state_av_act_2_rdm_spin_trace_mo(ii,jj,kk,ll) * integral + enddo + enddo + enddo + enddo + print*,'accu = ',accu(1) + +end diff --git a/src/casscf_cipsi/grad_old.irp.f b/src/casscf_cipsi/grad_old.irp.f new file mode 100644 index 00000000..d60a60c8 --- /dev/null +++ b/src/casscf_cipsi/grad_old.irp.f @@ -0,0 +1,74 @@ + +BEGIN_PROVIDER [real*8, gradvec_old, (nMonoEx)] + BEGIN_DOC + ! calculate the orbital gradient by hand, i.e. for + ! each determinant I we determine the string E_pq |I> (alpha and beta + ! separately) and generate + ! sum_I c_I is then the pq component of the orbital + ! gradient + ! E_pq = a^+_pa_q + a^+_Pa_Q + END_DOC + implicit none + integer :: ii,tt,aa,indx,ihole,ipart,istate + real*8 :: res + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + call calc_grad_elem(ihole,ipart,res) + gradvec_old(indx)=res + end do + + real*8 :: norm_grad + norm_grad=0.d0 + do indx=1,nMonoEx + norm_grad+=gradvec_old(indx)*gradvec_old(indx) + end do + norm_grad=sqrt(norm_grad) + if (bavard) then + write(6,*) + write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad + write(6,*) + endif + + +END_PROVIDER + +subroutine calc_grad_elem(ihole,ipart,res) + BEGIN_DOC + ! eq 18 of Siegbahn et al, Physica Scripta 1980 + ! we calculate 2 , q=hole, p=particle + END_DOC + implicit none + integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) + real*8 :: i_H_psi_array(N_states),phase + allocate(det_mu(N_int,2)) + allocate(det_mu_ex(N_int,2)) + + res=0.D0 + + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation on it + call det_copy(det_mu,det_mu_ex,N_int) + call do_signed_mono_excitation(det_mu,det_mu_ex,nu & + ,ihole,ipart,ispin,phase,ierr) + if (ierr.eq.1) then + call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase + end do + end if + end do + end do + + ! state-averaged gradient + res*=2.D0/dble(N_states) + +end subroutine calc_grad_elem + diff --git a/src/casscf_cipsi/gradient.irp.f b/src/casscf_cipsi/gradient.irp.f new file mode 100644 index 00000000..a1c5e947 --- /dev/null +++ b/src/casscf_cipsi/gradient.irp.f @@ -0,0 +1,215 @@ +use bitmasks + +BEGIN_PROVIDER [ integer, nMonoEx ] + BEGIN_DOC + ! Number of single excitations + END_DOC + implicit none + 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, n_c_a_prov] +&BEGIN_PROVIDER [integer, n_c_v_prov] +&BEGIN_PROVIDER [integer, n_a_v_prov] + implicit none + n_c_a_prov = n_core_inact_orb * n_act_orb + n_c_v_prov = n_core_inact_orb * n_virt_orb + n_a_v_prov = n_act_orb * n_virt_orb + END_PROVIDER + + BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] +&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] +&BEGIN_PROVIDER [integer, list_idx_c_a, (3,n_c_a_prov) ] +&BEGIN_PROVIDER [integer, list_idx_c_v, (3,n_c_v_prov) ] +&BEGIN_PROVIDER [integer, list_idx_a_v, (3,n_a_v_prov) ] +&BEGIN_PROVIDER [integer, mat_idx_c_a, (n_core_inact_orb,n_act_orb) +&BEGIN_PROVIDER [integer, mat_idx_c_v, (n_core_inact_orb,n_virt_orb) +&BEGIN_PROVIDER [integer, mat_idx_a_v, (n_act_orb,n_virt_orb) + BEGIN_DOC + ! a list of the orbitals involved in the excitation + END_DOC + + implicit none + integer :: i,t,a,ii,tt,aa,indx,indx_tmp + indx=0 + indx_tmp = 0 + do ii=1,n_core_inact_orb + i=list_core_inact(ii) + do tt=1,n_act_orb + t=list_act(tt) + indx+=1 + excit(1,indx)=i + excit(2,indx)=t + excit_class(indx)='c-a' + indx_tmp += 1 + list_idx_c_a(1,indx_tmp) = indx + list_idx_c_a(2,indx_tmp) = ii + list_idx_c_a(3,indx_tmp) = tt + mat_idx_c_a(ii,tt) = indx + end do + end do + + indx_tmp = 0 + do ii=1,n_core_inact_orb + i=list_core_inact(ii) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=i + excit(2,indx)=a + excit_class(indx)='c-v' + indx_tmp += 1 + list_idx_c_v(1,indx_tmp) = indx + list_idx_c_v(2,indx_tmp) = ii + list_idx_c_v(3,indx_tmp) = aa + mat_idx_c_v(ii,aa) = indx + end do + end do + + indx_tmp = 0 + do tt=1,n_act_orb + t=list_act(tt) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=t + excit(2,indx)=a + excit_class(indx)='a-v' + indx_tmp += 1 + list_idx_a_v(1,indx_tmp) = indx + list_idx_a_v(2,indx_tmp) = tt + list_idx_a_v(3,indx_tmp) = aa + mat_idx_a_v(tt,aa) = indx + end do + end do + + if (bavard) then + write(6,*) ' Filled the table of the Monoexcitations ' + do indx=1,nMonoEx + write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' & + ,excit(2,indx),' ',excit_class(indx) + end do + end if + +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 + ! eqs 14 a,b,c + END_DOC + implicit none + integer :: i,t,a,indx + 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 + + do i=1,n_core_inact_orb + 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 + + do t=1,n_act_orb + 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 + + norm_grad_vec2=0.d0 + 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,*) + write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad_vec2 + write(6,*) + endif + +END_PROVIDER + +real*8 function gradvec_it(i,t) + BEGIN_DOC + ! the orbital gradient core/inactive -> active + ! we assume natural orbitals + END_DOC + implicit none + integer :: i,t + + integer :: ii,tt,v,vv,x,y + integer :: x3,y3 + + 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 ! active + vv=list_act(v) + do x=1,n_act_orb ! active + x3=x+n_core_inact_orb ! list_act(x) + do y=1,n_act_orb ! active + y3=y+n_core_inact_orb ! list_act(y) + ! Gamma(2) a a a a 1/r12 i a a a + gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3) + end do + end do + end do + gradvec_it*=2.D0 +end function gradvec_it + +real*8 function gradvec_ia(i,a) + BEGIN_DOC + ! the orbital gradient core/inactive -> virtual + END_DOC + implicit none + integer :: i,a,ii,aa + + ii=list_core_inact(i) + aa=list_virt(a) + gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii)) + gradvec_ia*=2.D0 + +end function gradvec_ia + +real*8 function gradvec_ta(t,a) + BEGIN_DOC + ! the orbital gradient active -> virtual + ! we assume natural orbitals + END_DOC + implicit none + integer :: t,a,tt,aa,v,vv,x,y + + tt=list_act(t) + aa=list_virt(a) + gradvec_ta=0.D0 + gradvec_ta+=occnum(tt)*Fipq(aa,tt) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa) + end do + end do + end do + gradvec_ta*=2.D0 + +end function gradvec_ta + diff --git a/src/casscf_cipsi/hessian.irp.f b/src/casscf_cipsi/hessian.irp.f new file mode 100644 index 00000000..458c6aa6 --- /dev/null +++ b/src/casscf_cipsi/hessian.irp.f @@ -0,0 +1,539 @@ +use bitmasks + +real*8 function hessmat_itju(i,t,j,u) + BEGIN_DOC + ! 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 + END_DOC + implicit none + integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj + real*8 :: term,t2 + + ii=list_core_inact(i) + tt=list_act(t) + if (i.eq.j) then + if (t.eq.u) then + ! diagonal element + term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt)) & + -2.D0*(Fipq(ii,ii)+Fapq(ii,ii)) + term+=2.D0*(3.D0*bielec_pxxq_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i)) + term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt) & + -bielec_pqxx_no(tt,tt,i,i)) + term-=occnum(tt)*Fipq(tt,tt) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(vv,xx,i,i) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & + bielec_pxxq_no(vv,i,i,xx)) + do y=1,n_act_orb + term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) + end do + end do + end do + else + ! it/iu, t != u + uu=list_act(u) + term=2.D0*(Fipq(tt,uu)+Fapq(tt,uu)) + term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & + -bielec_PQxx_no(tt,uu,i,j)) + term-=occnum(tt)*Fipq(uu,tt) + term-=(occnum(tt)+occnum(uu)) & + *(3.D0*bielec_PxxQ_no(tt,i,i,uu)-bielec_PQxx_no(uu,tt,i,i)) + do v=1,n_act_orb + vv=list_act(v) + ! term-=D0tu(u,v)*Fipq(tt,vv) ! published, but inverting t and u seems more correct + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,i) & + +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & + *bielec_pxxq_no(vv,i,i,xx)) + do y=1,n_act_orb + term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(u,v,y,xx) + end do + end do + end do + end if + else + ! it/ju + jj=list_core_inact(j) + uu=list_act(u) + if (t.eq.u) then + term=occnum(tt)*Fipq(ii,jj) + term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) + else + term=0.D0 + end if + term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & + -bielec_PQxx_no(tt,uu,i,j)) + term-=(occnum(tt)+occnum(uu))* & + (4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & + -bielec_PQxx_no(uu,tt,i,j)) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,j) & + +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & + *bielec_pxxq_no(vv,i,j,xx)) + end do + end do + end if + + term*=2.D0 + hessmat_itju=term + +end function hessmat_itju + +real*8 function hessmat_itja(i,t,j,a) + BEGIN_DOC + ! 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_inact(i) + tt=list_act(t) + 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)) + term-=occnum(tt)*(4.D0*bielec_pxxq_no(aa,j,i,tt) & + -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt)) + if (i.eq.j) then + term+=2.D0*(Fipq(aa,tt)+Fapq(aa,tt)) + term-=0.5D0*occnum(tt)*Fipq(aa,tt) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa) + end do + end do + end do + end if + term*=2.D0 + hessmat_itja=term + +end function hessmat_itja + +real*8 function hessmat_itua(i,t,u,a) + BEGIN_DOC + ! 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_inact(i) + tt=list_act(t) + t3=t+n_core_inact_orb + uu=list_act(u) + u3=u+n_core_inact_orb + aa=list_virt(a) + if (t.eq.u) then + term=-occnum(tt)*Fipq(aa,ii) + else + term=0.D0 + end if + term-=occnum(uu)*(bielec_pqxx_no(aa,ii,t3,u3)-4.D0*bielec_pqxx_no(aa,uu,t3,i)& + +bielec_pxxq_no(aa,t3,u3,ii)) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_inact_orb + do x=1,n_act_orb + integer :: x3 + xx=list_act(x) + 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)) + end do + end do + if (t.eq.u) then + term+=Fipq(aa,ii)+Fapq(aa,ii) + end if + term*=2.D0 + hessmat_itua=term + +end function hessmat_itua + +real*8 function hessmat_iajb(i,a,j,b) + BEGIN_DOC + ! 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_inact(i) + aa=list_virt(a) + if (i.eq.j) then + if (a.eq.b) then + ! ia/ia + term=2.D0*(Fipq(aa,aa)+Fapq(aa,aa)-Fipq(ii,ii)-Fapq(ii,ii)) + term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,aa)-bielec_pqxx_no(aa,aa,i,i)) + else + bb=list_virt(b) + ! ia/ib + term=2.D0*(Fipq(aa,bb)+Fapq(aa,bb)) + term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,bb)-bielec_pqxx_no(aa,bb,i,i)) + end if + else + ! ia/jb + 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)) + if (a.eq.b) then + term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) + end if + end if + term*=2.D0 + hessmat_iajb=term + +end function hessmat_iajb + +real*8 function hessmat_iatb(i,a,t,b) + BEGIN_DOC + ! 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_inact(i) + aa=list_virt(a) + tt=list_act(t) + bb=list_virt(b) + 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 + term-=Fipq(tt,ii)+Fapq(tt,ii) + term-=0.5D0*occnum(tt)*Fipq(tt,ii) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,ii) + end do + end do + end do + end if + term*=2.D0 + hessmat_iatb=term + +end function hessmat_iatb + +real*8 function hessmat_taub(t,a,u,b) + BEGIN_DOC + ! the orbital hessian for act->virt,act->virt + END_DOC + implicit none + integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y + integer :: v3,x3 + real*8 :: term,t1,t2,t3 + + tt=list_act(t) + aa=list_virt(a) + if (t == u) then + if (a == b) then + ! ta/ta + t1=occnum(tt)*Fipq(aa,aa) + t2=0.D0 + t3=0.D0 + t1-=occnum(tt)*Fipq(tt,tt) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_inact_orb + do x=1,n_act_orb + xx=list_act(x) + 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)) + do y=1,n_act_orb + t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) + end do + end do + end do + term=t1+t2+t3 + else + bb=list_virt(b) + ! ta/tb b/=a + term=occnum(tt)*Fipq(aa,bb) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_inact_orb + do x=1,n_act_orb + xx=list_act(x) + 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)) + end do + end do + end if + else + ! ta/ub t/=u + uu=list_act(u) + bb=list_virt(b) + term=0.D0 + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_inact_orb + do x=1,n_act_orb + xx=list_act(x) + 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)) + end do + end do + if (a.eq.b) then + term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu)) + do v=1,n_act_orb + do y=1,n_act_orb + do x=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu) + term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt) + end do + end do + end do + end if + + end if + + term*=2.D0 + hessmat_taub=term + +end function hessmat_taub + +BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] + BEGIN_DOC + ! the diagonal of the Hessian, needed for the Davidson procedure + END_DOC + implicit none + integer :: i,t,a,indx,indx_shift + real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessdiag,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & + !$OMP PRIVATE(i,indx,t,a,indx_shift) + + !$OMP DO + do i=1,n_core_inact_orb + do t=1,n_act_orb + indx = t + (i-1)*n_act_orb + hessdiag(indx)=hessmat_itju(i,t,i,t) + end do + end do + !$OMP END DO NOWAIT + + indx_shift = n_core_inact_orb*n_act_orb + !$OMP DO + do a=1,n_virt_orb + do i=1,n_core_inact_orb + indx = a + (i-1)*n_virt_orb + indx_shift + hessdiag(indx)=hessmat_iajb(i,a,i,a) + end do + end do + !$OMP END DO NOWAIT + + indx_shift += n_core_inact_orb*n_virt_orb + !$OMP DO + do a=1,n_virt_orb + do t=1,n_act_orb + indx = a + (t-1)*n_virt_orb + indx_shift + hessdiag(indx)=hessmat_taub(t,a,t,a) + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] + implicit none + integer :: i,j,t,u,a,b + integer :: indx,indx_tmp, jndx, jndx_tmp + integer :: ustart,bstart + real*8 :: hessmat_itju + real*8 :: hessmat_itja + real*8 :: hessmat_itua + real*8 :: hessmat_iajb + real*8 :: hessmat_iatb + real*8 :: hessmat_taub + ! c-a c-v a-v + ! c-a | X X X + ! c-v | X X + ! a-v | X + + provide mo_two_e_integrals_in_map + + hessmat = 0.d0 + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_a_prov,list_idx_c_a,n_core_inact_orb,n_act_orb,mat_idx_c_a) & + !$OMP PRIVATE(indx_tmp,indx,i,t,j,u,ustart,jndx) + + !$OMP DO +!!!! < Core-active| H |Core-active > + ! Core-active excitations + do indx_tmp = 1, n_c_a_prov + indx = list_idx_c_a(1,indx_tmp) + i = list_idx_c_a(2,indx_tmp) + t = list_idx_c_a(3,indx_tmp) + ! Core-active excitations + do j = 1, n_core_inact_orb + if (i.eq.j) then + ustart=t + else + ustart=1 + end if + do u=ustart,n_act_orb + jndx = mat_idx_c_a(j,u) + hessmat(jndx,indx) = hessmat_itju(i,t,j,u) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_a_prov,n_c_v_prov,list_idx_c_a,list_idx_c_v) & + !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,t,j,a,jndx) + + !$OMP DO +!!!! < Core-active| H |Core-VIRTUAL > + ! Core-active excitations + do indx_tmp = 1, n_c_a_prov + indx = list_idx_c_a(1,indx_tmp) + i = list_idx_c_a(2,indx_tmp) + t = list_idx_c_a(3,indx_tmp) + ! Core-VIRTUAL excitations + do jndx_tmp = 1, n_c_v_prov + jndx = list_idx_c_v(1,jndx_tmp) + j = list_idx_c_v(2,jndx_tmp) + a = list_idx_c_v(3,jndx_tmp) + hessmat(jndx,indx) = hessmat_itja(i,t,j,a) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_a_prov,n_a_v_prov,list_idx_c_a,list_idx_a_v) & + !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,t,u,a,jndx) + + !$OMP DO +!!!! < Core-active| H |ACTIVE-VIRTUAL > + ! Core-active excitations + do indx_tmp = 1, n_c_a_prov + indx = list_idx_c_a(1,indx_tmp) + i = list_idx_c_a(2,indx_tmp) + t = list_idx_c_a(3,indx_tmp) + ! ACTIVE-VIRTUAL excitations + do jndx_tmp = 1, n_a_v_prov + jndx = list_idx_a_v(1,jndx_tmp) + u = list_idx_a_v(2,jndx_tmp) + a = list_idx_a_v(3,jndx_tmp) + hessmat(jndx,indx) = hessmat_itua(i,t,u,a) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + + !$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) + ! 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 + + !$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) & + !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,a,t,b,jndx) + + !$OMP DO +!!!! < Core-VIRTUAL | H |Active-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) + ! Active-VIRTUAL excitations + do jndx_tmp = 1, n_a_v_prov + jndx = list_idx_a_v(1,jndx_tmp) + t = list_idx_a_v(2,jndx_tmp) + b = list_idx_a_v(3,jndx_tmp) + hessmat(jndx,indx) = hessmat_iatb(i,a,t,b) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_a_v_prov,list_idx_a_v,n_act_orb,n_virt_orb,mat_idx_a_v) & + !$OMP PRIVATE(indx_tmp,indx,t,a,u,b,bstart,jndx) + + !$OMP DO +!!!! < Active-VIRTUAL | H |Active-VIRTUAL > + ! Active-VIRTUAL excitations + do indx_tmp = 1, n_a_v_prov + indx = list_idx_a_v(1,indx_tmp) + t = list_idx_a_v(2,indx_tmp) + a = list_idx_a_v(3,indx_tmp) + ! Active-VIRTUAL excitations + do u=t,n_act_orb + if (t.eq.u) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + jndx = mat_idx_a_v(u,b) + hessmat(jndx,indx) = hessmat_taub(t,a,u,b) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + +END_PROVIDER diff --git a/src/casscf_cipsi/hessian_old.irp.f b/src/casscf_cipsi/hessian_old.irp.f new file mode 100644 index 00000000..d17f1f0a --- /dev/null +++ b/src/casscf_cipsi/hessian_old.irp.f @@ -0,0 +1,310 @@ + +use bitmasks +BEGIN_PROVIDER [real*8, hessmat_old, (nMonoEx,nMonoEx)] + BEGIN_DOC + ! calculate the orbital hessian 2 + ! + + by hand, + ! determinant per determinant, as for the gradient + ! + ! we assume that we have natural active orbitals + END_DOC + implicit none + integer :: indx,ihole,ipart + integer :: jndx,jhole,jpart + character*3 :: iexc,jexc + real*8 :: res + + if (bavard) then + write(6,*) ' providing Hessian matrix hessmat_old ' + write(6,*) ' nMonoEx = ',nMonoEx + endif + + do indx=1,nMonoEx + do jndx=1,nMonoEx + hessmat_old(indx,jndx)=0.D0 + end do + end do + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + iexc=excit_class(indx) + do jndx=indx,nMonoEx + jhole=excit(1,jndx) + jpart=excit(2,jndx) + jexc=excit_class(jndx) + call calc_hess_elem(ihole,ipart,jhole,jpart,res) + hessmat_old(indx,jndx)=res + hessmat_old(jndx,indx)=res + end do + end do + +END_PROVIDER + +subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) + BEGIN_DOC + ! eq 19 of Siegbahn et al, Physica Scripta 1980 + ! we calculate 2 + ! + + + ! average over all states is performed. + ! no transition between states. + END_DOC + implicit none + integer :: ihole,ipart,ispin,mu,istate + integer :: jhole,jpart,jspin + integer :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:) + integer(bit_kind), allocatable :: det_nu(:,:) + integer(bit_kind), allocatable :: det_mu_pq(:,:) + integer(bit_kind), allocatable :: det_mu_rs(:,:) + integer(bit_kind), allocatable :: det_nu_rs(:,:) + integer(bit_kind), allocatable :: det_mu_pqrs(:,:) + integer(bit_kind), allocatable :: det_mu_rspq(:,:) + real*8 :: i_H_psi_array(N_states),phase,phase2,phase3 + real*8 :: i_H_j_element + allocate(det_mu(N_int,2)) + allocate(det_nu(N_int,2)) + allocate(det_mu_pq(N_int,2)) + allocate(det_mu_rs(N_int,2)) + allocate(det_nu_rs(N_int,2)) + allocate(det_mu_pqrs(N_int,2)) + allocate(det_mu_rspq(N_int,2)) + integer :: mu_pq_possible + integer :: mu_rs_possible + integer :: nu_rs_possible + integer :: mu_pqrs_possible + integer :: mu_rspq_possible + + res=0.D0 + + ! the terms <0|E E H |0> + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation pq on it + call det_copy(det_mu,det_mu_pq,N_int) + call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq & + ,ihole,ipart,ispin,phase,mu_pq_possible) + if (mu_pq_possible.eq.1) then + ! possible, but not necessarily in the list + ! do the second excitation + do jspin=1,2 + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& + ,jhole,jpart,jspin,phase2,mu_pqrs_possible) + ! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 + end do + end if + ! try the de-excitation with opposite sign + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& + ,jpart,jhole,jspin,phase2,mu_pqrs_possible) + phase2=-phase2 + ! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 + end do + end if + end do + end if + ! exchange the notion of pq and rs + ! do the monoexcitation rs on the initial determinant + call det_copy(det_mu,det_mu_rs,N_int) + call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs & + ,jhole,jpart,ispin,phase2,mu_rs_possible) + if (mu_rs_possible.eq.1) then + ! do the second excitation + do jspin=1,2 + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& + ,ihole,ipart,jspin,phase3,mu_rspq_possible) + ! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 + end do + end if + ! we may try the de-excitation, with opposite sign + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& + ,ipart,ihole,jspin,phase3,mu_rspq_possible) + phase3=-phase3 + ! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 + end do + end if + end do + end if + ! + ! the operator E H E, we have to do a double loop over the determinants + ! we still have the determinant mu_pq and the phase in memory + if (mu_pq_possible.eq.1) then + do nu=1,N_det + call det_extract(det_nu,nu,N_int) + do jspin=1,2 + call det_copy(det_nu,det_nu_rs,N_int) + call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs & + ,jhole,jpart,jspin,phase2,nu_rs_possible) + ! excitation possible ? + if (nu_rs_possible.eq.1) then + call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element) + do istate=1,N_states + res+=2.D0*i_H_j_element*psi_coef(mu,istate) & + *psi_coef(nu,istate)*phase*phase2 + end do + end if + end do + end do + end if + end do + end do + + ! state-averaged Hessian + res*=1.D0/dble(N_states) + +end subroutine calc_hess_elem + +BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] + BEGIN_DOC + ! explicit hessian matrix from density matrices and integrals + ! of course, this will be used for a direct Davidson procedure later + ! we will not store the matrix in real life + ! formulas are broken down as functions for the 6 classes of matrix elements + ! + END_DOC + implicit none + integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart,indx_shift + + real*8 :: hessmat_itju + real*8 :: hessmat_itja + real*8 :: hessmat_itua + real*8 :: hessmat_iajb + real*8 :: hessmat_iatb + real*8 :: hessmat_taub + + if (bavard) then + write(6,*) ' providing Hessian matrix hessmat_peter ' + write(6,*) ' nMonoEx = ',nMonoEx + endif + provide mo_two_e_integrals_in_map + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat_peter,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & + !$OMP PRIVATE(i,indx,jndx,j,ustart,t,u,a,bstart,indx_shift) + + !$OMP DO + ! (DOUBLY OCCUPIED ---> ACT ) + do i=1,n_core_inact_orb + do t=1,n_act_orb + indx = t + (i-1)*n_act_orb + jndx=indx + ! (DOUBLY OCCUPIED ---> ACT ) + do j=i,n_core_inact_orb + if (i.eq.j) then + ustart=t + else + ustart=1 + end if + do u=ustart,n_act_orb + hessmat_peter(jndx,indx)=hessmat_itju(i,t,j,u) + jndx+=1 + end do + end do + ! (DOUBLY OCCUPIED ---> VIRTUAL) + do j=1,n_core_inact_orb + do a=1,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_itja(i,t,j,a) + jndx+=1 + end do + end do + ! (ACTIVE ---> VIRTUAL) + do u=1,n_act_orb + do a=1,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_itua(i,t,u,a) + jndx+=1 + end do + end do + end do + end do + !$OMP END DO NOWAIT + + indx_shift = n_core_inact_orb*n_act_orb + !$OMP DO + ! (DOUBLY OCCUPIED ---> VIRTUAL) + do a=1,n_virt_orb + do i=1,n_core_inact_orb + indx = a + (i-1)*n_virt_orb + indx_shift + jndx=indx + ! (DOUBLY OCCUPIED ---> VIRTUAL) + do j=i,n_core_inact_orb + if (i.eq.j) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_iajb(i,a,j,b) + jndx+=1 + end do + end do + ! (ACT ---> VIRTUAL) + do t=1,n_act_orb + do b=1,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_iatb(i,a,t,b) + jndx+=1 + end do + end do + end do + end do + !$OMP END DO NOWAIT + + indx_shift += n_core_inact_orb*n_virt_orb + !$OMP DO + ! (ACT ---> VIRTUAL) + do a=1,n_virt_orb + do t=1,n_act_orb + indx = a + (t-1)*n_virt_orb + indx_shift + jndx=indx + ! (ACT ---> VIRTUAL) + do u=t,n_act_orb + if (t.eq.u) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_taub(t,a,u,b) + jndx+=1 + end do + end do + end do + end do + !$OMP END DO + + !$OMP END PARALLEL + + do jndx=1,nMonoEx + do indx=1,jndx-1 + hessmat_peter(indx,jndx) = hessmat_peter(jndx,indx) + enddo + enddo + + +END_PROVIDER + diff --git a/src/casscf_cipsi/mcscf_fock.irp.f b/src/casscf_cipsi/mcscf_fock.irp.f new file mode 100644 index 00000000..e4568405 --- /dev/null +++ b/src/casscf_cipsi/mcscf_fock.irp.f @@ -0,0 +1,80 @@ +BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] + BEGIN_DOC + ! the inactive Fock matrix, in molecular orbitals + END_DOC + implicit none + integer :: p,q,k,kk,t,tt,u,uu + + do q=1,mo_num + do p=1,mo_num + Fipq(p,q)=one_ints_no(p,q) + end do + end do + + ! the inactive Fock matrix + 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) + end do + end do + end do + + if (bavard) then + integer :: i + write(6,*) + write(6,*) ' the diagonal of the inactive effective Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) + write(6,*) + end if + + +END_PROVIDER + + +BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] + BEGIN_DOC + ! the active active Fock matrix, in molecular orbitals + ! we create them in MOs, quite expensive + ! + ! for an implementation in AOs we need first the natural orbitals + ! for forming an active density matrix in AOs + ! + END_DOC + implicit none + integer :: p,q,k,kk,t,tt,u,uu + + Fapq = 0.d0 + + ! the active Fock matrix, D0tu is diagonal + do t=1,n_act_orb + tt=list_act(t) + do q=1,mo_num + do p=1,mo_num + Fapq(p,q)+=occnum(tt) & + *(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q)) + end do + end do + end do + + if (bavard) then + integer :: i + write(6,*) + write(6,*) ' the effective Fock matrix over MOs' + write(6,*) + + write(6,*) + write(6,*) ' the diagonal of the inactive effective Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) + write(6,*) + write(6,*) + write(6,*) ' the diagonal of the active Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num) + write(6,*) + end if + + +END_PROVIDER + + diff --git a/src/casscf_cipsi/natorb.irp.f b/src/casscf_cipsi/natorb.irp.f new file mode 100644 index 00000000..9ce90304 --- /dev/null +++ b/src/casscf_cipsi/natorb.irp.f @@ -0,0 +1,231 @@ + BEGIN_PROVIDER [real*8, occnum, (mo_num)] + implicit none + BEGIN_DOC + ! MO occupation numbers + END_DOC + + integer :: i + occnum=0.D0 + do i=1,n_core_inact_orb + occnum(list_core_inact(i))=2.D0 + end do + + do i=1,n_act_orb + occnum(list_act(i))=occ_act(i) + end do + + if (bavard) then + write(6,*) ' occupation numbers ' + do i=1,mo_num + write(6,*) i,occnum(i) + end do + endif + +END_PROVIDER + + + BEGIN_PROVIDER [ real*8, natorbsCI, (n_act_orb,n_act_orb) ] +&BEGIN_PROVIDER [ real*8, occ_act, (n_act_orb) ] + implicit none + BEGIN_DOC + ! Natural orbitals of CI + END_DOC + integer :: i, j + double precision :: Vt(n_act_orb,n_act_orb) + +! call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb) + call svd(D0tu, size(D0tu,1), natorbsCI,size(natorbsCI,1), occ_act, Vt, size(Vt,1),n_act_orb,n_act_orb) + + if (bavard) then + write(6,*) ' found occupation numbers as ' + do i=1,n_act_orb + write(6,*) i,occ_act(i) + end do + + integer :: nmx + real*8 :: xmx + do i=1,n_act_orb + ! largest element of the eigenvector should be positive + xmx=0.D0 + nmx=0 + do j=1,n_act_orb + if (abs(natOrbsCI(j,i)).gt.xmx) then + nmx=j + xmx=abs(natOrbsCI(j,i)) + end if + end do + xmx=sign(1.D0,natOrbsCI(nmx,i)) + do j=1,n_act_orb + natOrbsCI(j,i)*=xmx + end do + + write(6,*) ' Eigenvector No ',i + write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb) + end do + end if + +END_PROVIDER + + +BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC + ! 4-index transformation of 2part matrices + END_DOC + integer :: i,j,k,l,p,q + real*8 :: d(n_act_orb) + + ! index per index + ! first quarter + P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:) + + 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 + end do + do p=1,n_act_orb + do q=1,n_act_orb + d(p)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(p,j,k,l)=d(p) + end do + end do + end do + end do + ! 2nd 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 + end do + do p=1,n_act_orb + do q=1,n_act_orb + d(p)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(j,p,k,l)=d(p) + end do + end do + end do + end do + ! 3rd 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 + end do + do p=1,n_act_orb + do q=1,n_act_orb + d(p)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(j,k,p,l)=d(p) + 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 + end do + do p=1,n_act_orb + do q=1,n_act_orb + d(p)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(j,k,l,p)=d(p) + end do + end do + end do + end do + +END_PROVIDER + + + +BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)] + implicit none + BEGIN_DOC + ! Transformed one-e integrals + END_DOC + integer :: i,j, p, q + real*8 :: d(n_act_orb) + one_ints_no(:,:)=mo_one_e_integrals(:,:) + + ! 1st half-trf + do j=1,mo_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + do q=1,n_act_orb + d(p)+=one_ints_no(list_act(q),j)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + one_ints_no(list_act(p),j)=d(p) + end do + end do + + ! 2nd half-trf + do j=1,mo_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + do q=1,n_act_orb + d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + one_ints_no(j,list_act(p))=d(p) + end do + end do +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, NatOrbsCI_mos, (mo_num, mo_num) ] + implicit none + BEGIN_DOC + ! Rotation matrix from current MOs to the CI natural MOs + END_DOC + integer :: p,q + + NatOrbsCI_mos(:,:) = 0.d0 + + do q = 1,mo_num + NatOrbsCI_mos(q,q) = 1.d0 + enddo + + do q = 1,n_act_orb + do p = 1,n_act_orb + NatOrbsCI_mos(list_act(p),list_act(q)) = natorbsCI(p,q) + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)] + implicit none + BEGIN_DOC +! FCI natural orbitals + END_DOC + + call dgemm('N','N', ao_num,mo_num,mo_num,1.d0, & + mo_coef, size(mo_coef,1), & + NatOrbsCI_mos, size(NatOrbsCI_mos,1), 0.d0, & + NatOrbsFCI, size(NatOrbsFCI,1)) +END_PROVIDER + diff --git a/src/casscf_cipsi/neworbs.irp.f b/src/casscf_cipsi/neworbs.irp.f new file mode 100644 index 00000000..a7cebbb2 --- /dev/null +++ b/src/casscf_cipsi/neworbs.irp.f @@ -0,0 +1,253 @@ + BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] +&BEGIN_PROVIDER [integer, n_guess_sx_mat ] + implicit none + BEGIN_DOC + ! Single-excitation matrix + END_DOC + + integer :: i,j + + do i=1,nMonoEx+1 + do j=1,nMonoEx+1 + SXmatrix(i,j)=0.D0 + end do + end do + + do i=1,nMonoEx + SXmatrix(1,i+1)=gradvec2(i) + SXmatrix(1+i,1)=gradvec2(i) + end do + if(diag_hess_cas)then + do i = 1, nMonoEx + SXmatrix(i+1,i+1) = hessdiag(i) + enddo + else + do i=1,nMonoEx + do j=1,nMonoEx + SXmatrix(i+1,j+1)=hessmat(i,j) + SXmatrix(j+1,i+1)=hessmat(i,j) + end do + end do + endif + + 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)] +&BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)] + implicit none + BEGIN_DOC + ! Eigenvectors/eigenvalues of the single-excitation matrix + END_DOC + 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 eigenvalues ' + write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1) + 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 + write(6,*) + write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1) + endif +END_PROVIDER + + BEGIN_PROVIDER [real*8, energy_improvement] + implicit none + if(state_following_casscf)then + energy_improvement = SXeigenval(best_vector_ovrlp_casscf) + else + energy_improvement = SXeigenval(1) + endif + END_PROVIDER + + + + BEGIN_PROVIDER [ integer, best_vector_ovrlp_casscf ] +&BEGIN_PROVIDER [ double precision, best_overlap_casscf ] + implicit none + integer :: i + double precision :: c0 + best_overlap_casscf = 0.D0 + best_vector_ovrlp_casscf = -1000 + do i=1,nMonoEx+1 + if (SXeigenval(i).lt.0.D0) then + 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 + end do + if(best_vector_ovrlp_casscf.lt.0)then + best_vector_ovrlp_casscf = minloc(SXeigenval,nMonoEx+1) + endif + 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,*) ' weight of the 1st element ',c0 + endif + END_PROVIDER + + BEGIN_PROVIDER [double precision, SXvector, (nMonoEx+1)] + implicit none + BEGIN_DOC + ! Best eigenvector of the single-excitation matrix + END_DOC + integer :: i + double precision :: c0 + c0=SXeigenvec(1,best_vector_ovrlp_casscf) + do i=1,nMonoEx+1 + SXvector(i)=SXeigenvec(i,best_vector_ovrlp_casscf)/c0 + end do + END_PROVIDER + + +BEGIN_PROVIDER [double precision, NewOrbs, (ao_num,mo_num) ] + implicit none + BEGIN_DOC + ! Updated orbitals + END_DOC + 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)) + + level_shift_casscf *= 0.5D0 + level_shift_casscf = max(level_shift_casscf,0.002d0) + !touch level_shift_casscf + else + if(best_vector_ovrlp_casscf.ne.1.and.n_orb_swap.ne.0)then + print*,'Taking the lowest root for the CASSCF' + print*,'!!! SWAPPING MOS !!!!!!' + level_shift_casscf *= 2.D0 + level_shift_casscf = min(level_shift_casscf,0.5d0) + print*,'level_shift_casscf = ',level_shift_casscf + NewOrbs = switch_mo_coef + !mo_coef = switch_mo_coef + !soft_touch mo_coef + !call save_mos_no_occ + !stop + else + level_shift_casscf *= 0.5D0 + level_shift_casscf = max(level_shift_casscf,0.002d0) + !touch level_shift_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)) + endif + endif + +END_PROVIDER + +BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! Orbital rotation matrix + END_DOC + integer :: i,j,indx,k,iter,t,a,ii,tt,aa + logical :: converged + + real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num) + real*8 :: Tmat(mo_num,mo_num) + real*8 :: f + + ! the orbital rotation matrix T + Tmat(:,:)=0.D0 + indx=1 + do i=1,n_core_inact_orb + ii=list_core_inact(i) + do t=1,n_act_orb + tt=list_act(t) + indx+=1 + Tmat(ii,tt)= SXvector(indx) + Tmat(tt,ii)=-SXvector(indx) + end do + end do + do i=1,n_core_inact_orb + ii=list_core_inact(i) + do a=1,n_virt_orb + aa=list_virt(a) + indx+=1 + Tmat(ii,aa)= SXvector(indx) + Tmat(aa,ii)=-SXvector(indx) + end do + end do + do t=1,n_act_orb + tt=list_act(t) + do a=1,n_virt_orb + aa=list_virt(a) + indx+=1 + Tmat(tt,aa)= SXvector(indx) + Tmat(aa,tt)=-SXvector(indx) + end do + end do + + ! Form the exponential + + Tpotmat(:,:)=0.D0 + Umat(:,:) =0.D0 + do i=1,mo_num + Tpotmat(i,i)=1.D0 + Umat(i,i) =1.d0 + end do + iter=0 + converged=.false. + do while (.not.converged) + iter+=1 + f = 1.d0 / dble(iter) + Tpotmat2(:,:) = Tpotmat(:,:) * f + call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, & + Tpotmat2, size(Tpotmat2,1), & + Tmat, size(Tmat,1), 0.d0, & + Tpotmat, size(Tpotmat,1)) + Umat(:,:) = Umat(:,:) + Tpotmat(:,:) + + converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30) + end do +END_PROVIDER + + + diff --git a/src/casscf_cipsi/reorder_orb.irp.f b/src/casscf_cipsi/reorder_orb.irp.f new file mode 100644 index 00000000..3cb90522 --- /dev/null +++ b/src/casscf_cipsi/reorder_orb.irp.f @@ -0,0 +1,70 @@ +subroutine reorder_orbitals_for_casscf + implicit none + BEGIN_DOC +! routine that reorders the orbitals of the CASSCF in terms block of core, active and virtual + END_DOC + integer :: i,j,iorb + integer, allocatable :: iorder(:),array(:) + allocate(iorder(mo_num),array(mo_num)) + do i = 1, n_core_orb + iorb = list_core(i) + array(iorb) = i + enddo + + do i = 1, n_inact_orb + iorb = list_inact(i) + array(iorb) = mo_num + i + enddo + + do i = 1, n_act_orb + iorb = list_act(i) + array(iorb) = 2 * mo_num + i + enddo + + do i = 1, n_virt_orb + iorb = list_virt(i) + array(iorb) = 3 * mo_num + i + enddo + + do i = 1, mo_num + iorder(i) = i + enddo + call isort(array,iorder,mo_num) + double precision, allocatable :: mo_coef_new(:,:) + allocate(mo_coef_new(ao_num,mo_num)) + do i = 1, mo_num + mo_coef_new(:,i) = mo_coef(:,iorder(i)) + enddo + mo_coef = mo_coef_new + touch mo_coef + + list_core_reverse = 0 + do i = 1, n_core_orb + list_core(i) = i + list_core_reverse(i) = i + mo_class(i) = "Core" + enddo + + list_inact_reverse = 0 + do i = 1, n_inact_orb + list_inact(i) = i + n_core_orb + list_inact_reverse(i+n_core_orb) = i + mo_class(i+n_core_orb) = "Inactive" + enddo + + list_act_reverse = 0 + do i = 1, n_act_orb + list_act(i) = n_core_inact_orb + i + list_act_reverse(n_core_inact_orb + i) = i + mo_class(n_core_inact_orb + i) = "Active" + enddo + + list_virt_reverse = 0 + do i = 1, n_virt_orb + list_virt(i) = n_core_inact_orb + n_act_orb + i + list_virt_reverse(n_core_inact_orb + n_act_orb + i) = i + mo_class(n_core_inact_orb + n_act_orb + i) = "Virtual" + enddo + touch list_core_reverse list_core list_inact list_inact_reverse list_act list_act_reverse list_virt list_virt_reverse + +end diff --git a/src/casscf_cipsi/save_energy.irp.f b/src/casscf_cipsi/save_energy.irp.f new file mode 100644 index 00000000..8729c5af --- /dev/null +++ b/src/casscf_cipsi/save_energy.irp.f @@ -0,0 +1,9 @@ +subroutine save_energy(E,pt2) + implicit none + BEGIN_DOC +! Saves the energy in |EZFIO|. + END_DOC + double precision, intent(in) :: E(N_states), pt2(N_states) + call ezfio_set_casscf_energy(E(1:N_states)) + call ezfio_set_casscf_energy_pt2(E(1:N_states)+pt2(1:N_states)) +end diff --git a/src/casscf_cipsi/superci_dm.irp.f b/src/casscf_cipsi/superci_dm.irp.f new file mode 100644 index 00000000..ee831c35 --- /dev/null +++ b/src/casscf_cipsi/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_cipsi/swap_orb.irp.f b/src/casscf_cipsi/swap_orb.irp.f new file mode 100644 index 00000000..49af207c --- /dev/null +++ b/src/casscf_cipsi/swap_orb.irp.f @@ -0,0 +1,132 @@ + BEGIN_PROVIDER [double precision, SXvector_lowest, (nMonoEx)] + implicit none + integer :: i + do i=2,nMonoEx+1 + SXvector_lowest(i-1)=SXeigenvec(i,1) + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, thresh_overlap_switch] + implicit none + thresh_overlap_switch = 0.5d0 + END_PROVIDER + + BEGIN_PROVIDER [integer, max_overlap, (nMonoEx)] +&BEGIN_PROVIDER [integer, n_max_overlap] +&BEGIN_PROVIDER [integer, dim_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 + iorder(i) = i + vec_tmp(i) = -dabs(SXvector_lowest(i)) + 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 + n_max_overlap += 1 + max_overlap(n_max_overlap) = iorder(i) + endif + enddo + dim_n_max_overlap = max(1,n_max_overlap) + END_PROVIDER + + BEGIN_PROVIDER [integer, orb_swap, (2,dim_n_max_overlap)] +&BEGIN_PROVIDER [integer, index_orb_swap, (dim_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" .and.hessmat(imono,imono).gt.0.d0)then ! core --> active rotation + n_orb_swap += 1 + orb_swap(1,n_orb_swap) = iorb ! core + orb_swap(2,n_orb_swap) = jorb ! active + index_orb_swap(n_orb_swap) = imono + else if (excit_class(imono) == "a-v" .and.hessmat(imono,imono).gt.0.d0)then ! active --> virtual rotation + n_orb_swap += 1 + orb_swap(1,n_orb_swap) = jorb ! virtual + orb_swap(2,n_orb_swap) = iorb ! active + index_orb_swap(n_orb_swap) = imono + endif + enddo + + integer,allocatable :: orb_swap_tmp(:,:) + allocate(orb_swap_tmp(2,dim_n_max_overlap)) + do i = 1, n_orb_swap + orb_swap_tmp(1,i) = orb_swap(1,i) + orb_swap_tmp(2,i) = orb_swap(2,i) + enddo + + 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)) + integer, allocatable :: index_orb_swap_tmp(:) + allocate(index_orb_swap_tmp(dim_n_max_overlap)) + index_orb_swap_tmp = index_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 + index_orb_swap(n_orb_swap) = index_orb_swap_tmp(i) + orb_swap(1,n_orb_swap) = orb_swap_tmp(1,i) + orb_swap(2,n_orb_swap) = orb_swap_tmp(2,i) + endif + enddo + + if(n_orb_swap.gt.0)then + print*,'n_orb_swap = ',n_orb_swap + endif + do i = 1, n_orb_swap + print*,'imono = ',index_orb_swap(i) + 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/casscf_cipsi/tot_en.irp.f b/src/casscf_cipsi/tot_en.irp.f new file mode 100644 index 00000000..1d70e087 --- /dev/null +++ b/src/casscf_cipsi/tot_en.irp.f @@ -0,0 +1,101 @@ + BEGIN_PROVIDER [real*8, etwo] +&BEGIN_PROVIDER [real*8, eone] +&BEGIN_PROVIDER [real*8, eone_bis] +&BEGIN_PROVIDER [real*8, etwo_bis] +&BEGIN_PROVIDER [real*8, etwo_ter] +&BEGIN_PROVIDER [real*8, ecore] +&BEGIN_PROVIDER [real*8, ecore_bis] + implicit none + integer :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3 + real*8 :: e_one_all,e_two_all + e_one_all=0.D0 + e_two_all=0.D0 + 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_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_inact_orb + do u=1,n_act_orb + uu=list_act(u) + 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 + end do + end do + do t=1,n_act_orb + tt=list_act(t) + do u=1,n_act_orb + 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_inact_orb + do x=1,n_act_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 + end do + end do + ecore =nuclear_repulsion + ecore_bis=nuclear_repulsion + 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_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 + end do + eone =0.D0 + eone_bis=0.D0 + etwo =0.D0 + etwo_bis=0.D0 + etwo_ter=0.D0 + do t=1,n_act_orb + tt=list_act(t) + t3=t+n_core_inact_orb + do u=1,n_act_orb + uu=list_act(u) + 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_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) & + -bielec_PxxQ(tt,i,i,uu)) + end do + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_inact_orb + do x=1,n_act_orb + xx=list_act(x) + 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) + h3=bielecCI(t,u,v,xx) + etwo +=P0tuvx(t,u,v,x)*h1 + etwo_bis+=P0tuvx(t,u,v,x)*h2 + etwo_ter+=P0tuvx(t,u,v,x)*h3 + if ((h1.ne.h2).or.(h1.ne.h3)) then + write(6,9901) t,u,v,x,h1,h2,h3 + 9901 format('aie: ',4I4,3E20.12) + end if + end do + end do + end do + end do + +END_PROVIDER + + From 1e0e06d9cd705553501af87b417e719760b918da Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 19 Jun 2023 15:26:14 +0200 Subject: [PATCH 03/21] fixed bug in bi_ort_ints/three_body_ijmk.irp.f: deallocate(tmp1) is wrong --- src/bi_ort_ints/three_body_ijmk.irp.f | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/bi_ort_ints/three_body_ijmk.irp.f b/src/bi_ort_ints/three_body_ijmk.irp.f index ee7e88ef..4a99fb1b 100644 --- a/src/bi_ort_ints/three_body_ijmk.irp.f +++ b/src/bi_ort_ints/three_body_ijmk.irp.f @@ -97,7 +97,7 @@ , tmp3, 3*n_points_final_grid, tmp1, 3*n_points_final_grid & , 0.d0, tmp_4d, mo_num*mo_num) - deallocate(tmp1) +! deallocate(tmp1) !$OMP PARALLEL DO PRIVATE(i,j,k,m) do i = 1, mo_num @@ -112,7 +112,7 @@ !$OMP END PARALLEL DO - +! allocate(tmp1(n_points_final_grid, 2, mo_num, mo_num)) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i, l, ipoint) & @@ -195,7 +195,7 @@ , tmp1, 3*n_points_final_grid, int2_grad1_u12_bimo_t, 3*n_points_final_grid & , 0.d0, tmp_4d, mo_num*mo_num) - deallocate(tmp1) +! deallocate(tmp1) !$OMP PARALLEL DO PRIVATE(i,j,k,m) do i = 1, mo_num From c3d257c7aceff4c8641fc751cf9996783c1ce96f Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 22 Jun 2023 11:58:32 +0200 Subject: [PATCH 04/21] added routines to rotate orbitals without touching core orbitals --- src/tc_bi_ortho/tc_natorb.irp.f | 13 +- src/tc_keywords/EZFIO.cfg | 6 + src/tc_scf/routines_rotates.irp.f | 16 +- src/utils/block_diag_degen_core.irp.f | 244 ++++++++++++++++++++++++++ 4 files changed, 271 insertions(+), 8 deletions(-) create mode 100644 src/utils/block_diag_degen_core.irp.f diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f index b7e5ae81..1b5a66f3 100644 --- a/src/tc_bi_ortho/tc_natorb.irp.f +++ b/src/tc_bi_ortho/tc_natorb.irp.f @@ -23,7 +23,7 @@ dm_tmp(1:mo_num,1:mo_num) = -tc_transition_matrix_mo(1:mo_num,1:mo_num,1,1) - print *, ' dm_tmp' + print *, ' Transition density matrix ' do i = 1, mo_num fock_diag(i) = fock_matrix_tc_mo_tot(i,i) write(*, '(100(F16.10,X))') -dm_tmp(:,i) @@ -32,8 +32,15 @@ thr_d = 1.d-6 thr_nd = 1.d-6 thr_deg = 1.d-3 - call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & - , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval) + if(n_core_orb.ne.0)then +! print*,'core orbitals' +! pause + call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & + , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval) + else + call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & + , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval) + endif ! call non_hrmt_bieig( mo_num, dm_tmp& ! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& ! , mo_num, natorb_tc_eigval ) diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index a69f5bac..f984d53a 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -208,6 +208,12 @@ doc: Threshold to determine if diagonal elements of the bi-orthogonal condition interface: ezfio,provider,ocaml default: 1.e-6 +[thresh_lr_angle] +type: double precision +doc: Maximum value of the angle between the couple of left and right orbital for the rotations +interface: ezfio,provider,ocaml +default: 20.0 + [thresh_biorthog_nondiag] type: Threshold doc: Threshold to determine if non-diagonal elements of L.T x R are close enouph to 0 diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f index 755c35b9..588382b5 100644 --- a/src/tc_scf/routines_rotates.irp.f +++ b/src/tc_scf/routines_rotates.irp.f @@ -140,7 +140,11 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) ! compute the overlap between the left and rescaled right call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat) ! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list) - call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list) + if(n_core_orb.ne.0)then + call give_degen_full_listcore(fock_diag, mo_num, list_core, n_core_orb, thr_deg, list_degen, n_degen_list) + else + call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list) + endif print *, ' fock_matrix_mo' do i = 1, mo_num print *, i, fock_diag(i), angle_left_right(i) @@ -152,6 +156,8 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) ! n_degen = ilast - ifirst +1 n_degen = list_degen(i,0) + if(n_degen .ge. 1000)n_degen = 1 ! convention for core orbitals + if(n_degen .eq. 1) cycle allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen)) @@ -279,7 +285,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) allocate(new_angles(mo_num)) new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num)) max_angle = maxval(new_angles) - good_angles = max_angle.lt.45.d0 + good_angles = max_angle.lt.thresh_lr_angle print *, ' max_angle = ', max_angle deallocate(new_angles) @@ -397,11 +403,11 @@ subroutine print_energy_and_mos(good_angles) print *, ' TC SCF energy gradient = ', grad_non_hermit print *, ' Max angle Left/right = ', max_angle_left_right - if(max_angle_left_right .lt. 45.d0) then + if(max_angle_left_right .lt. thresh_lr_angle) then print *, ' Maximum angle BELOW 45 degrees, everthing is OK !' good_angles = .true. - else if(max_angle_left_right .gt. 45.d0 .and. max_angle_left_right .lt. 75.d0) then - print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...' + else if(max_angle_left_right .gt. thresh_lr_angle .and. max_angle_left_right .lt. 75.d0) then + print *, ' Maximum angle between thresh_lr_angle and 75 degrees, this is not the best for TC-CI calculations ...' good_angles = .false. else if(max_angle_left_right .gt. 75.d0) then print *, ' Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...' diff --git a/src/utils/block_diag_degen_core.irp.f b/src/utils/block_diag_degen_core.irp.f new file mode 100644 index 00000000..5d46bd87 --- /dev/null +++ b/src/utils/block_diag_degen_core.irp.f @@ -0,0 +1,244 @@ + +subroutine diag_mat_per_fock_degen_core(fock_diag, mat_ref, listcore,ncore, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval) + + + BEGIN_DOC + ! + ! subroutine that diagonalizes a matrix mat_ref BY BLOCK + ! + ! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag" + ! + ! the elements of listcore are untouched + ! + ! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together + ! + ! : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together + ! + ! : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together + ! + ! etc... the advantage is to guarentee no spurious mixing because of numerical problems. + ! + END_DOC + + implicit none + integer, intent(in) :: n,ncore, listcore(ncore) + double precision, intent(in) :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg + double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n) + + integer :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen + integer :: ii, jj, i_good, j_good, n_real + integer :: icount_eigval + logical, allocatable :: is_ok(:) + integer, allocatable :: list_degen(:,:), list_same_degen(:) + integer, allocatable :: iorder(:), list_degen_sorted(:) + double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:) + double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:) + + allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n)) + leigvec_unsrtd = 0.d0 + reigvec_unsrtd = 0.d0 + eigval_unsrtd = 0.d0 + + ! obtain degeneracies + allocate(list_degen(n,0:n)) + call give_degen_full_listcore(fock_diag, n, listcore, ncore, thr_deg, list_degen, n_degen_list) + + allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list)) + do i = 1, n_degen_list + n_degen = list_degen(i,0) + list_degen_sorted(i) = n_degen + iorder(i) = i + enddo + + ! sort by number of degeneracies + call isort(list_degen_sorted, iorder, n_degen_list) + + allocate(is_ok(n_degen_list)) + is_ok = .True. + icount_eigval = 0 + + ! loop over degeneracies + do i = 1, n_degen_list + if(.not.is_ok(i)) cycle + + is_ok(i) = .False. + n_degen = list_degen_sorted(i) + + + if(n_degen.ge.1000)then + print*,'core orbital ' + else + print *, ' diagonalizing for n_degen = ', n_degen + endif + + k = 1 + + ! group all the entries having the same degeneracies +!! do while (list_degen_sorted(i+k)==n_degen) + do m = i+1, n_degen_list + if(list_degen_sorted(m)==n_degen) then + is_ok(i+k) = .False. + k += 1 + endif + enddo + + print *, ' number of identical degeneracies = ', k + if(n_degen.ge.1000)then + n_degen = 1 + endif + size_mat = k*n_degen + print *, ' size_mat = ', size_mat + allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat)) + allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat)) + ! group all the elements sharing the same degeneracy + icount = 0 + do j = 1, k ! jth set of degeneracy + index_degen = iorder(i+j-1) + do m = 1, n_degen + icount += 1 + list_same_degen(icount) = list_degen(index_degen,m) + enddo + enddo + + print *, ' list of elements ' + do icount = 1, size_mat + print *, icount, list_same_degen(icount) + enddo + + ! you copy subset of matrix elements having all the same degeneracy in mat_tmp + do ii = 1, size_mat + i_good = list_same_degen(ii) + do jj = 1, size_mat + j_good = list_same_degen(jj) + mat_tmp(jj,ii) = mat_ref(j_good,i_good) + enddo + enddo + + call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd & + , leigvec_tmp, reigvec_tmp & + , n_real, eigval_tmp ) + + do ii = 1, size_mat + icount_eigval += 1 + eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues + do jj = 1, size_mat ! copy the eigenvectors + j_good = list_same_degen(jj) + leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii) + reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii) + enddo + enddo + + deallocate(mat_tmp, list_same_degen) + deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp) + enddo + + if(icount_eigval .ne. n) then + print *, ' pb !! (icount_eigval.ne.n)' + print *, ' icount_eigval,n', icount_eigval, n + stop + endif + + deallocate(iorder) + allocate(iorder(n)) + do i = 1, n + iorder(i) = i + enddo + call dsort(eigval_unsrtd, iorder, n) + + do i = 1, n + print*,'sorted eigenvalues ' + i_good = iorder(i) + eigval(i) = eigval_unsrtd(i) + print*,'i,eigval(i) = ',i,eigval(i) + do j = 1, n + leigvec(j,i) = leigvec_unsrtd(j,i_good) + reigvec(j,i) = reigvec_unsrtd(j,i_good) + enddo + enddo + + deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd) + deallocate(list_degen) + deallocate(iorder, list_degen_sorted) + deallocate(is_ok) + +end + +! --- + +subroutine give_degen_full_listcore(A, n, listcore, ncore, thr, list_degen, n_degen_list) + + BEGIN_DOC + ! you enter with an array A(n) and spits out all the elements degenerated up to thr + ! + ! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL + ! + ! list_degen(i,0) = number of degenerate entries + ! + ! list_degen(i,1) = index of the first degenerate entry + ! + ! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries + ! + ! if list_degen(i,0) == 1 it means that there is no degeneracy for that element + ! + ! if list_degen(i,0) >= 1000 it means that it is core orbitals + END_DOC + + implicit none + + double precision, intent(in) :: A(n) + double precision, intent(in) :: thr + integer, intent(in) :: n,ncore, listcore(ncore) + integer, intent(out) :: list_degen(n,0:n), n_degen_list + integer :: i, j, icount, icheck,k + logical, allocatable :: is_ok(:) + + + allocate(is_ok(n)) + n_degen_list = 0 + is_ok = .True. + ! you first exclude the "core" orbitals + do i = 1, ncore + j=listcore(i) + is_ok(j) = .False. + enddo + do i = 1, n + if(.not.is_ok(i)) cycle + n_degen_list +=1 + is_ok(i) = .False. + list_degen(n_degen_list,1) = i + icount = 1 + do j = i+1, n + if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then + is_ok(j) = .False. + icount += 1 + list_degen(n_degen_list,icount) = j + endif + enddo + + list_degen(n_degen_list,0) = icount + enddo + ! you set all the core orbitals as separate entities + icheck = 0 + do i = 1, n_degen_list + icheck += list_degen(i,0) + enddo + if(icheck.ne.(n-ncore))then + print *, ' pb ! :: icheck.ne.n-ncore' + print *, icheck, n-ncore + stop + endif + k=1000 + do i = 1, ncore + n_degen_list+= 1 + j=listcore(i) + list_degen(n_degen_list,1) = i + list_degen(n_degen_list,0) = k + k+=1 + enddo + + + +end + +! --- + From 8c4a7226cdb7b594528db269cdb39c1d556b8bc6 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 10 Jul 2023 11:32:06 +0200 Subject: [PATCH 05/21] minor changes --- src/casscf_cipsi/casscf.irp.f | 10 +++++----- src/casscf_cipsi/save_energy.irp.f | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f index a2f3c5a7..02954ebf 100644 --- a/src/casscf_cipsi/casscf.irp.f +++ b/src/casscf_cipsi/casscf.irp.f @@ -17,14 +17,14 @@ end subroutine run implicit none double precision :: energy_old, energy, pt2_max_before, ept2_before,delta_E - logical :: converged,state_following_casscf_save + logical :: converged,state_following_casscf_cipsi_save integer :: iteration converged = .False. energy = 0.d0 mo_label = "MCSCF" iteration = 1 - state_following_casscf_save = state_following_casscf + state_following_casscf_cipsi_save = state_following_casscf state_following_casscf = .True. touch state_following_casscf ept2_before = 0.d0 @@ -44,8 +44,8 @@ subroutine run call write_double(6,energy,'CAS-SCF energy = ') if(n_states == 1)then double precision :: E_PT2, PT2 - call ezfio_get_casscf_energy_pt2(E_PT2) - call ezfio_get_casscf_energy(PT2) + call ezfio_get_casscf_cipsi_energy_pt2(E_PT2) + call ezfio_get_casscf_cipsi_energy(PT2) PT2 -= E_PT2 call write_double(6,E_PT2,'E + PT2 energy = ') call write_double(6,PT2,' PT2 = ') @@ -98,7 +98,7 @@ subroutine run SOFT_TOUCH pt2_max endif if(iteration .gt. 3)then - state_following_casscf = state_following_casscf_save + state_following_casscf = state_following_casscf_cipsi_save soft_touch state_following_casscf endif endif diff --git a/src/casscf_cipsi/save_energy.irp.f b/src/casscf_cipsi/save_energy.irp.f index 8729c5af..18750c3d 100644 --- a/src/casscf_cipsi/save_energy.irp.f +++ b/src/casscf_cipsi/save_energy.irp.f @@ -4,6 +4,6 @@ subroutine save_energy(E,pt2) ! Saves the energy in |EZFIO|. END_DOC double precision, intent(in) :: E(N_states), pt2(N_states) - call ezfio_set_casscf_energy(E(1:N_states)) - call ezfio_set_casscf_energy_pt2(E(1:N_states)+pt2(1:N_states)) + call ezfio_set_casscf_cipsi_energy(E(1:N_states)) + call ezfio_set_casscf_cipsi_energy_pt2(E(1:N_states)+pt2(1:N_states)) end From 2d383d09c6381f4a421d39a45e27abe36fe99133 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 10 Jul 2023 12:10:53 +0200 Subject: [PATCH 06/21] routine htilde_mu_mat_opt_bi_ortho works --- src/bi_ort_ints/three_body_ijmk.irp.f | 7 ++ src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f | 93 +++++++++++++++++++++++++ src/tc_bi_ortho/slater_tc_3e_slow.irp.f | 2 +- src/tc_bi_ortho/test_tc_fock.irp.f | 9 +-- src/tc_scf/fock_tc.irp.f | 24 +++---- 5 files changed, 114 insertions(+), 21 deletions(-) diff --git a/src/bi_ort_ints/three_body_ijmk.irp.f b/src/bi_ort_ints/three_body_ijmk.irp.f index 0c02e4c5..669861b7 100644 --- a/src/bi_ort_ints/three_body_ijmk.irp.f +++ b/src/bi_ort_ints/three_body_ijmk.irp.f @@ -86,6 +86,13 @@ tmp_loc_1 = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) tmp_loc_2 = tmp_aux_2(ipoint,n) + tmp1(ipoint,1,n) = int2_grad1_u12_bimo_t(ipoint,1,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,1,k,i) * tmp_loc_2 + tmp1(ipoint,2,n) = int2_grad1_u12_bimo_t(ipoint,2,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,2,k,i) * tmp_loc_2 + tmp1(ipoint,3,n) = int2_grad1_u12_bimo_t(ipoint,3,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,3,k,i) * tmp_loc_2 + tmp1(ipoint,4,n) = int2_grad1_u12_bimo_t(ipoint,1,n,n) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & + + int2_grad1_u12_bimo_t(ipoint,2,n,n) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + + int2_grad1_u12_bimo_t(ipoint,3,n,n) * int2_grad1_u12_bimo_t(ipoint,3,k,i) + enddo enddo !$OMP END DO diff --git a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f index 1d1b26cc..8524253a 100644 --- a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f +++ b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f @@ -90,3 +90,96 @@ subroutine htcdag_bi_ortho_calc_tdav_slow(v, u, N_st, sze) end +subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_chi_array,i_H_phi_array) + use bitmasks + implicit none + BEGIN_DOC +! Computes $\langle i|H|Phi \rangle = \sum_J c^R_J \langle i | H | J \rangle$. +! +! AND $\langle Chi|H| i \rangle = \sum_J c^L_J \langle J | H | i \rangle$. +! +! CONVENTION: i_H_phi_array(0) = total matrix element, +! +! i_H_phi_array(1) = one-electron matrix element, +! +! i_H_phi_array(2) = two-electron matrix element, +! +! i_H_phi_array(3) = three-electron matrix element, +! +! Uses filter_connected_i_H_psi0 to get all the $|J \rangle$ to which $|i \rangle$ +! is connected. +! +! The i_H_psi_minilist is much faster but requires to build the +! minilists. + END_DOC + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef_l(Ndet_max,Nstate),coef_r(Ndet_max,Nstate) + double precision, intent(out) :: i_H_chi_array(0:3,Nstate),i_H_phi_array(0:3,Nstate) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hmono, htwoe, hthree, htot + integer, allocatable :: idx(:) + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + allocate(idx(0:Ndet)) + + i_H_chi_array = 0.d0 + i_H_phi_array = 0.d0 + + call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) + if (Nstate == 1) then + + do ii=1,idx(0) + i = idx(ii) + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot) + i_H_chi_array(0,1) = i_H_chi_array(0,1) + coef_l(i,1)*htot + i_H_chi_array(1,1) = i_H_chi_array(1,1) + coef_l(i,1)*hmono + i_H_chi_array(2,1) = i_H_chi_array(2,1) + coef_l(i,1)*htwoe + i_H_chi_array(3,1) = i_H_chi_array(3,1) + coef_l(i,1)*hthree + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot) + i_H_phi_array(0,1) = i_H_phi_array(0,1) + coef_r(i,1)*htot + i_H_phi_array(1,1) = i_H_phi_array(1,1) + coef_r(i,1)*hmono + i_H_phi_array(2,1) = i_H_phi_array(2,1) + coef_r(i,1)*htwoe + i_H_phi_array(3,1) = i_H_phi_array(3,1) + coef_r(i,1)*hthree + enddo + + else + + do ii=1,idx(0) + i = idx(ii) + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot) + do j = 1, Nstate + i_H_chi_array(0,j) = i_H_chi_array(0,j) + coef_l(i,j)*htot + i_H_chi_array(1,j) = i_H_chi_array(1,j) + coef_l(i,j)*hmono + i_H_chi_array(2,j) = i_H_chi_array(2,j) + coef_l(i,j)*htwoe + i_H_chi_array(3,j) = i_H_chi_array(3,j) + coef_l(i,j)*hthree + enddo + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot) + do j = 1, Nstate + i_H_phi_array(0,j) = i_H_phi_array(0,j) + coef_r(i,j)*htot + i_H_phi_array(1,j) = i_H_phi_array(1,j) + coef_r(i,j)*hmono + i_H_phi_array(2,j) = i_H_phi_array(2,j) + coef_r(i,j)*htwoe + i_H_phi_array(3,j) = i_H_phi_array(3,j) + coef_r(i,j)*hthree + enddo + enddo + + endif + +end + diff --git a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f index 35abbbc4..cb33d343 100644 --- a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f @@ -184,7 +184,7 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree) ii = occ(i,s1) do j = i+1, Ne(s1) jj = occ(j,s1) -! ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1) +! !ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1) hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR enddo enddo diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f index 182c03d7..f1a7cc0a 100644 --- a/src/tc_bi_ortho/test_tc_fock.irp.f +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -152,9 +152,7 @@ subroutine routine_tot() ! do i = 1, elec_num_tab(s1) ! do a = elec_num_tab(s1)+1, mo_num ! virtual do i = 1, elec_beta_num - do a = elec_beta_num+1, elec_alpha_num! virtual -! do i = elec_beta_num+1, elec_alpha_num -! do a = elec_alpha_num+1, mo_num! virtual + do a = elec_beta_num+1, mo_num! virtual print*,i,a det_i = ref_bitmask @@ -167,7 +165,7 @@ subroutine routine_tot() call htilde_mu_mat_bi_ortho_slow(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij) print*,htilde_ij - if(dabs(htilde_ij).lt.1.d-10)cycle +! if(dabs(htilde_ij).lt.1.d-10)cycle print*, ' excited det' call debug_det(det_i, N_int) @@ -184,9 +182,12 @@ subroutine routine_tot() ! endif err_ai = dabs(dabs(ref) - dabs(new)) if(err_ai .gt. 1d-7) then + print*,'---------' print*,'s1 = ',s1 print*, ' warning on', i, a print*, ref,new,err_ai + print*,hmono, htwoe, hthree + print*,'---------' endif print*, ref,new,err_ai err_tot += err_ai diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 0ae515bb..f4553f3e 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -208,10 +208,10 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] if(three_body_h_tc) then !call wall_time(tt0) - !PROVIDE fock_a_tot_3e_bi_orth - !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth - PROVIDE fock_3e_uhf_mo_a - Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a + PROVIDE fock_a_tot_3e_bi_orth + Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth +! PROVIDE fock_3e_uhf_mo_a +! Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a !call wall_time(tt1) !print*, ' 3-e term:', tt1-tt0 endif @@ -241,21 +241,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] if(bi_ortho) then - !allocate(tmp(ao_num,ao_num)) - !tmp = Fock_matrix_tc_ao_beta - !if(three_body_h_tc) then - ! tmp += fock_3e_uhf_ao_b - !endif - !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1)) - !deallocate(tmp) - call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) if(three_body_h_tc) then - !PROVIDE fock_b_tot_3e_bi_orth - !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth - PROVIDE fock_3e_uhf_mo_b - Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b + PROVIDE fock_b_tot_3e_bi_orth + Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth +! PROVIDE fock_3e_uhf_mo_b +! Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b endif else From 8729a7ca451bde897be57d61358643bb3cd84229 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 10 Jul 2023 23:24:12 +0200 Subject: [PATCH 07/21] inactive --> virtual one-e term gradient ok --- src/casscf_tc_bi/grad_dm.irp.f | 78 +++++++++++++++++++++++ src/casscf_tc_bi/grad_old.irp.f | 109 ++++++++++++++++++++++++++++++++ src/casscf_tc_bi/gradient.irp.f | 94 +++++++++++++++++++++++++++ 3 files changed, 281 insertions(+) create mode 100644 src/casscf_tc_bi/grad_dm.irp.f create mode 100644 src/casscf_tc_bi/grad_old.irp.f create mode 100644 src/casscf_tc_bi/gradient.irp.f diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f new file mode 100644 index 00000000..0fc2e4eb --- /dev/null +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -0,0 +1,78 @@ + BEGIN_PROVIDER [real*8, gradvec_tc_r, (0:3,nMonoEx)] +&BEGIN_PROVIDER [real*8, gradvec_tc_l, (0:3,nMonoEx)] + implicit none + integer :: ii,tt,aa,indx + integer :: i,t,a,fff + double precision :: res_l(0:3), res_r(0:3) + gradvec_tc_l = 0.d0 + gradvec_tc_r = 0.d0 + do i=1,n_core_inact_orb + ii=list_core_inact(i) + do t=1,n_act_orb + tt=list_act(t) + indx = mat_idx_c_a(i,t) + call gradvec_tc_it(ii,tt,res_l) + call gradvec_tc_it(tt,ii,res_r) + do fff = 0,3 + gradvec_tc_l(fff,indx)=res_l(fff) + gradvec_tc_r(fff,indx)=res_r(fff) + enddo + end do + end do + + do i=1,n_core_inact_orb + ii=list_core_inact(i) + do a=1,n_virt_orb + indx = mat_idx_c_v(i,a) + aa=list_virt(a) + call gradvec_tc_ia(ii,aa,res_l) + call gradvec_tc_ia(aa,ii,res_r) + do fff = 0,3 + gradvec_tc_l(fff,indx)=res_l(fff) + gradvec_tc_r(fff,indx)=res_r(fff) + enddo + end do + end do + + do t=1,n_act_orb + do a=1,n_virt_orb + indx = mat_idx_a_v(i,a) +! gradvec_tc_l(indx)=gradvec_ta(t,a) + end do + end do +END_PROVIDER + +subroutine gradvec_tc_ia(i,a,res) + implicit none + BEGIN_DOC +! doubly occupied --> virtual TC gradient +! +! Corresponds to + END_DOC + integer, intent(in) :: i,a + double precision, intent(out) :: res(0:3) + res = 0.d0 + res(1) = -2 * mo_bi_ortho_tc_one_e(i,a) + +end + +subroutine gradvec_tc_it(i,t,res) + implicit none + BEGIN_DOC +! doubly occupied --> active TC gradient +! +! Corresponds to + END_DOC + integer, intent(in) :: i,t + double precision, intent(out) :: res(0:3) + integer :: rr,r,ss,s + double precision :: dm + res = 0.d0 + res(1) = -2 * mo_bi_ortho_tc_one_e(i,t) + do rr = 1, n_act_orb + r = list_act(rr) + dm = tc_transition_matrix_mo(t,r,1,1) + res(1) += mo_bi_ortho_tc_one_e(i,r) * dm + enddo + +end diff --git a/src/casscf_tc_bi/grad_old.irp.f b/src/casscf_tc_bi/grad_old.irp.f new file mode 100644 index 00000000..6610dee3 --- /dev/null +++ b/src/casscf_tc_bi/grad_old.irp.f @@ -0,0 +1,109 @@ + + BEGIN_PROVIDER [real*8, gradvec_detail_right_old, (0:3,nMonoEx)] +&BEGIN_PROVIDER [real*8, gradvec_detail_left_old, (0:3,nMonoEx)] + BEGIN_DOC + ! calculate the orbital gradient by hand, i.e. for + ! each determinant I we determine the string E_pq |I> (alpha and beta + ! separately) and generate + ! sum_I c_I is then the pq component of the orbital + ! gradient + ! E_pq = a^+_pa_q + a^+_Pa_Q + END_DOC + implicit none + integer :: ii,tt,aa,indx,ihole,ipart,istate,ll + real*8 :: res_l(0:3), res_r(0:3) + + do ii = 1, n_core_inact_orb + ihole = list_core_inact(ii) + do aa = 1, n_virt_orb + ipart = list_virt(aa) + indx = mat_idx_c_v(ii,aa) + call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) + do ll = 0, 3 + gradvec_detail_left_old (ll,indx)=res_l(ll) + gradvec_detail_right_old(ll,indx)=res_r(ll) + enddo + enddo + enddo +! do indx=1,nMonoEx +! ihole=excit(1,indx) +! ipart=excit(2,indx) +! call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) +! do ll = 0, 3 +! gradvec_detail_left_old (ll,indx)=res_l(ll) +! gradvec_detail_right_old(ll,indx)=res_r(ll) +! enddo +! end do + + real*8 :: norm_grad_left, norm_grad_right + norm_grad_left=0.d0 + norm_grad_right=0.d0 + do indx=1,nMonoEx + norm_grad_left+=gradvec_detail_left_old(0,indx)*gradvec_detail_left_old(0,indx) + norm_grad_right+=gradvec_detail_right_old(0,indx)*gradvec_detail_right_old(0,indx) + end do + norm_grad_left=sqrt(norm_grad_left) + norm_grad_right=sqrt(norm_grad_right) +! if (bavard) then + write(6,*) + write(6,*) ' Norm of the LEFT orbital gradient (via <0|EH|0>) : ', norm_grad_left + write(6,*) ' Norm of the RIGHT orbital gradient (via <0|HE|0>) : ', norm_grad_right + write(6,*) +! endif + + +END_PROVIDER + +subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) + BEGIN_DOC + ! eq 18 of Siegbahn et al, Physica Scripta 1980 + ! we calculate res_l = , and res_r = + ! q=hole, p=particle + ! res_l(0) = total matrix element + ! res_l(1) = one-electron part + ! res_l(2) = two-electron part + ! res_l(3) = three-electron part + END_DOC + implicit none + integer, intent(in) :: ihole,ipart + double precision, intent(out) :: res_l(0:3), res_r(0:3) + integer :: mu,iii,ispin,ierr,nu,istate,ll + integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) + real*8 :: i_H_chi_array(0:3,N_states),i_H_phi_array(0:3,N_states),phase + allocate(det_mu(N_int,2)) + allocate(det_mu_ex(N_int,2)) + + res_l=0.D0 + res_r=0.D0 + +! print*,'in i_h_psi' +! print*,ihole,ipart + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation on it + call det_copy(det_mu,det_mu_ex,N_int) + call do_signed_mono_excitation(det_mu,det_mu_ex,nu & + ,ihole,ipart,ispin,phase,ierr) + if (ierr.eq.1) then + + call i_H_tc_psi_phi(det_mu_ex,psi_det,psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,N_int & + ,N_det,N_det,N_states,i_H_chi_array,i_H_phi_array) +! print*,i_H_chi_array(1,1),i_H_phi_array(1,1) + do istate=1,N_states + do ll = 0,3 + res_l(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase + res_r(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase + enddo + end do + end if + end do + end do + + ! state-averaged gradient + res_l*=1.d0/dble(N_states) + res_r*=1.d0/dble(N_states) + +end + diff --git a/src/casscf_tc_bi/gradient.irp.f b/src/casscf_tc_bi/gradient.irp.f new file mode 100644 index 00000000..630bd891 --- /dev/null +++ b/src/casscf_tc_bi/gradient.irp.f @@ -0,0 +1,94 @@ +use bitmasks + +BEGIN_PROVIDER [ integer, nMonoEx ] + BEGIN_DOC + ! Number of single excitations + END_DOC + implicit none + 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, n_c_a_prov] +&BEGIN_PROVIDER [integer, n_c_v_prov] +&BEGIN_PROVIDER [integer, n_a_v_prov] + implicit none + n_c_a_prov = n_core_inact_orb * n_act_orb + n_c_v_prov = n_core_inact_orb * n_virt_orb + n_a_v_prov = n_act_orb * n_virt_orb + END_PROVIDER + + BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] +&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] +&BEGIN_PROVIDER [integer, list_idx_c_a, (3,n_c_a_prov) ] +&BEGIN_PROVIDER [integer, list_idx_c_v, (3,n_c_v_prov) ] +&BEGIN_PROVIDER [integer, list_idx_a_v, (3,n_a_v_prov) ] +&BEGIN_PROVIDER [integer, mat_idx_c_a, (n_core_inact_orb,n_act_orb) +&BEGIN_PROVIDER [integer, mat_idx_c_v, (n_core_inact_orb,n_virt_orb) +&BEGIN_PROVIDER [integer, mat_idx_a_v, (n_act_orb,n_virt_orb) + BEGIN_DOC + ! a list of the orbitals involved in the excitation + END_DOC + + implicit none + integer :: i,t,a,ii,tt,aa,indx,indx_tmp + indx=0 + indx_tmp = 0 + do ii=1,n_core_inact_orb + i=list_core_inact(ii) + do tt=1,n_act_orb + t=list_act(tt) + indx+=1 + excit(1,indx)=i + excit(2,indx)=t + excit_class(indx)='c-a' + indx_tmp += 1 + list_idx_c_a(1,indx_tmp) = indx + list_idx_c_a(2,indx_tmp) = ii + list_idx_c_a(3,indx_tmp) = tt + mat_idx_c_a(ii,tt) = indx + end do + end do + + indx_tmp = 0 + do ii=1,n_core_inact_orb + i=list_core_inact(ii) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=i + excit(2,indx)=a + excit_class(indx)='c-v' + indx_tmp += 1 + list_idx_c_v(1,indx_tmp) = indx + list_idx_c_v(2,indx_tmp) = ii + list_idx_c_v(3,indx_tmp) = aa + mat_idx_c_v(ii,aa) = indx + end do + end do + + indx_tmp = 0 + do tt=1,n_act_orb + t=list_act(tt) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=t + excit(2,indx)=a + excit_class(indx)='a-v' + indx_tmp += 1 + list_idx_a_v(1,indx_tmp) = indx + list_idx_a_v(2,indx_tmp) = tt + list_idx_a_v(3,indx_tmp) = aa + mat_idx_a_v(tt,aa) = indx + end do + end do + +! if (bavard) then + write(6,*) ' Filled the table of the Monoexcitations ' + do indx=1,nMonoEx + write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' & + ,excit(2,indx),' ',excit_class(indx) + end do +! end if + +END_PROVIDER From b3b080929b38c152ac1bd868bad4dd311108e7e9 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 11 Jul 2023 01:56:28 +0200 Subject: [PATCH 08/21] fixed stupid bug in TC 1-RDM and one-e gradient: o-v, o-a are ok --- src/casscf_tc_bi/grad_dm.irp.f | 28 ++++++++++++++++------------ src/casscf_tc_bi/grad_old.irp.f | 21 +++++++++++++++++---- src/tc_bi_ortho/tc_prop.irp.f | 8 +++++--- 3 files changed, 38 insertions(+), 19 deletions(-) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index 0fc2e4eb..7f6155ab 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -11,8 +11,7 @@ do t=1,n_act_orb tt=list_act(t) indx = mat_idx_c_a(i,t) - call gradvec_tc_it(ii,tt,res_l) - call gradvec_tc_it(tt,ii,res_r) + call gradvec_tc_it(ii,tt,res_l,res_r) do fff = 0,3 gradvec_tc_l(fff,indx)=res_l(fff) gradvec_tc_r(fff,indx)=res_r(fff) @@ -56,23 +55,28 @@ subroutine gradvec_tc_ia(i,a,res) end -subroutine gradvec_tc_it(i,t,res) +subroutine gradvec_tc_it(i,t,res_l, res_r) implicit none BEGIN_DOC ! doubly occupied --> active TC gradient ! -! Corresponds to +! Corresponds to res_r = +! +! res_l = END_DOC integer, intent(in) :: i,t - double precision, intent(out) :: res(0:3) - integer :: rr,r,ss,s + double precision, intent(out) :: res_l(0:3),res_r(0:3) + integer :: rr,r,ss,s,m double precision :: dm - res = 0.d0 - res(1) = -2 * mo_bi_ortho_tc_one_e(i,t) - do rr = 1, n_act_orb - r = list_act(rr) - dm = tc_transition_matrix_mo(t,r,1,1) - res(1) += mo_bi_ortho_tc_one_e(i,r) * dm + res_r = 0.d0 + do m = 1, mo_num + res_r(1) += mo_bi_ortho_tc_one_e(i,m) * tc_transition_matrix_mo(t,m,1,1) & + -mo_bi_ortho_tc_one_e(m,t) * tc_transition_matrix_mo(m,i,1,1) + enddo + res_l = 0.d0 + do m = 1, mo_num + res_l(1) += mo_bi_ortho_tc_one_e(t,m) * tc_transition_matrix_mo(i,m,1,1) & + -mo_bi_ortho_tc_one_e(m,i) * tc_transition_matrix_mo(m,t,1,1) enddo end diff --git a/src/casscf_tc_bi/grad_old.irp.f b/src/casscf_tc_bi/grad_old.irp.f index 6610dee3..ea6747b1 100644 --- a/src/casscf_tc_bi/grad_old.irp.f +++ b/src/casscf_tc_bi/grad_old.irp.f @@ -25,6 +25,19 @@ enddo enddo enddo + + do ii = 1, n_core_inact_orb + ihole = list_core_inact(ii) + do tt = 1, n_act_orb + ipart = list_act(tt) + indx = mat_idx_c_a(ii,tt) + call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) + do ll = 0, 3 + gradvec_detail_left_old (ll,indx)=res_l(ll) + gradvec_detail_right_old(ll,indx)=res_r(ll) + enddo + enddo + enddo ! do indx=1,nMonoEx ! ihole=excit(1,indx) ! ipart=excit(2,indx) @@ -57,7 +70,7 @@ END_PROVIDER subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) BEGIN_DOC ! eq 18 of Siegbahn et al, Physica Scripta 1980 - ! we calculate res_l = , and res_r = + ! we calculate res_r = , and res_r = ! q=hole, p=particle ! res_l(0) = total matrix element ! res_l(1) = one-electron part @@ -89,12 +102,12 @@ subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) if (ierr.eq.1) then call i_H_tc_psi_phi(det_mu_ex,psi_det,psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,N_int & - ,N_det,N_det,N_states,i_H_chi_array,i_H_phi_array) + ,N_det,psi_det_size,N_states,i_H_chi_array,i_H_phi_array) ! print*,i_H_chi_array(1,1),i_H_phi_array(1,1) do istate=1,N_states do ll = 0,3 - res_l(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase - res_r(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase + res_l(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase + res_r(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase enddo end do end if diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f index 5bb0e2c0..a13dc9a2 100644 --- a/src/tc_bi_ortho/tc_prop.irp.f +++ b/src/tc_bi_ortho/tc_prop.irp.f @@ -29,7 +29,7 @@ tc_transition_matrix_mo_alpha(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) enddo do p = 1, n_occ_ab(2) ! browsing the beta electrons - m = occ(p,1) + m = occ(p,2) tc_transition_matrix_mo_beta(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) enddo else @@ -38,12 +38,14 @@ ! Single alpha h = exc(1,1,1) ! hole in psi_det(1,1,j) p = exc(1,2,1) ! particle in psi_det(1,1,j) - tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= & + phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) else ! Single beta h = exc(1,1,2) ! hole in psi_det(1,1,j) p = exc(1,2,2) ! particle in psi_det(1,1,j) - tc_transition_matrix_mo_beta(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + tc_transition_matrix_mo_beta(p,h,istate,jstate)+= & + phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) endif endif enddo From 467f7563797c4c32fc9597542b9d761309e3565c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 16 Jul 2023 20:04:17 +0200 Subject: [PATCH 09/21] Optimized A1 in CCSD --- src/ccsd/ccsd_space_orb_sub_chol.irp.f | 48 ++++---------------------- 1 file changed, 6 insertions(+), 42 deletions(-) diff --git a/src/ccsd/ccsd_space_orb_sub_chol.irp.f b/src/ccsd/ccsd_space_orb_sub_chol.irp.f index 0ba46e56..ec6c2afb 100644 --- a/src/ccsd/ccsd_space_orb_sub_chol.irp.f +++ b/src/ccsd/ccsd_space_orb_sub_chol.irp.f @@ -1023,56 +1023,26 @@ subroutine compute_A1_chol(nO,nV,t1,t2,tau,A1) integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta - double precision, allocatable :: X_vooo(:,:,:,:), Y_oooo(:,:,:,:) - allocate(X_vooo(nV,nO,nO,nO), Y_oooo(nO,nO,nO,nO)) + double precision, allocatable :: Y_oooo(:,:,:,:) + allocate(Y_oooo(nO,nO,nO,nO)) ! A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j) - !$omp parallel & - !$omp shared(nO,nV,A1,cc_space_v_oooo,cc_space_v_ovoo,X_vooo) & - !$omp private(u,v,i,j) & - !$omp default(none) - !$omp do collapse(2) - do j = 1, nO - do i = 1, nO - do v = 1, nO - do u = 1, nO - A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j) - enddo - enddo - enddo - enddo - !$omp end do nowait - ! A1(u,v,i,j) += cc_space_v_ovoo(u,a,i,j) * t1(v,a) & - !$omp do collapse(2) - do j = 1, nO - do i = 1, nO - do u = 1, nO - do a = 1, nV - X_vooo(a,u,i,j) = cc_space_v_ovoo(u,a,i,j) - enddo - enddo - enddo - enddo - !$omp end do - !$omp end parallel - call dgemm('N','N', nO, nO*nO*nO, nV, & 1d0, t1 , size(t1,1), & - X_vooo, size(X_vooo,1), & + cc_space_v_vooo, size(cc_space_v_vooo,1), & 0d0, Y_oooo, size(Y_oooo,1)) !$omp parallel & - !$omp shared(nO,nV,A1,Y_oooo) & !$omp private(u,v,i,j) & - !$omp default(none) + !$omp default(shared) !$omp do collapse(2) do j = 1, nO do i = 1, nO do v = 1, nO do u = 1, nO - A1(u,v,i,j) = A1(u,v,i,j) + Y_oooo(v,u,i,j) + A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j) + Y_oooo(v,u,j,i) + Y_oooo(u,v,i,j) enddo enddo enddo @@ -1080,13 +1050,7 @@ subroutine compute_A1_chol(nO,nV,t1,t2,tau,A1) !$omp end do !$omp end parallel - deallocate(X_vooo,Y_oooo) - - ! A1(u,v,i,j) += cc_space_v_vooo(a,v,i,j) * t1(u,a) - call dgemm('N','N', nO, nO*nO*nO, nV, & - 1d0, t1 , size(t1,1), & - cc_space_v_vooo, size(cc_space_v_vooo,1), & - 1d0, A1 , size(A1,1)) + deallocate(Y_oooo) ! A1(u,v,i,j) += cc_space_v_vvoo(a,b,i,j) * tau(u,v,a,b) call dgemm('N','N', nO*nO, nO*nO, nV*nV, & From bd570b19c1d4ae0479b8ff4c4611ae1127605441 Mon Sep 17 00:00:00 2001 From: ydamour Date: Mon, 17 Jul 2023 17:05:48 +0200 Subject: [PATCH 10/21] fix bug restore_symmetry --- src/utils/linear_algebra.irp.f | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 65c57a76..314ad4f6 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1661,7 +1661,15 @@ subroutine restore_symmetry(m,n,A,LDA,thresh) ! Update i i = i + 1 enddo - copy(i:) = 0.d0 + + ! To nullify the remaining elements that are below the threshold + if (i == sze) then + if (-copy(i) <= thresh) then + copy(i) = 0d0 + endif + else + copy(i:) = 0.d0 + endif !$OMP PARALLEL if (sze>10000) & !$OMP SHARED(m,sze,copy_sign,copy,key,A,ii,jj) & From cc7b97c09b5f8a970319a3e247551c34401e731c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 5 Aug 2023 01:47:48 +0200 Subject: [PATCH 11/21] Cleaning in C --- external/ezfio | 2 +- external/irpf90 | 2 +- src/utils/fortran_mmap.c | 1 - 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/external/ezfio b/external/ezfio index ed1df9f3..d5805497 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c +Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93 diff --git a/external/irpf90 b/external/irpf90 index 33ca5e10..0007f72f 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271 +Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102 diff --git a/src/utils/fortran_mmap.c b/src/utils/fortran_mmap.c index 71426002..e8d85a2f 100644 --- a/src/utils/fortran_mmap.c +++ b/src/utils/fortran_mmap.c @@ -9,7 +9,6 @@ void* mmap_fortran(char* filename, size_t bytes, int* file_descr, int read_only) { - int i; int fd; int result; void* map; From c6b50d5f500faed66fdc13a2b2c4d2dd874fbca1 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 7 Aug 2023 16:39:53 +0200 Subject: [PATCH 12/21] found a bug in left-right coefficients in perturbation --- src/cipsi_tc_bi_ortho/selection.irp.f | 35 +++++++++++++++++++++++---- src/fci_tc_bi/selectors.irp.f | 2 ++ src/tc_bi_ortho/tc_h_eigvectors.irp.f | 5 ++++ 3 files changed, 37 insertions(+), 5 deletions(-) diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 77377554..06cf848b 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -893,20 +893,45 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function enddo else if(debug_tc_pt2 == 2)then !! debugging the new version +! psi_h_alpha_tmp = 0.d0 +! alpha_h_psi_tmp = 0.d0 +! do iii = 1, N_det_selectors ! old version +! call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) +! call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) +! psi_h_alpha_tmp += i_h_alpha * psi_selectors_coef_tc(iii,1,1) ! left function +! alpha_h_psi_tmp += alpha_h_i * psi_selectors_coef_tc(iii,2,1) ! right function +! enddo psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version psi_h_alpha = 0.d0 alpha_h_psi = 0.d0 - do iii = 1, N_det_selectors ! old version - call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) - call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) - psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function - alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function + do iii = 1, N_det ! old version + call htilde_mu_mat_opt_bi_ortho_no_3e(psi_det(1,1,iii), det, N_int, i_h_alpha) + call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_det(1,1,iii), N_int, alpha_h_i) + psi_h_alpha += i_h_alpha * psi_l_coef_bi_ortho(iii,1) ! left function + alpha_h_psi += alpha_h_i * psi_r_coef_bi_ortho(iii,1) ! right function enddo if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi) if(error.gt.1.d-2)then + call debug_det(det, N_int) print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E + print*,psi_h_alpha , alpha_h_psi + print*,psi_h_alpha_tmp , alpha_h_psi_tmp + print*,'selectors ' + do iii = 1, N_det_selectors ! old version + print*,'iii',iii,psi_selectors_coef_tc(iii,1,1),psi_selectors_coef_tc(iii,2,1) + call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) + call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) + print*,i_h_alpha,alpha_h_i + call debug_det(psi_selectors(1,1,iii),N_int) + enddo +! print*,'psi_det ' +! do iii = 1, N_det! old version +! print*,'iii',iii,psi_l_coef_bi_ortho(iii,1),psi_r_coef_bi_ortho(iii,1) +! call debug_det(psi_det(1,1,iii),N_int) +! enddo + stop endif endif else diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f index 4d3de7d0..7f93ae55 100644 --- a/src/fci_tc_bi/selectors.irp.f +++ b/src/fci_tc_bi/selectors.irp.f @@ -27,6 +27,8 @@ END_PROVIDER implicit none BEGIN_DOC ! Determinants on which we apply for perturbation. + ! psi_selectors_coef_tc(iii,1,istate) = left coefficient of the iii determinant + ! psi_selectors_coef_tc(iii,2,istate) = right coefficient of the iii determinant END_DOC integer :: i,k diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f index f027c38f..48257943 100644 --- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -328,6 +328,11 @@ end TOUCH psi_r_coef_bi_ortho call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(buffer) deallocate(buffer) +! print*,'After diag' +! do i = 1, N_det! old version +! print*,'i',i,psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1) +! call debug_det(psi_det(1,1,i),N_int) +! enddo END_PROVIDER From c945e027c03c283682f6bffdb933c588069450bf Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 7 Aug 2023 16:56:10 +0200 Subject: [PATCH 13/21] fixed a bug in psi_selectors_coef --- src/cipsi_tc_bi_ortho/get_d0_good.irp.f | 8 +++--- src/cipsi_tc_bi_ortho/get_d1_good.irp.f | 36 ++++++++++++------------- src/cipsi_tc_bi_ortho/get_d2_good.irp.f | 28 +++++++++---------- 3 files changed, 36 insertions(+), 36 deletions(-) diff --git a/src/cipsi_tc_bi_ortho/get_d0_good.irp.f b/src/cipsi_tc_bi_ortho/get_d0_good.irp.f index 4270e7b8..9bba162e 100644 --- a/src/cipsi_tc_bi_ortho/get_d0_good.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d0_good.irp.f @@ -53,7 +53,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (hij == (0.d0,0.d0)) cycle !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * hij ! HOTSPOT + mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij ! HOTSPOT enddo end do !!!!!!!!!! @@ -72,7 +72,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (hji == (0.d0,0.d0)) cycle !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,2) * hji ! HOTSPOT + mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji ! HOTSPOT enddo end do end do @@ -109,7 +109,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, end if !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij enddo end do @@ -128,7 +128,7 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, end if !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji enddo end do end do diff --git a/src/cipsi_tc_bi_ortho/get_d1_good.irp.f b/src/cipsi_tc_bi_ortho/get_d1_good.irp.f index bc19e7e4..b2a38e02 100644 --- a/src/cipsi_tc_bi_ortho/get_d1_good.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d1_good.irp.f @@ -76,7 +76,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,1) + tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2) enddo endif end do @@ -88,7 +88,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,1) + tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2) enddo endif end do @@ -114,7 +114,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,2) + tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1) enddo endif end do @@ -126,7 +126,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,2) + tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1) enddo endif end do @@ -169,7 +169,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,1) + tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2) enddo endif end if @@ -180,7 +180,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) do k=1,N_states - tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,1) + tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2) enddo endif end if @@ -211,7 +211,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,2) + tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1) enddo endif end if @@ -222,7 +222,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (hji /= 0.d0) then hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) do k=1,N_states - tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,2) + tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1) enddo endif end if @@ -265,7 +265,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij_cache(putj,1) - hij_cache(putj,2) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,1) + tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2) endif end do do putj=hfix+1,mo_num @@ -274,7 +274,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij_cache(putj,2) - hij_cache(putj,1) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,1) + tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2) endif end do @@ -293,7 +293,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji_cache(putj,1) - hji_cache(putj,2) if (hji /= 0.d0) then hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,2) + tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1) endif end do do putj=hfix+1,mo_num @@ -302,7 +302,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji_cache(putj,2) - hji_cache(putj,1) if (hji /= 0.d0) then hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,2) + tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1) endif end do @@ -342,7 +342,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,1) + tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2) enddo endif end if @@ -353,7 +353,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) do k=1,N_states - tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,1) + tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2) enddo endif end if @@ -385,7 +385,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,2) + tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1) enddo endif end if @@ -396,7 +396,7 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (hji /= 0.d0) then hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) do k=1,N_states - tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,2) + tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1) enddo endif end if @@ -445,8 +445,8 @@ subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, do k=1,N_states ! take conjugate to get contribution to instead of ! mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * dconjg(hij) - mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * hij - mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,2) * hji + mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij + mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji enddo end do end do diff --git a/src/cipsi_tc_bi_ortho/get_d2_good.irp.f b/src/cipsi_tc_bi_ortho/get_d2_good.irp.f index 0a08c808..d01ed433 100644 --- a/src/cipsi_tc_bi_ortho/get_d2_good.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d2_good.irp.f @@ -79,12 +79,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,1) * hij + mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij enddo else ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij enddo end if end do @@ -103,12 +103,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,2) * hji + mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji enddo else ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji enddo end if end do @@ -135,7 +135,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij enddo endif end do @@ -154,7 +154,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji enddo endif end do @@ -189,7 +189,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,1) * hij + mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,2) * hij enddo end do end do @@ -210,7 +210,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,2) * hji + mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,1) * hji enddo end do end do @@ -239,12 +239,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (puti < putj) then !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij enddo else !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,1) * hij + mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij enddo endif end do @@ -262,12 +262,12 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, if (puti < putj) then !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji enddo else !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,2) * hji + mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji enddo endif end do @@ -290,7 +290,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hij = hij * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij enddo end if !! @@ -299,7 +299,7 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, hji = hji * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji enddo end if end if From 0440def3637b4fc528ffd488f117247662e12c4d Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 9 Aug 2023 16:23:09 +0200 Subject: [PATCH 14/21] added lccsd --- src/cisd/EZFIO.cfg | 8 ++++ src/cisd/NEED | 1 + src/cisd/lccsd.irp.f | 95 +++++++++++++++++++++++++++++++++++++++ src/cisd/lccsd_prov.irp.f | 44 ++++++++++++++++++ 4 files changed, 148 insertions(+) create mode 100644 src/cisd/lccsd.irp.f create mode 100644 src/cisd/lccsd_prov.irp.f diff --git a/src/cisd/EZFIO.cfg b/src/cisd/EZFIO.cfg index 4565d2df..688f802a 100644 --- a/src/cisd/EZFIO.cfg +++ b/src/cisd/EZFIO.cfg @@ -5,3 +5,11 @@ interface: ezfio size: (determinants.n_states) + +[lcc_energy] +type: double precision +doc: lccsd energy +interface: ezfio +size: (determinants.n_states) + + diff --git a/src/cisd/NEED b/src/cisd/NEED index d9ad3efc..616d021e 100644 --- a/src/cisd/NEED +++ b/src/cisd/NEED @@ -1,3 +1,4 @@ selectors_full single_ref_method davidson_undressed +dav_general_mat diff --git a/src/cisd/lccsd.irp.f b/src/cisd/lccsd.irp.f new file mode 100644 index 00000000..919c5aaa --- /dev/null +++ b/src/cisd/lccsd.irp.f @@ -0,0 +1,95 @@ +program lccsd + implicit none + BEGIN_DOC +! Linerarized CCSD +! + ! This program takes a reference Slater determinant of ROHF-like occupancy, + ! + ! and performs all single and double excitations on top of it, disregarding + ! spatial symmetry and compute the "n_states" lowest eigenstates of that CI + ! matrix (see :option:`determinants n_states`). + ! + ! This program can be useful in many cases: + ! + ! * **Ground state calculation**: if even after a :c:func:`cis` calculation, natural + ! orbitals (see :c:func:`save_natorb`) and then :c:func:`scf` optimization, you are not sure to have the lowest scf + ! solution, + ! do the same strategy with the :c:func:`cisd` executable instead of the :c:func:`cis` exectuable to generate the natural + ! orbitals as a guess for the :c:func:`scf`. + ! + ! + ! + ! * **Excited states calculations**: the lowest excited states are much likely to + ! be dominanted by single- or double-excitations. + ! Therefore, running a :c:func:`cisd` will save the "n_states" lowest states within + ! the CISD space + ! in the |EZFIO| directory, which can afterward be used as guess wave functions + ! for a further multi-state fci calculation if you specify "read_wf" = True + ! before running the fci executable (see :option:`determinants read_wf`). + ! Also, if you specify "s2_eig" = True, the cisd will only retain states + ! having the good value :math:`S^2` value + ! (see :option:`determinants expected_s2` and :option:`determinants s2_eig`). + ! If "s2_eig" = False, it will take the lowest n_states, whatever + ! multiplicity they are. + ! + ! + ! + ! Note: if you would like to discard some orbitals, use + ! :ref:`qp_set_mo_class` to specify: + ! + ! * "core" orbitals which will be always doubly occupied + ! + ! * "act" orbitals where an electron can be either excited from or to + ! + ! * "del" orbitals which will be never occupied + ! + END_DOC + PROVIDE N_states + read_wf = .False. + TOUCH read_wf + call run +end + +subroutine run + implicit none + + if(pseudo_sym)then + call H_apply_cisd_sym + else + call H_apply_cisd + endif + call get_lccsd_2 +end + +subroutine get_lccsd_2 + implicit none + integer :: i,k + double precision :: cisdq(N_states), delta_e + double precision,external :: diag_h_mat_elem + psi_coef = lccsd_coef + SOFT_TOUCH psi_coef + call save_wavefunction_truncated(save_threshold) + call ezfio_set_cisd_lcc_energy(lccsd_energies) + + print *, 'N_det = ', N_det + print*,'' + print*,'******************************' + print *, 'LCCSD Energies' + do i = 1,N_states + print *, i, lccsd_energies(i) + enddo + if (N_states > 1) then + print*,'******************************' + print*,'Excitation energies (au) (LCCSD)' + do i = 2, N_states + print*, i ,lccsd_energies(i) - lccsd_energies(1) + enddo + print*,'' + print*,'******************************' + print*,'Excitation energies (eV) (LCCSD)' + do i = 2, N_states + print*, i ,(lccsd_energies(i) - lccsd_energies(1)) * ha_to_ev + enddo + endif + +end diff --git a/src/cisd/lccsd_prov.irp.f b/src/cisd/lccsd_prov.irp.f new file mode 100644 index 00000000..b071a8f8 --- /dev/null +++ b/src/cisd/lccsd_prov.irp.f @@ -0,0 +1,44 @@ + BEGIN_PROVIDER [ double precision, lccsd_coef, (N_det, N_states)] +&BEGIN_PROVIDER [ double precision, lccsd_energies, (N_states)] + implicit none + double precision, allocatable :: Dress_jj(:), H_jj(:), u_in(:,:) + double precision :: ebefore, eafter, ecorr, thresh + integer :: i,it + logical :: converged + external H_u_0_nstates_openmp + allocate(Dress_jj(N_det),H_jj(N_det),u_in(N_det,N_states_diag)) + thresh = 1.d-6 + converged = .False. + Dress_jj = 0.d0 + u_in = 0.d0 + it = 0 + ! initial guess + do i = 1, N_states_diag + u_in(i,i) = 1.d0 + enddo + do i = 1,N_det + call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,H_jj(i)) + enddo + ebefore = H_jj(1) + do while (.not.converged) + it += 1 + print*,'N_det = ',N_det + call davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,lccsd_energies,& + N_det,N_states,N_states_diag,converged,H_u_0_nstates_openmp) + ecorr = lccsd_energies(1) - H_jj(1) + print*,'---------------------' + print*,'it = ',it + print*,'ecorr = ',ecorr + Dress_jj(1) = 0.d0 + do i = 2, N_det + Dress_jj(i) = ecorr + enddo + eafter = lccsd_energies(1) + converged = (dabs(eafter - ebefore).lt.thresh) + ebefore = eafter + enddo + do i = 1, N_states + lccsd_coef(1:N_det,i) = u_in(1:N_det,i) + enddo + +END_PROVIDER From 1a2632c280bc74eca95912fc7eec4e4b9d1f8587 Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 9 Aug 2023 16:59:47 +0200 Subject: [PATCH 15/21] added condition in lccsd --- src/cisd/lccsd_prov.irp.f | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/cisd/lccsd_prov.irp.f b/src/cisd/lccsd_prov.irp.f index b071a8f8..38149ac9 100644 --- a/src/cisd/lccsd_prov.irp.f +++ b/src/cisd/lccsd_prov.irp.f @@ -31,7 +31,9 @@ print*,'ecorr = ',ecorr Dress_jj(1) = 0.d0 do i = 2, N_det - Dress_jj(i) = ecorr + if(ecorr + H_jj(i) .gt. H_jj(1))then + Dress_jj(i) = ecorr + endif enddo eafter = lccsd_energies(1) converged = (dabs(eafter - ebefore).lt.thresh) From 2e458e93bafcd652e996cd3067d652f61522a915 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 10 Aug 2023 08:38:05 +0200 Subject: [PATCH 16/21] Fixing CI --- .github/workflows/configuration.yml | 2 +- external/ezfio | 2 +- external/irpf90 | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/configuration.yml b/.github/workflows/configuration.yml index ba37f5dd..178b394e 100644 --- a/.github/workflows/configuration.yml +++ b/.github/workflows/configuration.yml @@ -22,7 +22,7 @@ jobs: - uses: actions/checkout@v3 - name: Install dependencies run: | - sudo apt install gfortran gcc liblapack-dev libblas-dev wget python3 make m4 pkg-config hdf5 + sudo apt install gfortran gcc liblapack-dev libblas-dev wget python3 make m4 pkg-config libhdf5-dev - name: zlib run: | ./configure -i zlib || echo OK diff --git a/external/ezfio b/external/ezfio index d5805497..ed1df9f3 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93 +Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c diff --git a/external/irpf90 b/external/irpf90 index 0007f72f..33ca5e10 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102 +Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271 From ee2c470054b76bffc680f54b960daae99ed9679d Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 10 Aug 2023 15:53:35 +0200 Subject: [PATCH 17/21] clarified the TC-CASSCF gradients --- src/casscf_tc_bi/grad_dm.irp.f | 17 +++++----- src/casscf_tc_bi/grad_old.irp.f | 42 ++++++++++++++----------- src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f | 22 ++++++------- 3 files changed, 45 insertions(+), 36 deletions(-) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index 7f6155ab..1618adc6 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -24,8 +24,7 @@ do a=1,n_virt_orb indx = mat_idx_c_v(i,a) aa=list_virt(a) - call gradvec_tc_ia(ii,aa,res_l) - call gradvec_tc_ia(aa,ii,res_r) + call gradvec_tc_ia(ii,aa,res_l,res_r) do fff = 0,3 gradvec_tc_l(fff,indx)=res_l(fff) gradvec_tc_r(fff,indx)=res_r(fff) @@ -41,17 +40,21 @@ end do END_PROVIDER -subroutine gradvec_tc_ia(i,a,res) +subroutine gradvec_tc_ia(i,a,res_l, res_r) implicit none BEGIN_DOC ! doubly occupied --> virtual TC gradient ! -! Corresponds to +! Corresponds to res_r = , +! +! res_l = END_DOC integer, intent(in) :: i,a - double precision, intent(out) :: res(0:3) - res = 0.d0 - res(1) = -2 * mo_bi_ortho_tc_one_e(i,a) + double precision, intent(out) :: res_l(0:3), res_r(0:3) + res_l = 0.d0 + res_r = 0.d0 + res_l(1) = -2 * mo_bi_ortho_tc_one_e(a,i) + res_r(1) = -2 * mo_bi_ortho_tc_one_e(i,a) end diff --git a/src/casscf_tc_bi/grad_old.irp.f b/src/casscf_tc_bi/grad_old.irp.f index ea6747b1..6c976d66 100644 --- a/src/casscf_tc_bi/grad_old.irp.f +++ b/src/casscf_tc_bi/grad_old.irp.f @@ -69,45 +69,51 @@ END_PROVIDER subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) BEGIN_DOC - ! eq 18 of Siegbahn et al, Physica Scripta 1980 - ! we calculate res_r = , and res_r = - ! q=hole, p=particle - ! res_l(0) = total matrix element - ! res_l(1) = one-electron part - ! res_l(2) = two-electron part - ! res_l(3) = three-electron part + ! Computes the gradient with respect to orbital rotation BRUT FORCE + ! + ! res_l = + ! + ! res_r = + ! + ! q=hole, p=particle. NOTE that on res_l it is E_qp and on res_r it is E_pq + ! + ! res_l(0) = total matrix element, res_l(1) = one-electron part, + ! + ! res_l(2) = two-electron part, res_l(3) = three-electron part + ! END_DOC implicit none integer, intent(in) :: ihole,ipart double precision, intent(out) :: res_l(0:3), res_r(0:3) integer :: mu,iii,ispin,ierr,nu,istate,ll integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) - real*8 :: i_H_chi_array(0:3,N_states),i_H_phi_array(0:3,N_states),phase + real*8 :: chi_H_mu_ex_array(0:3,N_states),mu_ex_H_phi_array(0:3,N_states),phase allocate(det_mu(N_int,2)) allocate(det_mu_ex(N_int,2)) res_l=0.D0 res_r=0.D0 -! print*,'in i_h_psi' -! print*,ihole,ipart do mu=1,n_det - ! get the string of the determinant + ! get the string of the determinant |mu> call det_extract(det_mu,mu,N_int) do ispin=1,2 - ! do the monoexcitation on it + ! do the monoexcitation on it: |det_mu_ex> = a^dagger_{p,ispin} a_{q,ispin} |mu> call det_copy(det_mu,det_mu_ex,N_int) call do_signed_mono_excitation(det_mu,det_mu_ex,nu & ,ihole,ipart,ispin,phase,ierr) + ! |det_mu_ex> = a^dagger_{p,ispin} a_{q,ispin} |mu> if (ierr.eq.1) then - call i_H_tc_psi_phi(det_mu_ex,psi_det,psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,N_int & - ,N_det,psi_det_size,N_states,i_H_chi_array,i_H_phi_array) -! print*,i_H_chi_array(1,1),i_H_phi_array(1,1) + ,N_det,psi_det_size,N_states,chi_H_mu_ex_array,mu_ex_H_phi_array) + ! chi_H_mu_ex_array = + ! mu_ex_H_phi_array = do istate=1,N_states - do ll = 0,3 - res_l(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase - res_r(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase + do ll = 0,3 ! loop over the body components (1e,2e,3e) + !res_l = \sum_mu c_mu^l = + res_l(ll)+= mu_ex_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase + !res_r = \sum_mu c_mu^r = + res_r(ll)+= chi_H_mu_ex_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase enddo end do end if diff --git a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f index 8524253a..e96e738e 100644 --- a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f +++ b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f @@ -90,7 +90,7 @@ subroutine htcdag_bi_ortho_calc_tdav_slow(v, u, N_st, sze) end -subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_chi_array,i_H_phi_array) +subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,chi_H_i_array,i_H_phi_array) use bitmasks implicit none BEGIN_DOC @@ -116,7 +116,7 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) integer(bit_kind), intent(in) :: key(Nint,2) double precision, intent(in) :: coef_l(Ndet_max,Nstate),coef_r(Ndet_max,Nstate) - double precision, intent(out) :: i_H_chi_array(0:3,Nstate),i_H_phi_array(0:3,Nstate) + double precision, intent(out) :: chi_H_i_array(0:3,Nstate),i_H_phi_array(0:3,Nstate) integer :: i, ii,j double precision :: phase @@ -131,7 +131,7 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c ASSERT (Ndet_max >= Ndet) allocate(idx(0:Ndet)) - i_H_chi_array = 0.d0 + chi_H_i_array = 0.d0 i_H_phi_array = 0.d0 call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) @@ -142,10 +142,10 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c ! computes !DIR$ FORCEINLINE call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot) - i_H_chi_array(0,1) = i_H_chi_array(0,1) + coef_l(i,1)*htot - i_H_chi_array(1,1) = i_H_chi_array(1,1) + coef_l(i,1)*hmono - i_H_chi_array(2,1) = i_H_chi_array(2,1) + coef_l(i,1)*htwoe - i_H_chi_array(3,1) = i_H_chi_array(3,1) + coef_l(i,1)*hthree + chi_H_i_array(0,1) = chi_H_i_array(0,1) + coef_l(i,1)*htot + chi_H_i_array(1,1) = chi_H_i_array(1,1) + coef_l(i,1)*hmono + chi_H_i_array(2,1) = chi_H_i_array(2,1) + coef_l(i,1)*htwoe + chi_H_i_array(3,1) = chi_H_i_array(3,1) + coef_l(i,1)*hthree ! computes !DIR$ FORCEINLINE call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot) @@ -163,10 +163,10 @@ subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_c !DIR$ FORCEINLINE call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot) do j = 1, Nstate - i_H_chi_array(0,j) = i_H_chi_array(0,j) + coef_l(i,j)*htot - i_H_chi_array(1,j) = i_H_chi_array(1,j) + coef_l(i,j)*hmono - i_H_chi_array(2,j) = i_H_chi_array(2,j) + coef_l(i,j)*htwoe - i_H_chi_array(3,j) = i_H_chi_array(3,j) + coef_l(i,j)*hthree + chi_H_i_array(0,j) = chi_H_i_array(0,j) + coef_l(i,j)*htot + chi_H_i_array(1,j) = chi_H_i_array(1,j) + coef_l(i,j)*hmono + chi_H_i_array(2,j) = chi_H_i_array(2,j) + coef_l(i,j)*htwoe + chi_H_i_array(3,j) = chi_H_i_array(3,j) + coef_l(i,j)*hthree enddo ! computes !DIR$ FORCEINLINE From cc20c97eda6bdf4240c7d689022f37a867bce81a Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 10 Aug 2023 17:07:40 +0200 Subject: [PATCH 18/21] all the one-body gradients are correct for TC-CASSCF --- src/casscf_tc_bi/grad_dm.irp.f | 59 +++++++++++++++++++++++++++++++-- src/casscf_tc_bi/grad_old.irp.f | 22 +++++++----- 2 files changed, 69 insertions(+), 12 deletions(-) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index 1618adc6..00b20d41 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -1,11 +1,29 @@ BEGIN_PROVIDER [real*8, gradvec_tc_r, (0:3,nMonoEx)] &BEGIN_PROVIDER [real*8, gradvec_tc_l, (0:3,nMonoEx)] + BEGIN_DOC +! gradvec_tc_r(0:3,i) = +! +! gradvec_tc_l(0:3,i) = +! +! where the indices "i" corresponds to E_q^p(i) +! +! i = mat_idx_c_a(q,p) +! +! and gradvec_tc_r/l(0) = full matrix element +! +! gradvec_tc_r/l(1) = one-body part + +! gradvec_tc_r/l(2) = two-body part + +! gradvec_tc_r/l(3) = three-body part + END_DOC implicit none integer :: ii,tt,aa,indx integer :: i,t,a,fff double precision :: res_l(0:3), res_r(0:3) gradvec_tc_l = 0.d0 gradvec_tc_r = 0.d0 + ! computing the core/inactive --> virtual orbitals gradients do i=1,n_core_inact_orb ii=list_core_inact(i) do t=1,n_act_orb @@ -33,9 +51,15 @@ end do do t=1,n_act_orb + tt=list_act(t) do a=1,n_virt_orb - indx = mat_idx_a_v(i,a) -! gradvec_tc_l(indx)=gradvec_ta(t,a) + aa=list_virt(a) + indx = mat_idx_a_v(t,a) + call gradvec_tc_ta(tt,aa,res_l, res_r) + do fff = 0,3 + gradvec_tc_l(fff,indx)=res_l(fff) + gradvec_tc_r(fff,indx)=res_r(fff) + enddo end do end do END_PROVIDER @@ -65,7 +89,7 @@ subroutine gradvec_tc_it(i,t,res_l, res_r) ! ! Corresponds to res_r = ! -! res_l = +! res_l = END_DOC integer, intent(in) :: i,t double precision, intent(out) :: res_l(0:3),res_r(0:3) @@ -83,3 +107,32 @@ subroutine gradvec_tc_it(i,t,res_l, res_r) enddo end + +subroutine gradvec_tc_ta(t,a,res_l, res_r) + implicit none + BEGIN_DOC +! active --> virtual TC gradient +! +! Corresponds to res_r = +! +! res_l = + END_DOC + integer, intent(in) :: t,a + double precision, intent(out) :: res_l(0:3),res_r(0:3) + integer :: rr,r,m + double precision :: dm + res_r = 0.d0 + res_l = 0.d0 +! do rr = 1, n_act_orb +! r = list_act(rr) +! res_l(1) += mo_bi_ortho_tc_one_e(a,r) * tc_transition_matrix_mo(t,r,1,1) +! res_r(1) += -mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1) +! enddo + do m = 1, mo_num + res_r(1) += mo_bi_ortho_tc_one_e(t,m) * tc_transition_matrix_mo(a,m,1,1) & + -mo_bi_ortho_tc_one_e(m,a) * tc_transition_matrix_mo(m,t,1,1) + res_l(1) += mo_bi_ortho_tc_one_e(a,m) * tc_transition_matrix_mo(t,m,1,1) & + -mo_bi_ortho_tc_one_e(m,t) * tc_transition_matrix_mo(m,a,1,1) + enddo + +end diff --git a/src/casscf_tc_bi/grad_old.irp.f b/src/casscf_tc_bi/grad_old.irp.f index 6c976d66..e8440513 100644 --- a/src/casscf_tc_bi/grad_old.irp.f +++ b/src/casscf_tc_bi/grad_old.irp.f @@ -38,15 +38,19 @@ enddo enddo enddo -! do indx=1,nMonoEx -! ihole=excit(1,indx) -! ipart=excit(2,indx) -! call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) -! do ll = 0, 3 -! gradvec_detail_left_old (ll,indx)=res_l(ll) -! gradvec_detail_right_old(ll,indx)=res_r(ll) -! enddo -! end do + + do tt = 1, n_act_orb + ihole = list_act(tt) + do aa = 1, n_virt_orb + ipart = list_virt(aa) + indx = mat_idx_a_v(tt,aa) + call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r) + do ll = 0, 3 + gradvec_detail_left_old (ll,indx)=res_l(ll) + gradvec_detail_right_old(ll,indx)=res_r(ll) + enddo + enddo + enddo real*8 :: norm_grad_left, norm_grad_right norm_grad_left=0.d0 From a15055e9648d65fd4c95513ccb9f7ed0b765199c Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 10 Aug 2023 17:26:33 +0200 Subject: [PATCH 19/21] optimization of one-body TC-CASSCF gradients --- src/casscf_tc_bi/grad_dm.irp.f | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index 00b20d41..6e06f6ce 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -93,17 +93,17 @@ subroutine gradvec_tc_it(i,t,res_l, res_r) END_DOC integer, intent(in) :: i,t double precision, intent(out) :: res_l(0:3),res_r(0:3) - integer :: rr,r,ss,s,m + integer :: rr,r,ss,s,m,mm double precision :: dm res_r = 0.d0 - do m = 1, mo_num - res_r(1) += mo_bi_ortho_tc_one_e(i,m) * tc_transition_matrix_mo(t,m,1,1) & - -mo_bi_ortho_tc_one_e(m,t) * tc_transition_matrix_mo(m,i,1,1) - enddo res_l = 0.d0 - do m = 1, mo_num - res_l(1) += mo_bi_ortho_tc_one_e(t,m) * tc_transition_matrix_mo(i,m,1,1) & - -mo_bi_ortho_tc_one_e(m,i) * tc_transition_matrix_mo(m,t,1,1) + res_r(1) += -2.d0 * mo_bi_ortho_tc_one_e(i,t) + res_l(1) += 2.D0 * mo_bi_ortho_tc_one_e(t,i) + + do rr = 1, n_act_orb + r = list_act(rr) + res_r(1) += mo_bi_ortho_tc_one_e(i,r) * tc_transition_matrix_mo(t,r,1,1) + res_l(1) += -mo_bi_ortho_tc_one_e(r,i) * tc_transition_matrix_mo(r,t,1,1) enddo end @@ -123,16 +123,10 @@ subroutine gradvec_tc_ta(t,a,res_l, res_r) double precision :: dm res_r = 0.d0 res_l = 0.d0 -! do rr = 1, n_act_orb -! r = list_act(rr) -! res_l(1) += mo_bi_ortho_tc_one_e(a,r) * tc_transition_matrix_mo(t,r,1,1) -! res_r(1) += -mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1) -! enddo - do m = 1, mo_num - res_r(1) += mo_bi_ortho_tc_one_e(t,m) * tc_transition_matrix_mo(a,m,1,1) & - -mo_bi_ortho_tc_one_e(m,a) * tc_transition_matrix_mo(m,t,1,1) - res_l(1) += mo_bi_ortho_tc_one_e(a,m) * tc_transition_matrix_mo(t,m,1,1) & - -mo_bi_ortho_tc_one_e(m,t) * tc_transition_matrix_mo(m,a,1,1) + do rr = 1, n_act_orb + r = list_act(rr) + res_l(1) += mo_bi_ortho_tc_one_e(a,r) * tc_transition_matrix_mo(t,r,1,1) + res_r(1) += -mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1) enddo end From 5dc4fb29284c496cf3b6921e9c9422e94b972a97 Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 16 Aug 2023 14:06:29 +0200 Subject: [PATCH 20/21] naive two rdm in tc works for He in cisd and bi ortho orbitals --- src/tc_bi_ortho/test_tc_two_rdm.irp.f | 60 ++++++++++++ src/tc_bi_ortho/two_rdm_naive.irp.f | 132 ++++++++++++++++++++++++++ 2 files changed, 192 insertions(+) create mode 100644 src/tc_bi_ortho/test_tc_two_rdm.irp.f create mode 100644 src/tc_bi_ortho/two_rdm_naive.irp.f diff --git a/src/tc_bi_ortho/test_tc_two_rdm.irp.f b/src/tc_bi_ortho/test_tc_two_rdm.irp.f new file mode 100644 index 00000000..ecdeef43 --- /dev/null +++ b/src/tc_bi_ortho/test_tc_two_rdm.irp.f @@ -0,0 +1,60 @@ +program test_tc_rdm + + BEGIN_DOC + ! + ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together + ! with the energy. Saves the left-right wave functions at the end. + ! + END_DOC + + my_grid_becke = .True. + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + read_wf = .True. + touch read_wf + + print*, ' nb of states = ', N_states + print*, ' nb of det = ', N_det + + call test() + +end + +subroutine test + implicit none + integer :: h1,p1,h2,p2,i,j,istate + double precision :: rdm, integral, accu,ref + double precision :: hmono, htwoe, hthree, htot + accu = 0.d0 + do h1 = 1, mo_num + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + rdm = tc_two_rdm(p1,h1,p2,h2) + accu += integral * rdm + enddo + enddo + enddo + enddo + accu *= 0.5d0 + print*,'accu = ',accu +! print*,tc_two_rdm(1,1,1,1),mo_bi_ortho_tc_two_e(1,1,1,1) + ref = 0.d0 + do i = 1, N_det + do j = 1, N_det +! if(i.ne.j)cycle + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot) + do istate = 1,N_states +! print*,'i,j',i,j,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe +! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) , htwoe + ref += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe + enddo + enddo + enddo + print*,' ref = ',ref + +end diff --git a/src/tc_bi_ortho/two_rdm_naive.irp.f b/src/tc_bi_ortho/two_rdm_naive.irp.f new file mode 100644 index 00000000..9694c653 --- /dev/null +++ b/src/tc_bi_ortho/two_rdm_naive.irp.f @@ -0,0 +1,132 @@ +BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)] + implicit none + BEGIN_DOC + ! tc_two_rdm(p,s,q,r) = psi_det(i) + call get_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1) + do istate = 2, N_states + contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate) + enddo + if(degree == 2)then + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + else if(degree==1)then +! cycle + ! occupation of the determinant psi_det(j) + call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) + + ! run over the electrons of opposite spin than the excitation + s2 = other_spin(s1) + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + ! run over the electrons of same spin than the excitation + s2 = s1 + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + if(h2.le.h1)cycle + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + endif + else if(degree == 0)then + contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * state_average_weight(1) +! print*,'contrib',contrib + do istate = 2, N_states + contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) + enddo + ! occupation of the determinant psi_det(j) + call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) + s1 = 1 ! alpha electrons + do nn = 1, n_occ_ab(s1) + h1 = occ(nn,s1) + p1 = occ(nn,s1) + ! run over the couple of alpha-beta electrons + s2 = other_spin(s1) + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + ! run over the couple of alpha-alpha electrons + s2 = s1 + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + if(h2.le.h1)cycle + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + enddo + s1 = 2 + do nn = 1, n_occ_ab(s1) + h1 = occ(nn,s1) + p1 = occ(nn,s1) + ! run over the couple of beta-beta electrons + s2 = s1 + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + if(h2.le.h1)cycle + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + enddo + endif + enddo + enddo + +END_PROVIDER + +subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib) + implicit none + integer, intent(in) :: h1,p1,h2,p2,s1,s2,sze + double precision, intent(in) :: contrib + double precision, intent(inout) :: array(sze, sze, sze, sze) + integer :: istate + if(s1.ne.s2)then + array(p1,h1,p2,h2) += contrib + ! permutation for particle symmetry + array(p2,h2,p1,h1) += contrib + else ! same spin double excitation + array(p1,h1,p2,h2) += contrib + ! exchange + ! exchanging the holes + array(p2,h1,p1,h2) -= contrib + ! exchanging the particles + array(p1,h2,p2,h1) -= contrib + + ! permutation for particle symmetry + array(p2,h2,p1,h1) += contrib + ! exchange + ! exchanging the holes + array(p1,h2,p2,h1) -= contrib + ! exchanging the particles + array(p2,h1,p1,h2) -= contrib + endif + +end From 0ffaf820a2bc156e021e0f68bddf119e04ce82b5 Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 16 Aug 2023 18:07:50 +0200 Subject: [PATCH 21/21] tc two rdm works for O CISD in biorthonormal basis --- src/tc_bi_ortho/test_tc_two_rdm.irp.f | 19 ++++++++++++++++--- src/tc_bi_ortho/two_rdm_naive.irp.f | 14 +++----------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/tc_bi_ortho/test_tc_two_rdm.irp.f b/src/tc_bi_ortho/test_tc_two_rdm.irp.f index ecdeef43..3e556312 100644 --- a/src/tc_bi_ortho/test_tc_two_rdm.irp.f +++ b/src/tc_bi_ortho/test_tc_two_rdm.irp.f @@ -35,6 +35,10 @@ subroutine test do p2 = 1, mo_num integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) rdm = tc_two_rdm(p1,h1,p2,h2) +! if(dabs(rdm).gt.1.d-10)then +! print*,h1,p1,h2,p2 +! print*,rdm,integral,rdm*integral +! endif accu += integral * rdm enddo enddo @@ -42,19 +46,28 @@ subroutine test enddo accu *= 0.5d0 print*,'accu = ',accu -! print*,tc_two_rdm(1,1,1,1),mo_bi_ortho_tc_two_e(1,1,1,1) +! print*,mo_bi_ortho_tc_two_e(2,15,2,1) +! print*,mo_bi_ortho_tc_two_e(15,2,2,1) +! print*,mo_bi_ortho_tc_two_e(2,1,2,15) +! print*,mo_bi_ortho_tc_two_e(2,1,15,2) ref = 0.d0 do i = 1, N_det do j = 1, N_det -! if(i.ne.j)cycle +! if(i.eq.j)cycle call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot) do istate = 1,N_states -! print*,'i,j',i,j,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe +! print*,'i,j',i,j +! print*,psi_l_coef_bi_ortho(i,istate) , psi_r_coef_bi_ortho(j,istate) , htwoe +! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe +! if(i.ne.j)then ! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) , htwoe +! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe +! endif ref += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe enddo enddo enddo print*,' ref = ',ref + print*,'delta= ',ref-accu end diff --git a/src/tc_bi_ortho/two_rdm_naive.irp.f b/src/tc_bi_ortho/two_rdm_naive.irp.f index 9694c653..8fd34975 100644 --- a/src/tc_bi_ortho/two_rdm_naive.irp.f +++ b/src/tc_bi_ortho/two_rdm_naive.irp.f @@ -30,7 +30,6 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)] if(degree == 2)then call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) else if(degree==1)then -! cycle ! occupation of the determinant psi_det(j) call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) @@ -48,13 +47,12 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)] m = occ(mm,s2) h2 = m p2 = m - if(h2.le.h1)cycle call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) enddo endif else if(degree == 0)then +! cycle contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * state_average_weight(1) -! print*,'contrib',contrib do istate = 2, N_states contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) enddo @@ -115,18 +113,12 @@ subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib) else ! same spin double excitation array(p1,h1,p2,h2) += contrib ! exchange - ! exchanging the holes - array(p2,h1,p1,h2) -= contrib ! exchanging the particles + array(p2,h1,p1,h2) -= contrib + ! exchanging the array(p1,h2,p2,h1) -= contrib - ! permutation for particle symmetry array(p2,h2,p1,h1) += contrib - ! exchange - ! exchanging the holes - array(p1,h2,p2,h1) -= contrib - ! exchanging the particles - array(p2,h1,p1,h2) -= contrib endif end