From 2bda93a97b74f15ca670c11f4a817343aec3461f Mon Sep 17 00:00:00 2001 From: Emmanuel Giner LCT Date: Sat, 26 Oct 2019 16:20:43 +0200 Subject: [PATCH] pert 2 rdm works for a single determinant --- src/casscf/50.casscf.bats | 49 ++++++++++++++++++++++++++++++ src/casscf/test_pert_2rdm.irp.f | 15 ++++++++- src/cipsi/pert_rdm_providers.irp.f | 37 ++++++++++++++++++++-- src/cipsi/update_2rdm.irp.f | 13 ++++---- 4 files changed, 104 insertions(+), 10 deletions(-) create mode 100644 src/casscf/50.casscf.bats diff --git a/src/casscf/50.casscf.bats b/src/casscf/50.casscf.bats new file mode 100644 index 00000000..a0db725d --- /dev/null +++ b/src/casscf/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/test_pert_2rdm.irp.f b/src/casscf/test_pert_2rdm.irp.f index 7c40de0f..890e8cbe 100644 --- a/src/casscf/test_pert_2rdm.irp.f +++ b/src/casscf/test_pert_2rdm.irp.f @@ -2,7 +2,13 @@ program test_pert_2rdm implicit none read_wf = .True. touch read_wf -!call get_pert_2rdm + +!provide is_pert_2rdm_provided + call test + end + + subroutine test + implicit none integer :: i,j,k,l,ii,jj,kk,ll double precision :: accu , get_two_e_integral, integral accu = 0.d0 @@ -25,5 +31,12 @@ program test_pert_2rdm enddo enddo enddo + do ii = 1, n_orb_pert_rdm + i = list_orb_pert_rdm(ii) + do jj = 1, n_orb_pert_rdm + j = list_orb_pert_rdm(jj) + accu += pert_1rdm_provider(j,i) * mo_one_e_integrals(j,i) + enddo + enddo print*,'accu = ',accu end diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index 85bea747..5a18e683 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -10,7 +10,23 @@ END_PROVIDER BEGIN_PROVIDER [logical , pert_2rdm ] implicit none - pert_2rdm = .False. + pert_2rdm = .True. +END_PROVIDER + +BEGIN_PROVIDER [logical, is_pert_2rdm_provided ] + implicit none + is_pert_2rdm_provided = .True. + provide pert_2rdm_provider + if(.True.)then + + double precision :: pt2(N_states),relative_error, error(N_states),variance(N_states),norm(N_states) + relative_error = 0.d0 + pert_2rdm_provider = 0.d0 + call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & + norm,0) ! Stochastic PT2 + print*,'is_pert_2rdm_provided = ',is_pert_2rdm_provided + print*,'pt2 = ',pt2 + endif END_PROVIDER BEGIN_PROVIDER [integer, n_orb_pert_rdm] @@ -30,9 +46,24 @@ BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)] END_PROVIDER -BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)] +BEGIN_PROVIDER [double precision, pert_1rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm)] + implicit none + integer :: i,j,k,l + if(is_pert_2rdm_provided)then + pert_1rdm_provider = 0.d0 + do i = 1, n_orb_pert_rdm + do j = 1, n_orb_pert_rdm + do k = 1, n_orb_pert_rdm + pert_1rdm_provider(j,i) += 2.d0 * pert_2rdm_provider(i,k,j,k) + enddo + enddo + enddo + endif +END_PROVIDER + + BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)] +&BEGIN_PROVIDER [double precision, pert_1rdm_provider_bis, (n_orb_pert_rdm,n_orb_pert_rdm)] implicit none - pert_2rdm_provider = 0.d0 END_PROVIDER diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f index 260c48fd..900be8e7 100644 --- a/src/cipsi/update_2rdm.irp.f +++ b/src/cipsi/update_2rdm.irp.f @@ -28,12 +28,12 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio endif call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff) else - !! case of double excitations - ! if (nkeys + 4 .ge. sze_buff)then - ! call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) - ! nkeys = 0 - ! endif - ! call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) + ! case of double excitations + if (nkeys + 4 .ge. sze_buff)then + call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) + nkeys = 0 + endif + call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) endif enddo !call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) @@ -69,6 +69,7 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v ispin = 2 other_spin = 1 endif + print*,'h1,p1,s1',h1,p1,ispin if(list_orb_reverse_pert_rdm(h1).lt.0)return h1 = list_orb_reverse_pert_rdm(h1) if(list_orb_reverse_pert_rdm(p1).lt.0)return