LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 20:35:19 +01:00
Code Issues Releases Wiki Activity

Merge pull request #119 from scemama/master

Browse Source 
Acceleration of PT2 and extra basis sets
This commit is contained in:
Anthony Scemama 2015-11-24 17:36:03 +01:00
commit 9fbc5afc50
11 changed files with 618 additions and 463 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
ocaml/qp_create_ezfio_from_xyz.ml
View File

@ -91,7 +91,10 @@ let run ?o b c d m p xyz_file =
in
let basis_table = Hashtbl.Poly.create () in
let basis_table =
Hashtbl.Poly.create ()
in
(* Open basis set channels *)
let basis_channel element =
let key =
@ -115,31 +118,46 @@ let run ?o b c d m p xyz_file =
Sys.remove temp_filename
in
let fetch_channel basis =
let command =
if (p) then
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "." ^ basis ^ "\" \"" ^ basis ^"\" pseudo"
else
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "." ^ basis ^ "\" \"" ^ basis ^"\""
in
match Sys.is_file basis with
| `Yes ->
In_channel.create basis
| _ ->
begin
let filename =
Unix.open_process_in command
|> In_channel.input_all
|> String.strip
in
let new_channel =
In_channel.create filename
in
Unix.unlink filename;
new_channel
end
in
let rec build_basis = function
| [] -> ()
| elem_and_basis_name :: rest ->
begin
match (String.lsplit2 ~on:':' elem_and_basis_name) with
| None -> (* Principal basis *)
let basis = elem_and_basis_name in
let command =
if (p) then
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "." ^ basis ^ "\" \"" ^ basis ^"\" pseudo"
else
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "." ^ basis ^ "\" \"" ^ basis ^"\""
in
begin
let filename =
Unix.open_process_in command
|> In_channel.input_all
|> String.strip
let basis =
elem_and_basis_name
in
let new_channel =
In_channel.create filename
fetch_channel basis
in
Unix.unlink filename;
List.iter nuclei ~f:(fun elem->
let key =
Element.to_string elem.Atom.element
@ -151,26 +169,18 @@ let run ?o b c d m p xyz_file =
end
| Some (key, basis) -> (*Aux basis *)
begin
let elem = Element.of_string key
and basis = String.lowercase basis
let elem =
Element.of_string key
and basis =
String.lowercase basis
in
let key =
Element.to_string elem
in
let command =
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "." ^ basis ^ "\" \"" ^ basis ^ "\" "
let new_channel =
fetch_channel basis
in
begin
let filename =
Unix.open_process_in command
|> In_channel.input_all
|> String.strip
in
let new_channel =
In_channel.create filename
in
Unix.unlink filename;
match Hashtbl.add basis_table ~key:key ~data:new_channel with
| `Ok -> ()
| `Duplicate -> failwith ("Duplicate definition of basis for "^(Element.to_long_string elem))
@ -335,13 +345,15 @@ let command =
Command.basic
~summary: "Quantum Package command"
~readme:(fun () -> "
Creates an EZFIO directory from a standard xyz file.
The basis set is defined as a single string if all the
atoms are taken from the same basis set, otherwise specific
elements can be defined as follows:
Creates an EZFIO directory from a standard xyz file. The basis set is defined
as a single string if all the atoms are taken from the same basis set,
otherwise specific elements can be defined as follows:
-b \"cc-pcvdz | H:cc-pvdz | C:6-31g\"
If a file with the same name as the basis set exists, this file will be read.
Otherwise, the basis set is obtained from the database.
")
spec
(fun o b c d m p xyz_file () ->

50
plugins/Perturbation/Moller_plesset.irp.f
View File

@ -1,50 +0,0 @@
subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/(difference of orbital energies)
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies)
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase,delta_e,h
integer :: h1,h2,p1,p2,s1,s2
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
if (degree == 2) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
(Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
delta_e = 1.d0/delta_e
else if (degree == 1) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)
delta_e = 1.d0/delta_e
else
delta_e = 0.d0
endif
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
H_pert_diag(i) = h
c_pert(i) = i_H_psi_array(i) *delta_e
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo
end

106
plugins/Perturbation/epstein_nesbet.irp.f
View File

@ -1,106 +0,0 @@
subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer, intent(in) :: N_minilist
integer, intent(in) :: idx_minilist(0:N_det_selectors)
integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors)
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem, h
PROVIDE selection_criterion
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = h
endif
enddo
end
subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer, intent(in) :: N_minilist
integer, intent(in) :: idx_minilist(0:N_det_selectors)
integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various N_st states.
!
! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
!
! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,delta_e, h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
PROVIDE CI_electronic_energy
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
if (i_H_psi_array(i) /= 0.d0) then
delta_e = h - CI_electronic_energy(i)
if (delta_e > 0.d0) then
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
else
e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
endif
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
endif
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
else
e_2_pert(i) = 0.d0
c_pert(i) = 0.d0
H_pert_diag(i) = 0.d0
endif
enddo
end

214
plugins/Perturbation/pert_sc2.irp.f
View File

@ -1,166 +1,3 @@
subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
integer, intent(in) :: N_minilist
integer, intent(in) :: idx_minilist(0:N_det_selectors)
integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! In addition, for the perturbative energetic contribution you have the standard second order
!
! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
!
! and also the purely projected contribution
!
! H_pert_diag = <HF|H|det_pert> c_pert
END_DOC
integer :: i,j,degree,l
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
!$IVDEP
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
e_2_pert(i) = -dabs(i_H_psi_array(i))
endif
enddo
degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + &
popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
!DEC$ NOUNROLL
do l=2,Nint
degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + &
popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
enddo
if(degree==4)then
! <psi|delta_H|psi>
e_2_pert_fonda = e_2_pert(1)
H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1)
do i = 1, N_st
do j = 1, idx_repeat(0)
e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i)
enddo
enddo
endif
end
subroutine pt2_epstein_nesbet_SC2_no_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
integer, intent(in) :: N_minilist
integer, intent(in) :: idx_minilist(0:N_det_selectors)
integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! In addition, for the perturbative energetic contribution you have the standard second order
!
! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
!
! and also the purely projected contribution
!
! H_pert_diag = <HF|H|det_pert> c_pert
END_DOC
integer :: i,j,degree,l
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
!$IVDEP
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
e_2_pert(i) = -dabs(i_H_psi_array(i))
endif
enddo
end
double precision function repeat_all_e_corr(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
@ -190,54 +27,3 @@ double precision function repeat_all_e_corr(key_in)
end
subroutine pt2_epstein_nesbet_sc2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer, intent(in) :: N_minilist
integer, intent(in) :: idx_minilist(0:N_det_selectors)
integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors)
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states, but with the CISD_SC2 energies and coefficients
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem, h
PROVIDE selection_criterion
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = h
endif
enddo
end

39
plugins/Perturbation/perturbation.template.f
View File

@ -2,7 +2,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
import perturbation
END_SHELL
subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask)
subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
implicit none
BEGIN_DOC
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
@ -12,6 +12,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
integer, intent(in) :: Nint, N_st, buffer_size, i_generator
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
integer(bit_kind),intent(in) :: key_mask(Nint,2)
double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
@ -43,24 +44,8 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
end if
buffer_loop : do i = 1,buffer_size
do i=1,buffer_size
! do k=1,N_minilist_gen
! ex = 0
! do ni=1,Nint
! ex += popcnt(xor(minilist_gen(ni,1,k), buffer(ni,1,i))) + popcnt(xor(minilist_gen(ni,2,k), buffer(ni,2,i)))
! end do
! if(ex <= 4) then
! cycle buffer_loop
! end if
! end do
! c_ref = connected_to_ref(buffer(1,1,i),miniList_gen,Nint,N_minilist_gen+1,N_minilist_gen)
!
! if (c_ref /= 0) then
! cycle
! endif
if(is_connected_to(buffer(1,1,i), miniList_gen, Nint, N_minilist_gen)) then
cycle
end if
@ -71,20 +56,18 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
integer :: degree
call get_excitation_degree(HF_bitmask,buffer(1,1,i),degree,N_int)
! call pt2_$PERT(buffer(1,1,i), &
! c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st,minilist,idx_minilist)
call pt2_$PERT(buffer(1,1,i), &
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) !!!!!!!!!!!!!!!!! MAUVAISE SIGNATURE PR LES AUTRES PT2_* !!!!!
call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
do k = 1,N_st
e_2_pert_buffer(k,i) = e_2_pert(k)
coef_pert_buffer(k,i) = c_pert(k)
sum_norm_pert(k) += c_pert(k) * c_pert(k)
sum_e_2_pert(k) += e_2_pert(k)
sum_H_pert_diag(k) += H_pert_diag(k)
e_2_pert_buffer(k,i) = e_2_pert(k)
coef_pert_buffer(k,i) = c_pert(k)
sum_norm_pert(k) = sum_norm_pert(k) + c_pert(k) * c_pert(k)
sum_e_2_pert(k) = sum_e_2_pert(k) + e_2_pert(k)
sum_H_pert_diag(k) = sum_H_pert_diag(k) + H_pert_diag(k)
enddo
enddo buffer_loop
enddo
end

367
plugins/Perturbation/pt2_equations.irp.f Normal file
View File

@ -0,0 +1,367 @@
BEGIN_TEMPLATE
subroutine pt2_epstein_nesbet ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem_fock, h
double precision :: i_H_psi_array(N_st)
PROVIDE selection_criterion
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = h
endif
enddo
end
subroutine pt2_epstein_nesbet_2x2 ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various N_st states.
!
! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
!
! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem_fock,delta_e, h
double precision :: i_H_psi_array(N_st)
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
PROVIDE CI_electronic_energy
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if (i_H_psi_array(i) /= 0.d0) then
delta_e = h - CI_electronic_energy(i)
if (delta_e > 0.d0) then
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
else
e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
endif
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
endif
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
else
e_2_pert(i) = 0.d0
c_pert(i) = 0.d0
H_pert_diag(i) = 0.d0
endif
enddo
end
subroutine pt2_moller_plesset ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/(difference of orbital energies)
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies)
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem_fock
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase,delta_e,h
double precision :: i_H_psi_array(N_st)
integer :: h1,h2,p1,p2,s1,s2
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
if (degree == 2) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
(Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
delta_e = 1.d0/delta_e
else if (degree == 1) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)
delta_e = 1.d0/delta_e
else
delta_e = 0.d0
endif
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,n_st
H_pert_diag(i) = h
c_pert(i) = i_H_psi_array(i) *delta_e
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo
end
subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! In addition, for the perturbative energetic contribution you have the standard second order
!
! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
!
! and also the purely projected contribution
!
! H_pert_diag = <HF|H|det_pert> c_pert
END_DOC
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
integer :: i,j,degree,l
double precision :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
!$IVDEP
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
h = h + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
e_2_pert(i) = -dabs(i_H_psi_array(i))
endif
enddo
degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + &
popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
!DEC$ NOUNROLL
do l=2,Nint
degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + &
popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
enddo
if(degree==4)then
! <psi|delta_H|psi>
e_2_pert_fonda = e_2_pert(1)
H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1)
do i = 1, N_st
do j = 1, idx_repeat(0)
e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i)
enddo
enddo
endif
end
subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
!
! that could be repeated to this determinant.
!
! In addition, for the perturbative energetic contribution you have the standard second order
!
! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
!
! and also the purely projected contribution
!
! H_pert_diag = <HF|H|det_pert> c_pert
END_DOC
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
integer :: i,j,degree,l
double precision :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
!$IVDEP
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
h = h + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
e_2_pert(i) = -dabs(i_H_psi_array(i))
endif
enddo
end
subroutine pt2_epstein_nesbet_sc2 ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states, but with the CISD_SC2 energies and coefficients
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: i_H_psi_array(N_st)
double precision :: diag_H_mat_elem_fock, h
PROVIDE selection_criterion
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = h
endif
enddo
end
SUBST [ arguments, declarations ]
det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
integer, intent(in) :: Nint
integer, intent(in) :: ndet
integer, intent(in) :: N_st
integer, intent(in) :: N_minilist
integer(bit_kind), intent(in) :: det_ref (Nint,2)
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
double precision , intent(out) :: c_pert(N_st)
double precision , intent(out) :: e_2_pert(N_st)
double precision, intent(out) :: H_pert_diag(N_st)
integer, intent(in) :: idx_minilist(0:N_det_selectors)
integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors)
;;
END_TEMPLATE
! Note : If the arguments are changed here, they should also be changed accordingly in
! the perturbation.template.f file.

11
scripts/generate_h_apply.py
View File

@ -131,10 +131,10 @@ class H_apply(object):
def filter_vvvv_excitation(self):
self["filter_vvvv_excitation"] = """
key_union_hole_part = 0_bit_kind
call set_bite_to_integer(i_a,key_union_hole_part,N_int)
call set_bite_to_integer(j_a,key_union_hole_part,N_int)
call set_bite_to_integer(i_b,key_union_hole_part,N_int)
call set_bite_to_integer(j_b,key_union_hole_part,N_int)
call set_bit_to_integer(i_a,key_union_hole_part,N_int)
call set_bit_to_integer(j_a,key_union_hole_part,N_int)
call set_bit_to_integer(i_b,key_union_hole_part,N_int)
call set_bit_to_integer(j_b,key_union_hole_part,N_int)
do jtest_vvvv = 1, N_int
if(iand(key_union_hole_part(jtest_vvvv),virt_bitmask(jtest_vvvv,1).ne.key_union_hole_part(jtest_vvvv)))then
b_cycle = .False.
@ -157,7 +157,6 @@ class H_apply(object):
def set_filter_2h_2p(self):
self["filter2h2p"] = """
! ! DIR$ FORCEINLINE
if (is_a_two_holes_two_particles(key)) cycle
"""
@ -205,7 +204,7 @@ class H_apply(object):
"""
self.data["keys_work"] = """
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask)
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp)
"""%(pert,)
self.data["finalization"] = """
"""

84
src/Determinants/Fock_diag.irp.f Normal file
View File

@ -0,0 +1,84 @@
subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Build the diagonal of the Fock matrix corresponding to a generator
! determinant. F_00 is <i|H|i> = E0.
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det_ref(Nint,2)
double precision, intent(out) :: fock_diag_tmp(2,mo_tot_num+1)
integer :: occ(Nint*bit_kind_size,2)
integer :: ne(2), i, j, ii, jj
double precision :: E0
! Compute Fock matrix diagonal elements
call bitstring_to_list_ab(det_ref,occ,Ne,Nint)
fock_diag_tmp = 0.d0
E0 = 0.d0
! Occupied MOs
do ii=1,elec_alpha_num
i = occ(ii,1)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i)
E0 = E0 + mo_mono_elec_integral(i,i)
do jj=1,elec_alpha_num
j = occ(jj,1)
if (i==j) cycle
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j)
E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j)
enddo
do jj=1,elec_beta_num
j = occ(jj,2)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j)
E0 = E0 + mo_bielec_integral_jj(i,j)
enddo
enddo
do ii=1,elec_beta_num
i = occ(ii,2)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i)
E0 = E0 + mo_mono_elec_integral(i,i)
do jj=1,elec_beta_num
j = occ(jj,2)
if (i==j) cycle
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j)
E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j)
enddo
do jj=1,elec_alpha_num
j = occ(jj,1)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j)
enddo
enddo
! Virtual MOs
do i=1,mo_tot_num
if (fock_diag_tmp(1,i) /= 0.d0) cycle
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i)
do jj=1,elec_alpha_num
j = occ(jj,1)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j)
enddo
do jj=1,elec_beta_num
j = occ(jj,2)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j)
enddo
enddo
do i=1,mo_tot_num
if (fock_diag_tmp(2,i) /= 0.d0) cycle
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i)
do jj=1,elec_beta_num
j = occ(jj,2)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j)
enddo
do jj=1,elec_alpha_num
j = occ(jj,1)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j)
enddo
enddo
fock_diag_tmp(1,mo_tot_num+1) = E0
fock_diag_tmp(2,mo_tot_num+1) = E0
end

57
src/Determinants/H_apply.template.f
View File

@ -1,6 +1,6 @@
subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters )
subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
integer(bit_kind), intent(in) :: key_in(N_int, 2), hole_1(N_int, 2), hole_2(N_int, 2)
integer(bit_kind), intent(in) :: particl_1(N_int, 2), particl_2(N_int, 2)
@ -8,8 +8,8 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
integer,intent(in) :: i_generator,iproc_in
integer(bit_kind) :: status(N_int*bit_kind_size, 2)
integer :: highest, p1,p2,sp,ni,i,mi,nt,ns
integer(bit_kind), intent(in) :: key_prev(N_int, 2, *)
double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
integer(bit_kind), intent(in) :: key_prev(N_int, 2, *)
PROVIDE N_int
PROVIDE N_det
@ -72,7 +72,7 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
if((status(p1, sp) == 1 .and. status(p2, sp) > 1) .or. &
(status(p1, sp) == 2 .and. status(p2, sp) == 3) .or. &
(status(p1, sp) == 3 .and. status(p2, sp) == 3 .and. p2 > p1)) then
call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, i_generator, iproc_in $parameters )
call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
end if
end do
end do
@ -89,16 +89,17 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
(status(p1, 1) == 1 .and. status(p2, 2) >= 2) .or. &
(status(p1, 1) == 2 .and. status(p2, 2) /= 2)) then
call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, i_generator, iproc_in $parameters )
call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
end if
end do
end do
end subroutine
subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, i_generator, iproc_in $parameters )
subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
integer(bit_kind), intent(in) :: key_in(N_int, 2), particl_1(N_int, 2), particl_2(N_int, 2)
double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
integer(bit_kind) :: p1_mask(N_int, 2), p2_mask(N_int, 2), key_mask(N_int, 2)
integer,intent(in) :: fh1,fh2,fs1,fs2,i_generator,iproc_in
integer(bit_kind) :: miniList(N_int, 2, N_det)
@ -115,11 +116,11 @@ subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2,
key_mask(ishft(fh1,-bit_kind_shift) + 1, fs1) -= ishft(1,iand(fh1-1,bit_kind_size-1))
key_mask(ishft(fh2,-bit_kind_shift) + 1, fs2) -= ishft(1,iand(fh2-1,bit_kind_size-1))
call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, i_generator, iproc_in $parameters )
call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
end subroutine
subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, i_generator, iproc_in $parameters )
subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
use omp_lib
use bitmasks
implicit none
@ -136,6 +137,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2)
integer, intent(in) :: iproc_in
double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
integer(bit_kind), allocatable :: hole_save(:,:)
integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:)
integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:)
@ -175,6 +177,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), &
occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),&
occ_hole_tmp(N_int*bit_kind_size,2),key_union_hole_part(N_int))
$init_thread
@ -362,17 +365,17 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
enddo ! ispin
$keys_work
$deinit_thread
deallocate (ia_ja_pairs, ib_jb_pairs, &
keys_out, hole_save, &
key,hole, particle, hole_tmp,&
particle_tmp, occ_particle, &
occ_hole, occ_particle_tmp,&
deallocate (ia_ja_pairs, ib_jb_pairs, &
keys_out, hole_save, &
key,hole, particle, hole_tmp, &
particle_tmp, occ_particle, &
occ_hole, occ_particle_tmp, &
occ_hole_tmp,array_pairs,key_union_hole_part)
$omp_end_parallel
$finalization
end
subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $parameters )
subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generator,iproc_in $parameters )
use omp_lib
use bitmasks
implicit none
@ -387,6 +390,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa
integer(bit_kind),intent(in) :: key_in(N_int,2)
integer(bit_kind),intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer, intent(in) :: iproc_in
double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
integer(bit_kind),allocatable :: keys_out(:,:,:)
integer(bit_kind),allocatable :: hole_save(:,:)
integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:)
@ -526,6 +530,7 @@ subroutine $subroutine($params_main)
integer(bit_kind), allocatable :: mask(:,:,:)
integer :: ispin, k
integer :: iproc
double precision, allocatable :: fock_diag_tmp(:,:)
$initialization
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators
@ -539,7 +544,7 @@ subroutine $subroutine($params_main)
call wall_time(wall_0)
iproc = 0
allocate( mask(N_int,2,6) )
allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) )
do i_generator=1,nmax
progress_bar(1) = i_generator
@ -549,6 +554,9 @@ subroutine $subroutine($params_main)
endif
$skip
! Compute diagonal of the Fock matrix
call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
! Create bit masks for holes and particles
do ispin=1,2
do k=1,N_int
@ -577,12 +585,12 @@ subroutine $subroutine($params_main)
psi_det_generators(1,1,1), &
mask(1,1,d_hole1), mask(1,1,d_part1), &
mask(1,1,d_hole2), mask(1,1,d_part2), &
i_generator, iproc $params_post)
fock_diag_tmp, i_generator, iproc $params_post)
endif
if($do_mono_excitations)then
call $subroutine_monoexc(psi_det_generators(1,1,i_generator), &
mask(1,1,s_hole ), mask(1,1,s_part ), &
i_generator, iproc $params_post)
fock_diag_tmp, i_generator, iproc $params_post)
endif
call wall_time(wall_1)
$printout_always
@ -592,13 +600,13 @@ subroutine $subroutine($params_main)
endif
enddo
deallocate( mask )
deallocate( mask, fock_diag_tmp )
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc)
!$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc,fock_diag_tmp)
call wall_time(wall_0)
!$ iproc = omp_get_thread_num()
allocate( mask(N_int,2,6) )
allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) )
!$OMP DO SCHEDULE(dynamic,1)
do i_generator=nmax+1,N_det_generators
if (iproc == 0) then
@ -609,6 +617,9 @@ subroutine $subroutine($params_main)
endif
$skip
! Compute diagonal of the Fock matrix
call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
! Create bit masks for holes and particles
do ispin=1,2
do k=1,N_int
@ -638,12 +649,12 @@ subroutine $subroutine($params_main)
psi_det_generators(1,1,1), &
mask(1,1,d_hole1), mask(1,1,d_part1), &
mask(1,1,d_hole2), mask(1,1,d_part2), &
i_generator, iproc $params_post)
fock_diag_tmp, i_generator, iproc $params_post)
endif
if($do_mono_excitations)then
call $subroutine_monoexc(psi_det_generators(1,1,i_generator), &
mask(1,1,s_hole ), mask(1,1,s_part ), &
i_generator, iproc $params_post)
fock_diag_tmp, i_generator, iproc $params_post)
endif
!$ call omp_set_lock(lck)
call wall_time(wall_1)
@ -655,7 +666,7 @@ subroutine $subroutine($params_main)
!$ call omp_unset_lock(lck)
enddo
!$OMP END DO
deallocate( mask )
deallocate( mask, fock_diag_tmp )
!$OMP END PARALLEL
!$ call omp_destroy_lock(lck)

2
src/Determinants/create_excitations.irp.f
View File

@ -35,7 +35,7 @@ subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok)
endif
end
subroutine set_bite_to_integer(i_physical,key,Nint)
subroutine set_bit_to_integer(i_physical,key,Nint)
use bitmasks
implicit none
integer, intent(in) :: i_physical,Nint

71
src/Determinants/slater_rules.irp.f
View File

@ -1264,6 +1264,75 @@ end
double precision function diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes <i|H|i> when i is at most a double excitation from
! a reference.
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: det_ref(Nint,2), det_pert(Nint,2)
double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
integer :: degree
double precision :: phase, E0
integer :: exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
call get_excitation_degree(det_ref,det_pert,degree,Nint)
E0 = fock_diag_tmp(1,mo_tot_num+1)
if (degree == 2) then
call get_double_excitation(det_ref,det_pert,exc,phase,Nint)
call decode_exc(exc,2,h1,p1,h2,p2,s1,s2)
if ( (s1 == 1).and.(s2 == 1) ) then ! alpha/alpha
diag_H_mat_elem_fock = E0 &
- fock_diag_tmp(1,h1) &
+ ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) &
- ( fock_diag_tmp(1,h2) - mo_bielec_integral_jj_anti(h1,h2) &
+ mo_bielec_integral_jj_anti(p1,h2) ) &
+ ( fock_diag_tmp(1,p2) - mo_bielec_integral_jj_anti(h1,p2) &
+ mo_bielec_integral_jj_anti(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) )
else if ( (s1 == 2).and.(s2 == 2) ) then ! beta/beta
diag_H_mat_elem_fock = E0 &
- fock_diag_tmp(2,h1) &
+ ( fock_diag_tmp(2,p1) - mo_bielec_integral_jj_anti(h1,p1) ) &
- ( fock_diag_tmp(2,h2) - mo_bielec_integral_jj_anti(h1,h2) &
+ mo_bielec_integral_jj_anti(p1,h2) ) &
+ ( fock_diag_tmp(2,p2) - mo_bielec_integral_jj_anti(h1,p2) &
+ mo_bielec_integral_jj_anti(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) )
else ! alpha/beta
diag_H_mat_elem_fock = E0 &
- fock_diag_tmp(1,h1) &
+ ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) &
- ( fock_diag_tmp(2,h2) - mo_bielec_integral_jj(h1,h2) &
+ mo_bielec_integral_jj(p1,h2) ) &
+ ( fock_diag_tmp(2,p2) - mo_bielec_integral_jj(h1,p2) &
+ mo_bielec_integral_jj(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) )
endif
else if (degree == 1) then
call get_mono_excitation(det_ref,det_pert,exc,phase,Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
if (s1 == 1) then
diag_H_mat_elem_fock = E0 - fock_diag_tmp(1,h1) &
+ ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) )
else
diag_H_mat_elem_fock = E0 - fock_diag_tmp(2,h1) &
+ ( fock_diag_tmp(2,p1) - mo_bielec_integral_jj_anti(h1,p1) )
endif
else if (degree == 0) then
diag_H_mat_elem_fock = E0
else
STOP 'Bug in diag_H_mat_elem_fock'
endif
end
double precision function diag_H_mat_elem(det_in,Nint)
implicit none
BEGIN_DOC
@ -1541,8 +1610,8 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0)
do i=shortcut(sh),shortcut(sh+1)-1
local_threshold = threshold_davidson - dabs(u_0(org_i))
org_i = sort_idx(i)
local_threshold = threshold_davidson - dabs(u_0(org_i))
do j=shortcut(sh),i-1
org_j = sort_idx(j)
if ( dabs(u_0(org_j)) > local_threshold ) then