mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2025-01-03 18:15:58 +01:00
Plugin for reading AO integrals
This commit is contained in:
parent
ce2f6f4a67
commit
5e60f68b89
5
devel/density_matrix/.gitignore
vendored
5
devel/density_matrix/.gitignore
vendored
@ -1,5 +0,0 @@
|
||||
IRPF90_temp/
|
||||
IRPF90_man/
|
||||
irpf90.make
|
||||
irpf90_entities
|
||||
tags
|
||||
@ -1 +0,0 @@
|
||||
determinants
|
||||
@ -1,5 +0,0 @@
|
||||
====================
|
||||
DensityMatrix Module
|
||||
====================
|
||||
|
||||
Prints the 1- and 2- body density matrices
|
||||
@ -1,115 +0,0 @@
|
||||
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
|
||||
|
||||
|
||||
@ -1,43 +0,0 @@
|
||||
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
|
||||
@ -1,78 +0,0 @@
|
||||
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
|
||||
59
devel/import_integrals_ao/.gitignore
vendored
Normal file
59
devel/import_integrals_ao/.gitignore
vendored
Normal 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
|
||||
1
devel/import_integrals_ao/NEED
Normal file
1
devel/import_integrals_ao/NEED
Normal file
@ -0,0 +1 @@
|
||||
nuclei ao_one_e_ints ao_two_e_ints
|
||||
28
devel/import_integrals_ao/README.rst
Normal file
28
devel/import_integrals_ao/README.rst
Normal file
@ -0,0 +1,28 @@
|
||||
===================
|
||||
import_integrals_ao
|
||||
===================
|
||||
|
||||
Module to read all the integrals in the |AO| basis from files (all in atomic units).
|
||||
|
||||
|
||||
The following files are required:
|
||||
|
||||
- :file:`S.qp` : overlap integrals
|
||||
- :file:`T.qp` : kinetic integrals
|
||||
- :file:`V.qp` : electron-nucleus potential integrals
|
||||
- :file:`P.qp` : pseudo-potential integrals
|
||||
- :file:`W.qp` : electron repulsion integrals
|
||||
|
||||
If present, the :file:`E.qp` file, should contain the nuclear repulsion energy.
|
||||
|
||||
In all the other files, there is one integral per line and for the one-electron integral
|
||||
$\int \chi_i(r) \hat{O} \chi_j(r) dr$, the format is
|
||||
|
||||
i j value
|
||||
|
||||
and for two electron integral the format uses the physicists' convention,
|
||||
$\int \chi_i(r_1) \chi_j(r_2) \hat{O} \chi_k(r_1) \chi_l(r_2) dr$:
|
||||
|
||||
i j k l value
|
||||
|
||||
|
||||
108
devel/import_integrals_ao/import_integrals_ao.irp.f
Normal file
108
devel/import_integrals_ao/import_integrals_ao.irp.f
Normal file
@ -0,0 +1,108 @@
|
||||
program print_integrals
|
||||
print *, 'Number of AOs?'
|
||||
read(*,*) ao_num
|
||||
TOUCH ao_num
|
||||
call run
|
||||
end
|
||||
|
||||
subroutine run
|
||||
use map_module
|
||||
implicit none
|
||||
|
||||
integer :: iunit
|
||||
integer :: getunitandopen
|
||||
|
||||
integer ::i,j,k,l
|
||||
double precision :: integral
|
||||
double precision, allocatable :: A(:,:)
|
||||
|
||||
integer :: n_integrals
|
||||
integer(key_kind), allocatable :: buffer_i(:)
|
||||
real(integral_kind), allocatable :: buffer_values(:)
|
||||
|
||||
call ezfio_set_ao_basis_ao_num(ao_num)
|
||||
|
||||
allocate (A(ao_num,ao_num))
|
||||
|
||||
A(1,1) = huge(1.d0)
|
||||
iunit = getunitandopen('E.qp','r')
|
||||
read (iunit,*,end=9) A(1,1)
|
||||
9 continue
|
||||
close(iunit)
|
||||
if (A(1,1) /= huge(1.d0)) then
|
||||
call ezfio_set_nuclei_nuclear_repulsion(A(1,1))
|
||||
call ezfio_set_nuclei_io_nuclear_repulsion("Read")
|
||||
endif
|
||||
|
||||
A = 0.d0
|
||||
iunit = getunitandopen('T.qp','r')
|
||||
do
|
||||
read (iunit,*,end=10) i,j, integral
|
||||
A(i,j) = integral
|
||||
enddo
|
||||
10 continue
|
||||
close(iunit)
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A(1:ao_num, 1:ao_num))
|
||||
call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic("Read")
|
||||
|
||||
A = 0.d0
|
||||
iunit = getunitandopen('S.qp','r')
|
||||
do
|
||||
read (iunit,*,end=11) i,j, integral
|
||||
A(i,j) = integral
|
||||
enddo
|
||||
11 continue
|
||||
close(iunit)
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A(1:ao_num, 1:ao_num))
|
||||
call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap("Read")
|
||||
|
||||
A = 0.d0
|
||||
iunit = getunitandopen('P.qp','r')
|
||||
do
|
||||
read (iunit,*,end=14) i,j, integral
|
||||
A(i,j) = integral
|
||||
enddo
|
||||
14 continue
|
||||
close(iunit)
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A(1:ao_num,1:ao_num))
|
||||
call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo("Read")
|
||||
|
||||
A = 0.d0
|
||||
iunit = getunitandopen('V.qp','r')
|
||||
do
|
||||
read (iunit,*,end=12) i,j, integral
|
||||
A(i,j) = integral
|
||||
enddo
|
||||
12 continue
|
||||
close(iunit)
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_e_n(A(1:ao_num, 1:ao_num))
|
||||
call ezfio_set_ao_one_e_ints_io_ao_integrals_e_n("Read")
|
||||
|
||||
allocate(buffer_i(ao_num**3), buffer_values(ao_num**3))
|
||||
iunit = getunitandopen('W.qp','r')
|
||||
n_integrals=0
|
||||
buffer_values = 0.d0
|
||||
do
|
||||
read (iunit,*,end=13) i,j,k,l, integral
|
||||
n_integrals += 1
|
||||
call two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) )
|
||||
buffer_values(n_integrals) = integral
|
||||
if (n_integrals == size(buffer_i)) then
|
||||
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
|
||||
n_integrals = 0
|
||||
endif
|
||||
enddo
|
||||
13 continue
|
||||
close(iunit)
|
||||
|
||||
if (n_integrals > 0) then
|
||||
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
|
||||
endif
|
||||
|
||||
call map_sort(ao_integrals_map)
|
||||
call map_unique(ao_integrals_map)
|
||||
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
|
||||
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
|
||||
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user