LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-25 17:54:44 +02:00
Code Releases Activity

Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

Browse Source
This commit is contained in:
eginer 2025-03-11 18:37:12 +01:00
commit 8770bf2bf3
17 changed files with 1386 additions and 1057 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/irpf90 vendored

@ -1 +1 @@
Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6
Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b

12
plugins/local/basis_correction/print_routine.irp.f
View File

@ -22,7 +22,10 @@ subroutine print_basis_correction
print*, '****************************************'
print*, '****************************************'
print*, 'mu_of_r_potential = ',mu_of_r_potential
if(mu_of_r_potential.EQ."hf".or.mu_of_r_potential.EQ."hf_old".or.mu_of_r_potential.EQ."hf_sparse")then
if(mu_of_r_potential.EQ."hf".or. &
mu_of_r_potential.EQ."hf_old".or.&
mu_of_r_potential.EQ."hf_sparse".or.&
mu_of_r_potential.EQ."proj")then
print*, ''
print*,'Using a HF-like two-body density to define mu(r)'
print*,'This assumes that HF is a qualitative representation of the wave function '
@ -38,7 +41,9 @@ subroutine print_basis_correction
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo
else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
else if(mu_of_r_potential.EQ."cas_full".or. &
mu_of_r_potential.EQ."cas_truncated".or. &
mu_of_r_potential.EQ."pure_act") then
print*, ''
print*,'Using a CAS-like two-body density to define mu(r)'
print*,'This assumes that the CAS is a qualitative representation of the wave function '
@ -74,7 +79,8 @@ subroutine print_basis_correction
print*,''
print*,'**************'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) , state ',istate,' = ',mu_average_prov(istate)
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) [rho ], state ',istate,' = ',mu_average_prov(istate)
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) [rho^2], state ',istate,' = ',mu_average_prov2(istate)
enddo
end

1
scripts/qp_cipsi_rsh Symbolic link
View File

@ -0,0 +1 @@
/home/scemama/qp2/plugins/qp_plugins_lct/stable/rsdft_cipsi/qp_cipsi_rsh

1
scripts/qp_cipsi_rsh_mu_of_r Symbolic link
View File

@ -0,0 +1 @@
/home/scemama/qp2/plugins/qp_plugins_lct/stable/rsdft_cipsi/qp_cipsi_rsh_mu_of_r

151
scripts/qp_geom_opt.py Executable file
View File

@ -0,0 +1,151 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Usage:
qp_geom_opt [-s state] [-r executable] [-f] [-t tolerance] <EZFIO_FILE>
Options:
-s --state=<state> Excited state to optimize
-f --scf Perform an SCF after each geomety change
-r --qp_run=executable Excited state to optimize
-t --tol=tolerance Convergence criterion on the energy
"""
try:
from docopt import docopt
from module_handler import ModuleHandler, get_dict_child
from module_handler import get_l_module_descendant
from qp_path import QP_SRC, QP_PLUGINS, QP_DATA, QP_ROOT
except ImportError:
print("Please check if you have sourced the ${QP_ROOT}/quantum_package.rc")
print("(`source ${QP_ROOT}/quantum_package.rc`)")
print(sys.exit(1))
import numpy as np
import subprocess
from scipy.optimize import minimize
from ezfio import ezfio
import sys
def set_unbuffered_output():
"""Ensure sys.stdout is unbuffered or line-buffered in a portable way."""
if hasattr(sys.stdout, "reconfigure"): # Python 3.7+
sys.stdout.reconfigure(line_buffering=True)
else:
sys.stdout = open(sys.stdout.fileno(), mode='w', buffering=1)
set_unbuffered_output()
def get_energy(file, state, arguments):
"""Compute the energy of the given state by calling Quantum Package."""
if not arguments["--qp_run"]:
raise ValueError("--qp_run option missing")
if arguments["--scf"]:
executable = "scf"
else:
executable = "save_ortho_mos"
result = subprocess.run(f"qp_run {executable} {file} > {file}.energy.out",
shell=True, capture_output=True, text=True, check=True
)
executable = arguments["--qp_run"]
result = subprocess.run( f"qp_run {executable} {file} > {file}.energy.out",
shell=True)
energy = None
with open(f"{file}.energy.out", 'r') as f:
for line in f:
if "Energy of state" in line and f"{state}" in line:
energy = float(line.split()[-1]) # Extracts the energy value
return energy
raise ValueError("Energy not found in Quantum Package output. Update script {sys.argv[0]}")
def set_coordinates(coord):
"""Update the nuclear coordinates in EZFIO."""
ezfio.set_nuclei_nucl_coord(coord)
def get_coordinates():
"""Retrieve the current nuclear coordinates from EZFIO."""
return np.array(ezfio.get_nuclei_nucl_coord())
memo_energy = {}
def energy_function(coord, file, state, arguments):
"""Wrapper for the energy calculation, ensuring coordinates are updated."""
h = np.array_str(coord)
if h in memo_energy:
return memo_energy[h]
set_coordinates(coord)
energy = get_energy(file, state, arguments)
memo_energy[h] = energy
label = ezfio.get_nuclei_nucl_label()
num_atoms = len(label)
coord = coord.reshape(3, num_atoms).T # Reshape into (num_atoms, 3)
coord_angstrom = coord * 0.529177 # Convert atomic units to angstroms
print(num_atoms)
print(f"Energy: {energy:15.10f}")
for i, (x, y, z) in enumerate(coord_angstrom):
print(f"{label[i]:3s} {x:15.8f} {y:15.8f} {z:15.8f}") # Replace 'X' with actual atomic symbols
return energy
def optimize_geometry(file, state, arguments):
"""Perform geometry optimization using SciPy's minimize function."""
x0 = get_coordinates().flatten()
if arguments["--tol"]:
tolerance = float(tol=arguments["--tol"])
else:
tolerance = 1.e-3
result = minimize(energy_function, x0, args=(file, state, arguments),
method='Powell',
tol=tolerance,
options={'xtol': tolerance, 'ftol': tolerance})
# result = minimize(energy_function, x0, args=(file, state, arguments),
# method='BFGS',
# jac=None,
# tol=tolerance,
# options={'eps': 1.e-3})
if result.success:
print("Optimization successful!")
print("Final energy:", result.fun)
print("Optimized coordinates:", result.x)
else:
print("Optimization failed:", result.message)
set_coordinates(result.x) # Store the optimized geometry
return result
def main(arguments):
if arguments["--state"]:
state=arguments["--state"]
else:
state=1
ezfio_filename = arguments["<EZFIO_FILE>"]
ezfio.set_file(ezfio_filename)
optimize_geometry(ezfio_filename, state, arguments)
if __name__ == "__main__":
ARG = docopt(__doc__)
main(ARG)

2
src/ao_two_e_ints/cholesky.irp.f
View File

@ -178,7 +178,7 @@ END_PROVIDER
rank_max = np
! Avoid too large arrays when there are many electrons
if (elec_num > 10) then
rank_max = min(np,20*elec_num*elec_num)
rank_max = min(np,25*elec_num*elec_num)
endif
call mmap_create_d('', (/ ndim8, rank_max /), .False., .True., map)

1
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -54,6 +54,7 @@ double precision function ao_two_e_integral(i, j, k, l)
else if (use_only_lr) then
ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l)
return
else if (do_schwartz_accel(i,j,k,l)) then

37
src/bitmask/bitmasks_routines.irp.f
View File

@ -283,33 +283,16 @@ subroutine print_det_one_dimension(string,Nint)
end
logical function is_integer_in_string(bite,string,Nint)
logical function is_integer_in_string(orb,bitmask,Nint)
use bitmasks
implicit none
integer, intent(in) :: bite,Nint
integer(bit_kind), intent(in) :: string(Nint)
integer(bit_kind) :: string_bite(Nint)
integer :: i,itot,itot_and
character*(2048) :: output(1)
string_bite = 0_bit_kind
call set_bit_to_integer(bite,string_bite,Nint)
itot = 0
itot_and = 0
is_integer_in_string = .False.
!print*,''
!print*,''
!print*,'bite = ',bite
!call bitstring_to_str( output(1), string_bite, Nint )
! print *, trim(output(1))
!call bitstring_to_str( output(1), string, Nint )
! print *, trim(output(1))
do i = 1, Nint
itot += popcnt(string(i))
itot_and += popcnt(ior(string(i),string_bite(i)))
enddo
!print*,'itot,itot_and',itot,itot_and
if(itot == itot_and)then
is_integer_in_string = .True.
endif
!pause
BEGIN_DOC
! Checks is the orbital orb is set to 1 in the bit string
END_DOC
integer, intent(in) :: orb, Nint
integer(bit_kind), intent(in) :: bitmask(Nint)
integer :: j, k
k = ishft(orb-1,-bit_kind_shift)+1
j = orb-ishft(k-1,bit_kind_shift)-1
is_integer_in_string = iand(bitmask(k), ibset(0_bit_kind, j)) /= 0_bit_kind
end

11
src/cas_based_on_top/on_top_cas_prov.irp.f
View File

@ -15,14 +15,17 @@
pure_act_on_top_of_r = 0.d0
do l = 1, n_act_orb
phi_l = act_mos_in_r_array(l,ipoint)
if (dabs(phi_l) < 1.d-12) cycle
do k = 1, n_act_orb
phi_k = act_mos_in_r_array(k,ipoint)
phi_k = act_mos_in_r_array(k,ipoint) * phi_l
if (dabs(phi_k) < 1.d-12) cycle
do j = 1, n_act_orb
phi_j = act_mos_in_r_array(j,ipoint)
phi_j = act_mos_in_r_array(j,ipoint) * phi_k
if (dabs(phi_j) < 1.d-12) cycle
do i = 1, n_act_orb
phi_i = act_mos_in_r_array(i,ipoint)
phi_i = act_mos_in_r_array(i,ipoint) * phi_j
! 1 2 1 2
pure_act_on_top_of_r += act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i * phi_j * phi_k * phi_l
pure_act_on_top_of_r = pure_act_on_top_of_r + act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i !* phi_j * phi_k * phi_l
enddo
enddo
enddo

88
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -8,21 +8,14 @@ BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)]
END_DOC
implicit none
integer :: i, j
double precision :: tmp_array(ao_num), r(3)
integer :: i
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,tmp_array,j) &
!$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points)
!$OMP PRIVATE (i) &
!$OMP SHARED(aos_in_r_array,n_points_final_grid,final_grid_points)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_at_r(r, tmp_array)
do j = 1, ao_num
aos_in_r_array(j,i) = tmp_array(j)
enddo
call give_all_aos_at_r(final_grid_points(1,i), aos_in_r_array(1,i))
enddo
!$OMP END PARALLEL DO
@ -62,25 +55,27 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_gri
implicit none
integer :: i, j, m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
double precision :: r(3)
double precision, allocatable :: aos_grad_array(:,:), aos_array(:)
!$OMP PARALLEL DO &
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) &
!$OMP SHARED(aos_grad_in_r_array,n_points_final_grid,ao_num,final_grid_points)
allocate(aos_grad_array(3,ao_num), aos_array(ao_num))
!$OMP DO
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
call give_all_aos_and_grad_at_r(final_grid_points(1,i),aos_array,aos_grad_array)
do m = 1, 3
do j = 1, ao_num
aos_grad_in_r_array(j,i,m) = aos_grad_array(m,j)
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END DO
deallocate(aos_grad_array,aos_array)
!$OMP END PARALLEL
END_PROVIDER
@ -126,25 +121,25 @@ END_PROVIDER
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: i,j,m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
double precision :: aos_lapl_array(3,ao_num)
!$OMP PARALLEL DO &
double precision, allocatable :: aos_lapl_array(:,:), aos_grad_array(:,:), aos_array(:)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) &
!$OMP PRIVATE (i,aos_array,aos_grad_array,aos_lapl_array,j,m) &
!$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points)
allocate( aos_array(ao_num), aos_grad_array(3,ao_num), aos_lapl_array(3,ao_num))
!$OMP DO
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
call give_all_aos_and_grad_and_lapl_at_r(final_grid_points(1,i),aos_array,aos_grad_array,aos_lapl_array)
do j = 1, ao_num
do m = 1, 3
aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j)
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END DO
deallocate( aos_array, aos_grad_array, aos_lapl_array)
!$OMP END PARALLEL
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)]
@ -189,20 +184,12 @@ END_PROVIDER
BEGIN_DOC
! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point of the EXTRA grid
END_DOC
integer :: i,j
double precision :: aos_array(ao_num), r(3)
integer :: i
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,aos_array,j) &
!$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra)
!$OMP DEFAULT (NONE) PRIVATE (i) &
!$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,final_grid_points_extra)
do i = 1, n_points_extra_final_grid
r(1) = final_grid_points_extra(1,i)
r(2) = final_grid_points_extra(2,i)
r(3) = final_grid_points_extra(3,i)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
aos_in_r_array_extra(j,i) = aos_array(j)
enddo
call give_all_aos_at_r(final_grid_points_extra(1,i),aos_in_r_array_extra(1,i))
enddo
!$OMP END PARALLEL DO
@ -235,25 +222,26 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array_extra, (ao_num,n_points_ext
implicit none
integer :: i, j, m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
double precision, allocatable :: aos_array(:), aos_grad_array(:,:)
!$OMP PARALLEL DO &
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) &
!$OMP PRIVATE (i,j,m,aos_array,aos_grad_array) &
!$OMP SHARED(aos_grad_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra)
allocate(aos_array(ao_num), aos_grad_array(3,ao_num))
!$OMP DO
do i = 1, n_points_extra_final_grid
r(1) = final_grid_points_extra(1,i)
r(2) = final_grid_points_extra(2,i)
r(3) = final_grid_points_extra(3,i)
call give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array)
call give_all_aos_and_grad_at_r(final_grid_points_extra(1,i), aos_array, aos_grad_array)
do m = 1, 3
do j = 1, ao_num
aos_grad_in_r_array_extra(j,i,m) = aos_grad_array(m,j)
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END DO
deallocate(aos_array,aos_grad_array)
!$OMP END PARALLEL
END_PROVIDER

15
src/dft_utils_in_r/mo_in_r.irp.f
View File

@ -21,20 +21,11 @@
BEGIN_DOC
! mos_in_r_array(i,j) = value of the ith mo on the jth grid point
END_DOC
integer :: i,j
double precision :: mos_array(mo_num), r(3)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,mos_array,j) &
integer :: i
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE (i) &
!$OMP SHARED(mos_in_r_array_omp,n_points_final_grid,mo_num,final_grid_points)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_mos_at_r(r,mos_array)
do j = 1, mo_num
mos_in_r_array_omp(j,i) = mos_array(j)
enddo
call give_all_mos_at_r(final_grid_points(1,i),mos_in_r_array_omp(1,i))
enddo
!$OMP END PARALLEL DO
END_PROVIDER

6
src/mo_two_e_ints/EZFIO.cfg
View File

@ -1,3 +1,9 @@
[do_mo_cholesky]
type: logical
doc: Use Cholesky decomposition of MO integrals in CI calculations
interface: ezfio,provider,ocaml
default: False
[io_mo_cholesky]
type: Disk_access
doc: Read/Write |MO| Cholesky integrals from/to disk [ Write | Read | None ]

9
src/mo_two_e_ints/cholesky.irp.f
View File

@ -1,12 +1,3 @@
BEGIN_PROVIDER [ logical, do_mo_cholesky ]
implicit none
BEGIN_DOC
! If True, use Cholesky vectors for MO integrals
END_DOC
do_mo_cholesky = do_ao_cholesky
! do_mo_cholesky = .False.
END_PROVIDER
BEGIN_PROVIDER [ integer, cholesky_mo_num ]
&BEGIN_PROVIDER [ integer, cholesky_mo_num_split, (1:5)]
implicit none

107
src/mu_of_r/mu_of_r_conditions.irp.f
View File

@ -22,15 +22,26 @@
endif
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
if(mu_of_r_potential.EQ."hf")then
do ipoint = 1, n_points_final_grid
mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint)
enddo
else if(mu_of_r_potential.EQ."hf_old")then
do ipoint = 1, n_points_final_grid
mu_of_r_prov(ipoint,istate) = mu_of_r_hf_old(ipoint)
enddo
else if(mu_of_r_potential.EQ."hf_sparse")then
do ipoint = 1, n_points_final_grid
mu_of_r_prov(ipoint,istate) = mu_of_r_hf_sparse(ipoint)
enddo
else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
do ipoint = 1, n_points_final_grid
mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate)
enddo
else if(mu_of_r_potential.EQ."proj")then
do ipoint = 1, n_points_final_grid
mu_of_r_prov(ipoint,istate) = mu_of_r_projector_mo(ipoint)
enddo
else
print*,'you requested the following mu_of_r_potential'
print*,mu_of_r_potential
@ -38,7 +49,6 @@
stop
endif
enddo
enddo
if (write_mu_of_r) then
print*,'Writing mu(r) on disk ...'
@ -201,7 +211,7 @@
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)]
BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)]
implicit none
BEGIN_DOC
! average value of mu(r) weighted with the total one-e density and divided by the number of electrons
@ -223,5 +233,94 @@
enddo
mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate)
enddo
END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_average_prov2, (N_states)]
implicit none
BEGIN_DOC
! average value of mu(r) weighted with square of the total one-e density
!
! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
!
! in the one- and two-body density matrix are excluded
END_DOC
integer :: ipoint,istate
double precision :: weight,density,norm
mu_average_prov2 = 0.d0
do istate = 1, N_states
norm = 0.d0
do ipoint = 1, n_points_final_grid
weight =final_weight_at_r_vector(ipoint)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) &
+ one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
if(mu_of_r_prov(ipoint,istate).gt.1.d+09)cycle
mu_average_prov2(istate) += mu_of_r_prov(ipoint,istate) * weight * density*density
norm = norm + density*density*weight
enddo
mu_average_prov2(istate) = mu_average_prov2(istate) / norm
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_of_r_projector_mo, (n_points_final_grid) ]
implicit none
BEGIN_DOC
! mu(r) computed with the projector onto the atomic basis
! P_B(\mathbf{r},\mathbf{r}') = \sum_{ij} |
! \chi_{i} \rangle \left[S^{-1}\right]_{ij} \langle \chi_{j} |
! \] where $i$ and $j$ denote all atomic orbitals.
END_DOC
double precision, parameter :: factor = dsqrt(2.d0*dacos(-1.d0))
double precision, allocatable :: tmp(:,:)
integer :: ipoint
do ipoint=1,n_points_final_grid
mu_of_r_projector_mo(ipoint) = 0.d0
integer :: i,j
do j=1,n_inact_act_orb
i = list_inact_act(j)
mu_of_r_projector_mo(ipoint) = mu_of_r_projector_mo(ipoint) + &
mos_in_r_array_omp(i,ipoint) * mos_in_r_array_omp(i,ipoint)
enddo
do j=1,n_virt_orb
i = list_virt(j)
mu_of_r_projector_mo(ipoint) = mu_of_r_projector_mo(ipoint) + &
mos_in_r_array_omp(i,ipoint) * mos_in_r_array_omp(i,ipoint)
enddo
enddo
do ipoint=1,n_points_final_grid
! epsilon
mu_of_r_projector_mo(ipoint) = 1.d0/(2.d0*dacos(-1.d0) * mu_of_r_projector_mo(ipoint)**(2.d0/3.d0))
! mu
mu_of_r_projector_mo(ipoint) = 1.d0/dsqrt( 2.d0*mu_of_r_projector_mo(ipoint) )
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_average_proj, (N_states)]
implicit none
BEGIN_DOC
! average value of mu(r) weighted with the total one-e density and divided by the number of electrons
!
! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
!
! in the one- and two-body density matrix are excluded
END_DOC
integer :: ipoint,istate
double precision :: weight,density
do istate = 1, N_states
mu_average_proj(istate) = 0.d0
do ipoint = 1, n_points_final_grid
weight =final_weight_at_r_vector(ipoint)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) &
+ one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
mu_average_proj(istate) += mu_of_r_projector_mo(ipoint) * weight * density
enddo
mu_average_proj(istate) = mu_average_proj(istate) / elec_num_grid_becke(istate)
enddo
END_PROVIDER

1
src/two_body_rdm/act_2_rdm.irp.f
View File

@ -145,6 +145,7 @@
print*,''
print*,'Providing act_2_rdm_spin_trace_mo '
character*(128) :: name_file
PROVIDE all_mo_integrals
name_file = 'act_2_rdm_spin_trace_mo'
ispin = 4
act_2_rdm_spin_trace_mo = 0.d0

164
src/two_rdm_routines/davidson_like_2rdm.irp.f
View File

@ -13,7 +13,7 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
END_DOC
integer, intent(in) :: N_st,sze
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st)
double precision, intent(in) :: u_0(sze,N_st)
integer :: k
@ -50,7 +50,7 @@ subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st)
double precision, intent(in) :: u_t(N_st,N_det)
integer :: k
@ -91,7 +91,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st)
integer(omp_lock_kind) :: lock_2rdm
integer :: i,j,k,l
@ -139,6 +139,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60
sze_buff = sze_buff*100
list_orb_reverse = -1000
do i = 1, norb
list_orb_reverse(list_orb(i)) = i
@ -154,6 +155,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
! Prepare the array of all alpha single excitations
! -------------------------------------------------
double precision, allocatable :: big_array_local(:,:,:,:,:)
PROVIDE N_int nthreads_davidson elec_alpha_num
!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,&
@ -173,7 +176,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
!$OMP buffer, doubles, n_doubles, &
!$OMP tmp_det2, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, &
!$OMP n_singles_b, nkeys, keys, values)
!$OMP n_singles_b, nkeys, keys, values, big_array_local)
! Alpha/Beta double excitations
! =============================
@ -184,6 +187,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
allocate( big_array_local(N_states,dim1, dim1, dim1, dim1) )
big_array_local(:,:,:,:,:) = 0.d0
kcol_prev=-1
@ -191,8 +196,9 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(dynamic,64)
!$OMP DO SCHEDULE(dynamic)
do k_a=istart+ishift,iend,istep
!print *, 'aa', k_a, '/', iend
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
@ -254,33 +260,36 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
do l= 1, N_states
c_1(l) = u_t(l,l_a) * u_t(l,k_a)
enddo
if(alpha_beta)then
! only ONE contribution
if (nkeys+1 .ge. sze_buff) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
nkeys = 0
endif
else if (spin_trace)then
! TWO contributions
! if(alpha_beta)then
! ! only ONE contribution
! if (nkeys+1 .ge. sze_buff) then
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! nkeys = 0
! endif
! else if (spin_trace)then
! ! TWO contributions
if (nkeys+2 .ge. sze_buff) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
endif
! endif
call orb_range_off_diag_double_to_all_states_ab_dm_buffer(tmp_det,tmp_det2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
enddo
endif
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
enddo
enddo
!$OMP END DO NOWAIT
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,64)
!$OMP DO SCHEDULE(dynamic)
do k_a=istart+ishift,iend,istep
!print *, 'ab', k_a, '/', iend
! Single and double alpha exitations
@ -331,6 +340,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if(alpha_beta.or.spin_trace.or.alpha_alpha)then
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
@ -342,25 +352,27 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
do l= 1, N_states
c_1(l) = u_t(l,l_a) * u_t(l,k_a)
enddo
if(alpha_beta.or.spin_trace.or.alpha_alpha)then
! increment the alpha/beta part for single excitations
if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
! increment the alpha/alpha part for single excitations
if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
endif
enddo
endif
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
nkeys = 0
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! nkeys = 0
! Compute Hij for all alpha doubles
! ----------------------------------
@ -377,14 +389,15 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
c_1(l) = u_t(l,l_a) * u_t(l,k_a)
enddo
if (nkeys+4 .ge. sze_buff) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
call orb_range_off_diag_double_to_all_states_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
enddo
endif
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
nkeys = 0
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! nkeys = 0
! Single and double beta excitations
@ -432,6 +445,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
if(alpha_beta.or.spin_trace.or.beta_beta)then
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
@ -444,23 +458,26 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
do l= 1, N_states
c_1(l) = u_t(l,l_a) * u_t(l,k_a)
enddo
if(alpha_beta.or.spin_trace.or.beta_beta)then
! increment the alpha/beta part for single excitations
if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
! increment the beta /beta part for single excitations
if (nkeys+4 * elec_alpha_num .ge. sze_buff) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
call orb_range_off_diag_single_to_all_states_bb_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
endif
enddo
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
nkeys = 0
endif
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! nkeys = 0
! Compute Hij for all beta doubles
! ----------------------------------
@ -478,7 +495,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
c_1(l) = u_t(l,l_a) * u_t(l,k_a)
enddo
if (nkeys+4 .ge. sze_buff) then
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
call orb_range_off_diag_double_to_all_states_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
@ -487,8 +505,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
enddo
endif
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
nkeys = 0
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! nkeys = 0
! Diagonal contribution
@ -514,15 +532,27 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
c_1(l) = u_t(l,k_a) * u_t(l,k_a)
enddo
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
if (nkeys+elec_alpha_num*elec_alpha_num .ge. sze_buff) then
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
endif
call orb_range_diag_to_all_states_2_rdm_dm_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
nkeys = 0
end do
!$OMP END DO
!$OMP END DO NOWAIT
deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
!$OMP CRITICAL
do i=1,N_states
big_array(:,:,:,:,i) = big_array(:,:,:,:,i) + big_array_local(i,:,:,:,:)
enddo
!$OMP END CRITICAL
deallocate(big_array_local)
!$OMP END PARALLEL
end
@ -550,22 +580,66 @@ subroutine update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,loc
integer :: istate
integer :: i,h1,h2,p1,p2
integer, allocatable :: iorder(:)
integer*8, allocatable :: to_sort(:)
allocate(iorder(nkeys))
do i=1,nkeys
iorder(i) = i
enddo
! If the lock is already taken, sort the keys while waiting for a faster access
if (.not.omp_test_lock(lock_2rdm)) then
allocate(to_sort(nkeys))
do i=1,nkeys
h1 = keys(1,iorder(i))
h2 = keys(2,iorder(i))-1
p1 = keys(3,iorder(i))-1
p2 = keys(4,iorder(i))-1
to_sort(i) = int(h1,8) + int(dim1,8)*(int(h2,8) + int(dim1,8)*(int(p1,8) + int(dim1,8)*int(p2,8)))
enddo
call i8sort(to_sort, iorder, nkeys)
deallocate(to_sort)
call omp_set_lock(lock_2rdm)
endif
! print*,'*************'
! print*,'updating'
! print*,'nkeys',nkeys
do istate = 1, N_st
do i = 1, nkeys
h1 = keys(1,iorder(i))
h2 = keys(2,iorder(i))
p1 = keys(3,iorder(i))
p2 = keys(4,iorder(i))
big_array(h1,h2,p1,p2,istate) = big_array(h1,h2,p1,p2,istate) + values(istate,iorder(i))
enddo
enddo
call omp_unset_lock(lock_2rdm)
deallocate(iorder)
end
subroutine update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
use omp_lib
implicit none
integer, intent(in) :: n_st,nkeys,dim1
integer, intent(in) :: keys(4,nkeys)
double precision, intent(in) :: values(n_st,nkeys)
double precision, intent(inout) :: big_array_local(n_st,dim1,dim1,dim1,dim1)
integer :: istate
integer :: i,h1,h2,p1,p2
do i = 1, nkeys
h1 = keys(1,i)
h2 = keys(2,i)
p1 = keys(3,i)
p2 = keys(4,i)
do istate = 1, N_st
! print*,h1,h2,p1,p2,values(istate,i)
big_array(h1,h2,p1,p2,istate) += values(istate,i)
big_array_local(istate,h1,h2,p1,p2) = big_array_local(istate,h1,h2,p1,p2) + values(istate,i)
enddo
enddo
call omp_unset_lock(lock_2rdm)
end

653
src/two_rdm_routines/update_rdm.irp.f
View File

@ -1,4 +1,15 @@
subroutine orb_range_diag_to_all_states_2_rdm_dm_buffer(det_1,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
logical function is_integer_in_string_local(orb,bitmask,Nint)
use bitmasks
implicit none
integer, intent(in) :: orb, Nint
integer(bit_kind), intent(in) :: bitmask(Nint)
integer :: j, k
k = ishft(orb-1,-bit_kind_shift)+1
j = orb-ishft(k-1,bit_kind_shift)-1
is_integer_in_string_local = iand(bitmask(k), ibset(0_bit_kind, j)) /= 0_bit_kind
end
subroutine orb_range_diag_to_all_states_2_rdm_dm_buffer(det_1,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
use bitmasks
BEGIN_DOC
! routine that update the DIAGONAL PART of the two body rdms in a specific range of orbitals for a given determinant det_1
@ -22,7 +33,7 @@
double precision, intent(in) :: c_1(N_st)
double precision, intent(out) :: values(N_st,sze_buff)
integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys
integer , intent(inout) :: nkeys
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
@ -47,15 +58,16 @@
else if(ispin == 4)then
spin_trace = .True.
endif
call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int)
logical :: is_integer_in_string
logical :: is_integer_in_string_local
integer :: i1,i2,istate
if(alpha_beta)then
do i = 1, n_occ_ab(1)
i1 = occ(i,1)
h1 = list_orb_reverse(i1)
do j = 1, n_occ_ab(2)
i2 = occ(j,2)
h1 = list_orb_reverse(i1)
h2 = list_orb_reverse(i2)
! If alpha/beta, electron 1 is alpha, electron 2 is beta
! Therefore you don't necessayr have symmetry between electron 1 and 2
@ -80,11 +92,12 @@
enddo
else if (alpha_alpha)then
do i = 1, n_occ_ab(1)
i1 = occ(i,1)
h1 = list_orb_reverse(i1)
do j = 1, n_occ_ab(1)
i2 = occ(j,1)
h1 = list_orb_reverse(i1)
h2 = list_orb_reverse(i2)
nkeys += 1
do istate = 1, N_st
@ -104,12 +117,14 @@
keys(4,nkeys) = h1
enddo
enddo
else if (beta_beta)then
do i = 1, n_occ_ab(2)
i1 = occ(i,2)
h1 = list_orb_reverse(i1)
do j = 1, n_occ_ab(2)
i2 = occ(j,2)
h1 = list_orb_reverse(i1)
h2 = list_orb_reverse(i2)
nkeys += 1
do istate = 1, N_st
@ -129,13 +144,15 @@
keys(4,nkeys) = h1
enddo
enddo
else if(spin_trace)then
! 0.5 * (alpha beta + beta alpha)
do i = 1, n_occ_ab(1)
i1 = occ(i,1)
h1 = list_orb_reverse(i1)
do j = 1, n_occ_ab(2)
i2 = occ(j,2)
h1 = list_orb_reverse(i1)
h2 = list_orb_reverse(i2)
nkeys += 1
do istate = 1, N_st
@ -154,12 +171,8 @@
keys(3,nkeys) = h2
keys(4,nkeys) = h1
enddo
enddo
do i = 1, n_occ_ab(1)
i1 = occ(i,1)
do j = 1, n_occ_ab(1)
i2 = occ(j,1)
h1 = list_orb_reverse(i1)
h2 = list_orb_reverse(i2)
nkeys += 1
do istate = 1, N_st
@ -181,9 +194,9 @@
enddo
do i = 1, n_occ_ab(2)
i1 = occ(i,2)
h1 = list_orb_reverse(i1)
do j = 1, n_occ_ab(2)
i2 = occ(j,2)
h1 = list_orb_reverse(i1)
h2 = list_orb_reverse(i2)
nkeys += 1
do istate = 1, N_st
@ -204,28 +217,28 @@
enddo
enddo
endif
end
end
subroutine orb_range_off_diag_double_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
subroutine orb_range_off_diag_double_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
!
! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another
!
! c_1 is the array of the contributions to the rdm for all states
!
! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
!
! ispin determines which spin-spin component of the two-rdm you will update
!
! ispin == 1 :: alpha/ alpha
! ispin == 2 :: beta / beta
! ispin == 3 :: alpha/ beta
! ispin == 4 :: spin traced <=> total two-rdm
!
! here, only ispin == 3 or 4 will do something
! routine that updates the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
!
! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another
!
! c_1 is the array of the contributions to the rdm for all states
!
! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
!
! ispin determines which spin-spin component of the two-rdm you will update
!
! ispin == 1 :: alpha/ alpha
! ispin == 2 :: beta / beta
! ispin == 3 :: alpha/ beta
! ispin == 4 :: spin traced <=> total two-rdm
!
! here, only ispin == 3 or 4 will do something
END_DOC
implicit none
integer, intent(in) :: ispin,sze_buff,N_st
@ -234,26 +247,29 @@
double precision, intent(in) :: c_1(N_st)
double precision, intent(out) :: values(N_st,sze_buff)
integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys
integer , intent(inout) :: nkeys
integer :: i,j,h1,h2,p1,p2,istate
integer :: exc(0:2,2,2)
double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
logical :: is_integer_in_string
alpha_alpha = .False.
beta_beta = .False.
alpha_beta = .False.
spin_trace = .False.
if( ispin == 1)then
alpha_alpha = .True.
else if(ispin == 2)then
beta_beta = .True.
else if(ispin == 3)then
alpha_beta = .True.
else if(ispin == 4)then
spin_trace = .True.
endif
logical :: is_integer_in_string_local
if (ispin <= 2) return
! alpha_alpha = .False.
! beta_beta = .False.
! alpha_beta = .False.
! spin_trace = .False.
! if( ispin == 1)then
! alpha_alpha = .True.
! else if(ispin == 2)then
! beta_beta = .True.
! else if(ispin == 3)then
! alpha_beta = .True.
! else if(ispin == 4)then
! spin_trace = .True.
! endif
call get_double_excitation(det_1,det_2,exc,phase,N_int)
h1 = exc(1,1,1)
if(list_orb_reverse(h1).lt.0)return
h1 = list_orb_reverse(h1)
@ -266,10 +282,11 @@
p2 = exc(1,2,2)
if(list_orb_reverse(p2).lt.0)return
p2 = list_orb_reverse(p2)
if(alpha_beta)then
! if(alpha_beta)then
nkeys += 1
phase = phase * 0.5d0
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -277,36 +294,36 @@
keys(4,nkeys) = p2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p2
keys(4,nkeys) = p1
else if(spin_trace)then
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = p2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p2
keys(4,nkeys) = p1
endif
end
! else if(spin_trace)then
! nkeys += 1
! do istate = 1, N_st
! values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
! enddo
! keys(1,nkeys) = h1
! keys(2,nkeys) = h2
! keys(3,nkeys) = p1
! keys(4,nkeys) = p2
! nkeys += 1
! do istate = 1, N_st
! values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
! enddo
! keys(1,nkeys) = h2
! keys(2,nkeys) = h1
! keys(3,nkeys) = p2
! keys(4,nkeys) = p1
! endif
end
subroutine orb_range_off_diag_single_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
subroutine orb_range_off_diag_single_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
! routine that updates the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
!
! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another
!
@ -333,7 +350,7 @@
double precision, intent(in) :: c_1(N_st)
double precision, intent(out) :: values(N_st,sze_buff)
integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys
integer , intent(inout) :: nkeys
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
@ -342,16 +359,12 @@
double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
logical :: is_integer_in_string
alpha_alpha = .False.
beta_beta = .False.
logical :: is_integer_in_string_local
if (ispin <= 2) return
alpha_beta = .False.
spin_trace = .False.
if( ispin == 1)then
alpha_alpha = .True.
else if(ispin == 2)then
beta_beta = .True.
else if(ispin == 3)then
if(ispin == 3)then
alpha_beta = .True.
else if(ispin == 4)then
spin_trace = .True.
@ -360,129 +373,145 @@
call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
call get_single_excitation(det_1,det_2,exc,phase,N_int)
if(alpha_beta)then
if (exc(0,1,1) == 1) then
! Mono alpha
h1 = exc(1,1,1)
if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = exc(1,2,1)
if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
p1 = list_orb_reverse(p1)
do i = 1, n_occ_ab(2)
h2 = occ(i,2)
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
else
! Mono beta
h1 = exc(1,1,2)
if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = exc(1,2,2)
if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
p1 = list_orb_reverse(p1)
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
endif
else if(spin_trace)then
if (exc(0,1,1) == 1) then
! Mono alpha
h1 = exc(1,1,1)
if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = exc(1,2,1)
if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
p1 = list_orb_reverse(p1)
do i = 1, n_occ_ab(2)
h2 = occ(i,2)
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
else
! Mono beta
h1 = exc(1,1,2)
if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = exc(1,2,2)
if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
p1 = list_orb_reverse(p1)
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
endif
endif
end
subroutine orb_range_off_diag_single_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
if (exc(0,1,1) == 1) then
! Mono alpha
h1 = exc(1,1,1)
if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return
p1 = exc(1,2,1)
if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = list_orb_reverse(p1)
phase = 0.5d0 * phase
do i = 1, n_occ_ab(2)
h2 = occ(i,2)
if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
else
! Mono beta
h1 = exc(1,1,2)
if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return
p1 = exc(1,2,2)
if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = list_orb_reverse(p1)
phase = 0.5d0 * phase
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
endif
else if(spin_trace)then
if (exc(0,1,1) == 1) then
! Mono alpha
h1 = exc(1,1,1)
if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return
p1 = exc(1,2,1)
if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = list_orb_reverse(p1)
phase = 0.5d0 * phase
do i = 1, n_occ_ab(2)
h2 = occ(i,2)
if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
else
! Mono beta
h1 = exc(1,1,2)
if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return
p1 = exc(1,2,2)
if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
p1 = list_orb_reverse(p1)
phase = 0.5d0 * phase
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
endif
endif
end
subroutine orb_range_off_diag_single_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
!
@ -512,7 +541,7 @@
double precision, intent(in) :: c_1(N_st)
double precision, intent(out) :: values(N_st,sze_buff)
integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys
integer , intent(inout) :: nkeys
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
@ -521,40 +550,42 @@
double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
logical :: is_integer_in_string
logical :: is_integer_in_string_local
alpha_alpha = .False.
beta_beta = .False.
alpha_beta = .False.
! beta_beta = .False.
! alpha_beta = .False.
spin_trace = .False.
if( ispin == 1)then
alpha_alpha = .True.
else if(ispin == 2)then
beta_beta = .True.
else if(ispin == 3)then
alpha_beta = .True.
else if(ispin == 4)then
spin_trace = .True.
else
return
endif
call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
! if(alpha_alpha.or.spin_trace)then
call get_single_excitation(det_1,det_2,exc,phase,N_int)
if(alpha_alpha.or.spin_trace)then
if (exc(0,1,1) == 1) then
! Mono alpha
h1 = exc(1,1,1)
if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return
p1 = exc(1,2,1)
if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return
call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
h1 = list_orb_reverse(h1)
p1 = list_orb_reverse(p1)
phase = 0.5d0*phase
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -563,7 +594,16 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -572,29 +612,20 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = h2
enddo
else
return
endif
endif
end
! endif
end
subroutine orb_range_off_diag_single_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
subroutine orb_range_off_diag_single_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
@ -624,7 +655,7 @@
double precision, intent(in) :: c_1(N_st)
double precision, intent(out) :: values(N_st,sze_buff)
integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys
integer , intent(inout) :: nkeys
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
@ -632,42 +663,43 @@
integer :: exc(0:2,2,2)
double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
logical :: is_integer_in_string
alpha_alpha = .False.
logical :: is_integer_in_string_local
! alpha_alpha = .False.
beta_beta = .False.
alpha_beta = .False.
! alpha_beta = .False.
spin_trace = .False.
if( ispin == 1)then
alpha_alpha = .True.
else if(ispin == 2)then
if(ispin == 2)then
beta_beta = .True.
else if(ispin == 3)then
alpha_beta = .True.
else if(ispin == 4)then
spin_trace = .True.
else
return
endif
call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
! if(beta_beta.or.spin_trace)then
call get_single_excitation(det_1,det_2,exc,phase,N_int)
if(beta_beta.or.spin_trace)then
if (exc(0,1,1) == 1) then
return
else
! Mono beta
h1 = exc(1,1,2)
if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return
h1 = list_orb_reverse(h1)
if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return
p1 = exc(1,2,2)
if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return
call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
h1 = list_orb_reverse(h1)
p1 = list_orb_reverse(p1)
phase = 0.5d0*phase
do i = 1, n_occ_ab(2)
h2 = occ(i,2)
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2)
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -676,7 +708,16 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -685,28 +726,19 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = h2
enddo
endif
endif
end
! endif
end
subroutine orb_range_off_diag_double_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
subroutine orb_range_off_diag_double_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
@ -735,7 +767,7 @@
double precision, intent(in) :: c_1(N_st)
double precision, intent(out) :: values(N_st,sze_buff)
integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys
integer , intent(inout) :: nkeys
integer :: i,j,h1,h2,p1,p2,istate
@ -743,38 +775,39 @@
double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
logical :: is_integer_in_string
logical :: is_integer_in_string_local
alpha_alpha = .False.
beta_beta = .False.
alpha_beta = .False.
! beta_beta = .False.
! alpha_beta = .False.
spin_trace = .False.
if( ispin == 1)then
alpha_alpha = .True.
else if(ispin == 2)then
beta_beta = .True.
else if(ispin == 3)then
alpha_beta = .True.
else if(ispin == 4)then
spin_trace = .True.
else
return
endif
call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
h1 =exc(1,1)
if(list_orb_reverse(h1).lt.0)return
h1 = list_orb_reverse(h1)
h2 =exc(2,1)
if(list_orb_reverse(h2).lt.0)return
h2 = list_orb_reverse(h2)
p1 =exc(1,2)
if(list_orb_reverse(p1).lt.0)return
p1 = list_orb_reverse(p1)
p2 =exc(2,2)
if(list_orb_reverse(p2).lt.0)return
h1 = list_orb_reverse(h1)
h2 = list_orb_reverse(h2)
p1 = list_orb_reverse(p1)
p2 = list_orb_reverse(p2)
if(alpha_alpha.or.spin_trace)then
phase = 0.5d0*phase
! if(alpha_alpha.or.spin_trace)then
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -783,7 +816,16 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = p2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -792,25 +834,16 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p2
keys(4,nkeys) = p1
! endif
end
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = p2
endif
end
subroutine orb_range_off_diag_double_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
subroutine orb_range_off_diag_double_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
@ -840,25 +873,23 @@
double precision, intent(in) :: c_1(N_st)
double precision, intent(out) :: values(N_st,sze_buff)
integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys
integer , intent(inout) :: nkeys
integer :: i,j,h1,h2,p1,p2,istate
integer :: exc(0:2,2)
double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
logical :: is_integer_in_string
alpha_alpha = .False.
logical :: is_integer_in_string_local
! alpha_alpha = .False.
beta_beta = .False.
alpha_beta = .False.
! alpha_beta = .False.
spin_trace = .False.
if( ispin == 1)then
alpha_alpha = .True.
else if(ispin == 2)then
if(ispin == 2)then
beta_beta = .True.
else if(ispin == 3)then
alpha_beta = .True.
else if(ispin == 4)then
spin_trace = .True.
else
return
endif
call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
@ -874,10 +905,12 @@
p2 =exc(2,2)
if(list_orb_reverse(p2).lt.0)return
p2 = list_orb_reverse(p2)
if(beta_beta.or.spin_trace)then
! if(beta_beta.or.spin_trace)then
phase = 0.5d0*phase
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -886,7 +919,16 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = p2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - c_1(istate) * phase
enddo
keys(1,nkeys) = h1
keys(2,nkeys) = h2
@ -895,21 +937,12 @@
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
values(istate,nkeys) = c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p2
keys(4,nkeys) = p1
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p1
keys(4,nkeys) = p2
endif
end
! endif
end