scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2025-04-10 18:51:28 +02:00
Code Issues Releases Wiki Activity

Added 2RDM plugin

Browse Source
This commit is contained in:
Anthony Scemama 2019-05-15 09:43:03 +02:00
parent 9f58127726
commit 39250742fd
8 changed files with 311 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
devel/density_matrix/.gitignore vendored Normal file
View File

@ -0,0 +1,5 @@
IRPF90_temp/
IRPF90_man/
irpf90.make
irpf90_entities
tags

1
devel/density_matrix/NEED Normal file
View File

@ -0,0 +1 @@
determinants

5
devel/density_matrix/README.rst Normal file
View File

@ -0,0 +1,5 @@
====================
DensityMatrix Module
====================
Prints the 1- and 2- body density matrices

115
devel/density_matrix/density_matrix_array.irp.f Normal file
View File

@ -0,0 +1,115 @@
use bitmasks
BEGIN_PROVIDER [ double precision, two_e_dm_aa, (mo_num,mo_num,mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, two_e_dm_bb, (mo_num,mo_num,mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, two_e_dm_ab, (mo_num,mo_num,mo_num,mo_num) ]
implicit none
use bitmasks
BEGIN_DOC
! Temporary files for 2-e dm calculation
END_DOC
integer :: getUnitAndOpen
integer :: k,l,degree, idx,i,j
integer :: exc(0:2,2,2),n_occ_alpha
double precision :: phase, coef
integer :: h1,h2,p1,p2,s1,s2, e1, e2
double precision :: ck, cl
character*(128), parameter :: f = '(i8,4(x,i5),x,d16.8)'
integer :: istate
two_e_dm_aa = 0.d0
two_e_dm_ab = 0.d0
two_e_dm_bb = 0.d0
istate = 1
! OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,ck,ckl,i,j,e1,e2,cl,phase,h1,p1,h2,p2,s1,s2,occ)
! OMP DO SCHEDULE(dynamic,64)
do k=1,N_det
ck = psi_coef(k,istate)
call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
ckl = psi_coef(k,istate) * psi_coef(k,istate)
do i = 1,elec_alpha_num
e1=occ(i,1)
do j = 1,elec_alpha_num
e2=occ(j,1)
! alpha-alpha
two_e_dm_aa(e1,e2,e1,e2) += 0.5d0*ckl
two_e_dm_aa(e1,e2,e2,e1) -= 0.5d0*ckl
enddo
do j = 1,elec_beta_num
e2=occ(j,2)
! alpha-beta
two_e_dm_ab(e1,e2,e1,e2) += ckl
enddo
enddo
do i = 1,elec_beta_num
e1=occ(i,2)
do j = 1,elec_beta_num
e2=occ(j,2)
! beta-beta
two_e_dm_bb(e1,e2,e1,e2) += 0.5d0*ckl
two_e_dm_bb(e1,e2,e2,e1) -= 0.5d0*ckl
enddo
enddo
do l=1,k-1
cl = 2.d0*psi_coef(l,istate)
call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
if (degree == 2) then
call get_double_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
ckl = phase*ck*cl
select case (s1+s2)
case(2) ! alpha alpha
two_e_dm_aa(h1,h2,p1,p2) += ckl
two_e_dm_aa(h1,h2,p2,p1) -= ckl
case(3) ! alpha beta
two_e_dm_ab(h1,h2,p1,p2) += ckl
case(4) ! beta beta
two_e_dm_bb(h1,h2,p1,p2) += ckl
two_e_dm_bb(h1,h2,p2,p1) -= ckl
end select
else if (degree == 1) then
call get_single_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
double precision :: ckl
ckl = phase*ck*cl
call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
select case (s1)
case (1) ! Alpha single excitation
integer :: occ(N_int*bit_kind_size,2)
do i = 1, elec_alpha_num
p2=occ(i,1)
h2=p2
two_e_dm_aa(h1,h2,p1,p2) += ckl
two_e_dm_aa(h1,h2,p2,p1) -= ckl
enddo
do i = 1, elec_beta_num
p2=occ(i,2)
h2=p2
two_e_dm_ab(h1,h2,p1,p2) += ckl
enddo
case (2) ! Beta single excitation
do i = 1, elec_alpha_num
p2=occ(i,1)
h2=p2
two_e_dm_ab(h1,h2,p1,p2) += ckl
enddo
do i = 1, elec_beta_num
p2=occ(i,2)
h2=p2
two_e_dm_bb(h1,h2,p1,p2) += ckl
two_e_dm_bb(h1,h2,p2,p1) -= ckl
enddo
end select
endif
enddo
enddo
! OMP END DO
! OMP END PARALLEL
END_PROVIDER

43
devel/density_matrix/print_2rdm.irp.f Normal file
View File

@ -0,0 +1,43 @@
program print_2rdm
implicit none
integer :: i,j,k,l
double precision, external :: get_two_e_integral
read_wf = .True.
TOUCH read_wf
double precision :: e(10)
double precision, parameter :: thr = 1.d-15
e = 0.d0
print *, '1RDM'
do i=1,mo_num
do j=1,mo_num
if (dabs(one_e_dm_mo_alpha(i,j,1) + one_e_dm_mo_beta(i,j,1)) > thr) then
print *, i, j, one_e_dm_mo_alpha(i,j,1) + one_e_dm_mo_beta(i,j,1)
endif
e(4) += one_e_dm_mo_alpha(i,j,1) * mo_one_e_integrals(i,j)
e(4) += one_e_dm_mo_beta(i,j,1) * mo_one_e_integrals(i,j)
enddo
enddo
print *, '2RDM'
do i=1,mo_num
do j=1,mo_num
do k=1,mo_num
do l=1,mo_num
if (dabs(two_e_dm_aa(i,j,k,l) + two_e_dm_bb(i,j,k,l) + two_e_dm_ab(i,j,k,l)) > thr) then
print *, i, j, k, l, two_e_dm_aa(i,j,k,l) + two_e_dm_bb(i,j,k,l) + two_e_dm_ab(i,j,k,l)
endif
e(1) += two_e_dm_aa(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
e(2) += two_e_dm_bb(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
e(3) += two_e_dm_ab(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, ''
print *, 'Energy ', sum(e(1:4)) + nuclear_repulsion
end

78
devel/density_matrix/print_2rdm_decomposed.irp.f Normal file
View File

@ -0,0 +1,78 @@
program print_2rdm_decomposed
implicit none
integer :: i,j,k,l
double precision, external :: get_two_e_integral
read_wf = .True.
TOUCH read_wf
double precision :: e(10)
double precision, parameter :: thr = 1.d-15
e = 0.d0
print *, '1RDM ALPHA'
do i=1,mo_num
do j=1,mo_num
if (dabs(one_e_dm_mo_alpha(i,j,1)) > thr) then
print *, i, j, one_e_dm_mo_alpha(i,j,1)
endif
e(4) += one_e_dm_mo_alpha(i,j,1) * mo_one_e_integrals(i,j)
enddo
enddo
print *, '1RDM BETA'
do i=1,mo_num
do j=1,mo_num
if (dabs(one_e_dm_mo_beta(i,j,1)) > thr) then
print *, i, j, one_e_dm_mo_beta(i,j,1)
endif
e(4) += one_e_dm_mo_beta(i,j,1) * mo_one_e_integrals(i,j)
enddo
enddo
print *, '2RDM ALPHA ALPHA'
do i=1,mo_num
do j=1,mo_num
do k=1,mo_num
do l=1,mo_num
if (dabs(two_e_dm_aa(i,j,k,l)) > thr) then
print *, i, j, k, l, two_e_dm_aa(i,j,k,l)
endif
e(1) += two_e_dm_aa(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, '2RDM BETA BETA'
do i=1,mo_num
do j=1,mo_num
do k=1,mo_num
do l=1,mo_num
if (dabs(two_e_dm_bb(i,j,k,l)) > thr) then
print *, i, j, k, l, two_e_dm_bb(i,j,k,l)
endif
e(2) += two_e_dm_bb(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, '2RDM ALPHA BETA'
do i=1,mo_num
do j=1,mo_num
do k=1,mo_num
do l=1,mo_num
if (dabs(two_e_dm_ab(i,j,k,l)) > thr) then
print *, i, j, k, l, two_e_dm_ab(i,j,k,l)
endif
e(3) += two_e_dm_ab(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, ''
print *, 'Energy ', sum(e(1:4)) + nuclear_repulsion
end

5
stable/champ/.gitignore vendored Normal file
View File

@ -0,0 +1,5 @@
IRPF90_temp/
IRPF90_man/
irpf90.make
irpf90_entities
tags

59
stable/mp2/.gitignore vendored Normal file
View File

@ -0,0 +1,59 @@
IRPF90_temp/
IRPF90_man/
build.ninja
irpf90.make
ezfio_interface.irp.f
irpf90_entities
tags
Makefile
ao_basis
ao_one_e_ints
ao_two_e_erf_ints
ao_two_e_ints
aux_quantities
becke_numerical_grid
bitmask
cis
cisd
cipsi
davidson
davidson_dressed
davidson_undressed
density_for_dft
determinants
dft_keywords
dft_utils_in_r
dft_utils_one_e
dft_utils_two_body
dressing
dummy
electrons
ezfio_files
fci
generators_cas
generators_full
hartree_fock
iterations
kohn_sham
kohn_sham_rs
mo_basis
mo_guess
mo_one_e_ints
mo_two_e_erf_ints
mo_two_e_ints
mpi
mrpt_utils
nuclei
perturbation
pseudo
psiref_cas
psiref_utils
scf_utils
selectors_cassd
selectors_full
selectors_utils
single_ref_method
slave
tools
utils
zmq