scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-07-25 20:27:35 +02:00
Code Issues Releases Wiki Activity

Merge branch 'master' of https://gitlab.com/scemama/qp_plugins_scemama

Browse Source
This commit is contained in:
Emmanuel Giner 2021-07-01 11:41:31 +02:00
commit d5192bcbc9
33 changed files with 4800 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
devel/fci_complete/EZFIO.cfg Normal file
View File

@ -0,0 +1,6 @@
[energy]
type: double precision
doc: Calculated Selected |FCI| energy
interface: ezfio
size: (determinants.n_states)

3
devel/fci_complete/NEED Normal file
View File

@ -0,0 +1,3 @@
davidson_undressed
hartree_fock
determinants

4
devel/fci_complete/README.rst Normal file
View File

@ -0,0 +1,4 @@
===
mpn
===

8
devel/fci_complete/fci_complete.irp.f Normal file
View File

@ -0,0 +1,8 @@
program fci_complete
call generate_fci_space
call diagonalize_ci
end

75
devel/fci_complete/generate_fci.irp.f Normal file
View File

@ -0,0 +1,75 @@
subroutine generate_fci_space
use bitmasks
implicit none
BEGIN_DOC
! Generates the complete FCI space
END_DOC
integer :: i, sze, ncore
integer(bit_kind) :: o(N_int,2)
integer(bit_kind) :: u, coremask
if (mo_num > 64) then
stop 'No more than 64 MOs'
endif
ncore = 0
coremask = 0_bit_kind
do i=1,mo_num
if (trim(mo_class(i)) == 'Core') then
ncore += 1
coremask = ibset(coremask,i-1)
endif
enddo
o(1,1) = iand(full_ijkl_bitmask(1),not(coremask))
o(1,2) = 0_bit_kind
call configuration_to_dets_size(o,n_det_alpha_unique,elec_alpha_num-ncore,N_int)
TOUCH n_det_alpha_unique
integer :: k,n,m, t, t1, t2
k=0
n = elec_alpha_num
m = mo_num - n
n = n
u = shiftl(1_bit_kind,n) -1
do while (u < shiftl(1_bit_kind,n+m))
if (iand(coremask, u) == coremask) then
k = k+1
psi_det_alpha_unique(1,k) = u
endif
t = ior(u,u-1)
t1 = t+1
IRP_IF WITHOUT_TRAILZ
t2 = shiftr((iand(not(t),t1)-1), popcnt(ieor(u,u-1)))
IRP_ELSE
t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
IRP_ENDIF
u = ior(t1,t2)
enddo
call configuration_to_dets_size(o,n_det_beta_unique,elec_beta_num-ncore,N_int)
TOUCH n_det_beta_unique
k=0
n = elec_beta_num
m = mo_num - n
u = shiftl(1_bit_kind,n) -1
do while (u < shiftl(1_bit_kind,n+m))
if (iand(coremask, u) == coremask) then
k = k+1
psi_det_beta_unique(1,k) = u
endif
t = ior(u,u-1)
t1 = t+1
t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
u = ior(t1,t2)
enddo
call generate_all_alpha_beta_det_products
print *, 'Ndet = ', N_det
end

1
devel/import_integrals/export_integrals_ao.irp.f
View File

@ -77,6 +77,7 @@ subroutine run
do k=1,ao_num
do j=l,ao_num
do i=k,ao_num
if (i==j .and. k<l) cycle
if (i>=j) then
integral = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
if (integral /= 0.d0) then

12
devel/import_integrals/export_integrals_mo.irp.f
View File

@ -1,4 +1,5 @@
program export_integrals_mo
PROVIDE mo_two_e_integrals_in_map
call run
end
@ -53,7 +54,8 @@ subroutine run
integer :: n_integrals
integer(key_kind) :: idx
double precision, external :: get_two_e_integral
double precision, external :: mo_two_e_integral
allocate (A(mo_num,mo_num))
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
@ -61,13 +63,15 @@ subroutine run
call write_2d('T_mo.qp',mo_kinetic_integrals)
call write_2d('P_mo.qp',mo_pseudo_integrals)
call write_2d('V_mo.qp',mo_integrals_n_e)
call write_2d('H_mo.qp',mo_one_e_integrals)
iunit = getunitandopen('W_mo.qp','w')
do l=1,mo_num
do k=1,mo_num
do j=1,mo_num
do i=1,mo_num
integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
do j=l,mo_num
do i=k,mo_num
if (i==j .and. k<l) cycle
integral = mo_two_e_integral(i,j,k,l)
if (integral /= 0.d0) then
write (iunit,'(4(I5,2X),E22.15)') i,j,k,l, integral
endif

207
devel/import_integrals/import_integrals_mo.irp.f
View File

@ -17,31 +17,36 @@ subroutine run
!
! E.qp : Contains the nuclear repulsion energy
!
! Tmo.qp : kinetic energy integrals
! T_mo.qp : kinetic energy integrals
!
! Smo.qp : overlap matrix
! P_mo.qp : pseudopotential integrals
!
! Pmo.qp : pseudopotential integrals
! V_mo.qp : electron-nucleus potential
!
! Vmo.qp : electron-nucleus potential
! H_mo.qp : core hamiltonian. If hmo exists, the other 1-e files are not read
!
! Wmo.qp : electron repulsion integrals
! W_mo.qp : electron repulsion integrals
!
END_DOC
integer :: iunit
integer :: getunitandopen
integer :: iunit
integer :: getunitandopen
integer ::i,j,k,l
double precision :: integral
double precision, allocatable :: A(:,:)
integer :: i,j,k,l
double precision :: integral
double precision, allocatable :: A(:,:)
integer :: n_integrals
integer(key_kind), allocatable :: buffer_i(:)
integer :: n_integrals
logical :: exists
integer(key_kind), allocatable :: buffer_i(:)
real(integral_kind), allocatable :: buffer_values(:)
allocate(buffer_i(mo_num**3), buffer_values(mo_num**3))
allocate (A(mo_num,mo_num))
PROVIDE mo_integrals_map
PROVIDE mo_integrals_threshold
! Nuclear repulsion
A(1,1) = huge(1.d0)
iunit = getunitandopen('E.qp','r')
@ -53,70 +58,109 @@ subroutine run
call ezfio_set_nuclei_io_nuclear_repulsion('Read')
endif
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('Tmo.qp','r')
do
read (iunit,*,end=10) i,j, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in Tmo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in Tmo.qp'
endif
A(i,j) = integral
enddo
10 continue
close(iunit)
call ezfio_set_mo_one_e_ints_mo_integrals_kinetic(A)
call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('Read')
! One-electron integrals
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('Pmo.qp','r')
do
read (iunit,*,end=14) i,j, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in Pmo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in Pmo.qp'
endif
A(i,j) = integral
enddo
14 continue
close(iunit)
call ezfio_set_mo_one_e_ints_mo_integrals_pseudo(A)
call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('Read')
exists = .False.
inquire(file='H_mo.qp',exist=exists)
if (exists) then
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('Vmo.qp','r')
do
read (iunit,*,end=12) i,j, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in Vmo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in Vmo.qp'
endif
A(i,j) = integral
enddo
12 continue
close(iunit)
call ezfio_set_mo_one_e_ints_mo_integrals_n_e(A)
call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('Read')
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('H_mo.qp','r')
do
read (iunit,*,end=8) i,j, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in H_mo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in H_mo.qp'
endif
A(i,j) = integral
enddo
8 continue
close(iunit)
call ezfio_set_mo_one_e_ints_mo_one_e_integrals(A)
call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('Read')
call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None')
call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None')
call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None')
iunit = getunitandopen('Wmo.qp','r')
else
call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None')
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('T_mo.qp','r')
do
read (iunit,*,end=10) i,j, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in T_mo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in T_mo.qp'
endif
A(i,j) = integral
enddo
10 continue
close(iunit)
call ezfio_set_mo_one_e_ints_mo_integrals_kinetic(A)
call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('Read')
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('P_mo.qp','r')
do
read (iunit,*,end=14) i,j, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in P_mo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in P_mo.qp'
endif
A(i,j) = integral
enddo
14 continue
close(iunit)
call ezfio_set_mo_one_e_ints_mo_integrals_pseudo(A)
call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('Read')
A = 0.d0
i = 1
j = 1
iunit = getunitandopen('V_mo.qp','r')
do
read (iunit,*,end=12) i,j, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in V_mo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in V_mo.qp'
endif
A(i,j) = integral
enddo
12 continue
close(iunit)
call ezfio_set_mo_one_e_ints_mo_integrals_n_e(A)
call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('Read')
end if
! Two-electron integrals
iunit = getunitandopen('W_mo.qp','r')
n_integrals=0
i = 1
j = 1
@ -127,26 +171,25 @@ subroutine run
read (iunit,*,end=13) i,j,k,l, integral
if (i<0 .or. i>mo_num) then
print *, i
stop 'i out of bounds in Wmo.qp'
stop 'i out of bounds in W_mo.qp'
endif
if (j<0 .or. j>mo_num) then
print *, j
stop 'j out of bounds in Wmo.qp'
stop 'j out of bounds in W_mo.qp'
endif
if (k<0 .or. k>mo_num) then
print *, k
stop 'k out of bounds in Wmo.qp'
stop 'k out of bounds in W_mo.qp'
endif
if (l<0 .or. l>mo_num) then
print *, l
stop 'l out of bounds in Wmo.qp'
stop 'l out of bounds in W_mo.qp'
endif
n_integrals += 1
call mo_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_mo_integrals_map(n_integrals, buffer_i, buffer_values, &
real(mo_integrals_threshold,integral_kind))
call map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals)
n_integrals = 0
endif
enddo
@ -154,13 +197,15 @@ subroutine run
close(iunit)
if (n_integrals > 0) then
call insert_into_mo_integrals_map(n_integrals, buffer_i, buffer_values, &
real(mo_integrals_threshold,integral_kind))
call map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals)
endif
call map_sort(mo_integrals_map)
call map_unique(mo_integrals_map)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
end

59
devel/mpn/.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

12
devel/mpn/EZFIO.cfg Normal file
View File

@ -0,0 +1,12 @@
[energy]
type: double precision
doc: Calculated Selected |FCI| energy
interface: ezfio
size: (determinants.n_states)
[mp_order]
type: integer
doc: Max order of MPn
interface: ezfio, provider, ocaml
default: 4

3
devel/mpn/NEED Normal file
View File

@ -0,0 +1,3 @@
davidson_undressed
hartree_fock
determinants

4
devel/mpn/README.rst Normal file
View File

@ -0,0 +1,4 @@
===
mpn
===

23
devel/mpn/energies.irp.f Normal file
View File

@ -0,0 +1,23 @@
BEGIN_PROVIDER [ double precision, energy_det_i, (N_det) ]
implicit none
BEGIN_DOC
! Fock Energy of determinant |I> (sum of epsilon_i)
END_DOC
integer :: i, k, n
integer :: list(elec_alpha_num)
do k=1,N_det
call bitstring_to_list(psi_det(1,1,k), list, n, N_int)
energy_det_i(k) = 0.d0
do i=1,n
energy_det_i(k) += fock_matrix_diag_mo(list(i))
enddo
call bitstring_to_list(psi_det(1,2,k), list, n, N_int)
do i=1,n
energy_det_i(k) += fock_matrix_diag_mo(list(i))
enddo
enddo
END_PROVIDER

80
devel/mpn/generate_fci.irp.f Normal file
View File

@ -0,0 +1,80 @@
subroutine generate_fci_space
use bitmasks
implicit none
BEGIN_DOC
! Generates the complete FCI space
END_DOC
integer :: i, sze, ncore
integer(bit_kind) :: o(N_int,2)
integer(bit_kind) :: u, coremask
if (mo_num > 64) then
stop 'No more than 64 MOs'
endif
ncore = 0
coremask = 0_bit_kind
do i=1,mo_num
if (trim(mo_class(i)) == 'Core') then
ncore += 1
coremask = ibset(coremask,i-1)
endif
enddo
o(1,1) = iand(full_ijkl_bitmask(1),not(coremask))
o(1,2) = 0_bit_kind
call configuration_to_dets_size(o,n_det_alpha_unique,elec_alpha_num-ncore,N_int)
TOUCH n_det_alpha_unique
integer :: k,n,m, t, t1, t2
k=0
n = elec_alpha_num
m = mo_num - n
n = n
u = shiftl(1_bit_kind,n) -1
do while (u < shiftl(1_bit_kind,n+m))
if (iand(coremask, u) == coremask) then
k = k+1
psi_det_alpha_unique(1,k) = u
endif
t = ior(u,u-1)
t1 = t+1
IRP_IF WITHOUT_TRAILZ
t2 = shiftr((iand(not(t),t1)-1), popcnt(ieor(u,u-1)))
IRP_ELSE
t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
IRP_ENDIF
u = ior(t1,t2)
enddo
call configuration_to_dets_size(o,n_det_beta_unique,elec_beta_num-ncore,N_int)
TOUCH n_det_beta_unique
k=0
n = elec_beta_num
m = mo_num - n
u = shiftl(1_bit_kind,n) -1
do while (u < shiftl(1_bit_kind,n+m))
if (iand(coremask, u) == coremask) then
k = k+1
psi_det_beta_unique(1,k) = u
endif
t = ior(u,u-1)
t1 = t+1
IRP_IF WITHOUT_TRAILZ
t2 = shiftr((iand(not(t),t1)-1), popcnt(ieor(u,u-1)))
IRP_ELSE
t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
IRP_ENDIF
u = ior(t1,t2)
enddo
call generate_all_alpha_beta_det_products
print *, N_det
end

53
devel/mpn/mpn.irp.f Normal file
View File

@ -0,0 +1,53 @@
program mpn
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
integer :: i, k, l
double precision, allocatable :: c_pert(:,:)
double precision, allocatable :: e_pert(:)
double precision, allocatable :: hc(:), s2(:)
n_states_diag = 1
TOUCH n_states_diag
call generate_fci_space
allocate(c_pert(N_det,0:mp_order))
allocate(s2(N_det))
allocate(e_pert(0:mp_order))
e_pert(0) = energy_det_i(1)
c_pert(:,:) = 0.d0
c_pert(1,0) = 1.d0
e_pert(1) = hf_energy - e_pert(0) - nuclear_repulsion
do k=1,mp_order
! H_ij C^(k-1)
if (distributed_davidson) then
call H_S2_u_0_nstates_zmq (c_pert(1,k),s2,c_pert(1,k-1),1,N_det)
else
call H_S2_u_0_nstates_openmp(c_pert(1,k),s2,c_pert(1,k-1),1,N_det)
endif
if (k>1) then
e_pert(k) += c_pert(1,k)
endif
print *, k, e_pert(k), sum(e_pert) + nuclear_repulsion
c_pert(:,k) = -c_pert(:,k)
c_pert(1,k) = 0.d0
do l=1,k
do i=2,N_det
c_pert(i,k) = c_pert(i,k) + e_pert(l) * c_pert(i,k-l)
enddo
enddo
do i=2,N_det
c_pert(i,k) = c_pert(i,k) + energy_det_i(i) * c_pert(i,k-1)
enddo
do i=2,N_det
c_pert(i,k) = c_pert(i,k) / (energy_det_i(i) - energy_det_i(1))
enddo
enddo
end

57
devel/qmcchem/qmc_e_curve.irp.f
View File

@ -6,16 +6,13 @@ program e_curve
integer, allocatable :: iorder(:)
double precision , allocatable :: norm_sort(:)
double precision :: e_0(N_states)
PROVIDE mo_two_e_integrals_in_map
PROVIDE mo_two_e_integrals_in_map mo_one_e_integrals
nab = n_det_alpha_unique+n_det_beta_unique
allocate ( norm_sort(0:nab), iorder(0:nab) )
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
double precision, allocatable :: u_0(:,:), v_0(:,:)
allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
allocate(u_0(N_states,N_det),v_0(N_states,N_det))
norm_sort(0) = 0.d0
@ -24,19 +21,20 @@ program e_curve
norm_sort(i) = det_alpha_norm(i)
iorder(i) = i
enddo
do i=1,n_det_beta_unique
norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
iorder(i+n_det_alpha_unique) = -i
enddo
call dsort(norm_sort(1),iorder(1),nab)
if (.not.read_wf) then
stop 'Please set read_wf to true'
endif
PROVIDE psi_bilinear_matrix_values nuclear_repulsion
PROVIDE psi_bilinear_matrix_values nuclear_repulsion
print *, ''
print *, '=============================='
print *, 'Energies at different cut-offs'
@ -47,16 +45,21 @@ program e_curve
print *, '=========================================================='
double precision :: thresh
integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
integer :: na, nb
thresh = 1.d-10
na = n_det_alpha_unique
nb = n_det_beta_unique
do j=0,nab
i = iorder(j)
if (i<0) then
nb -= 1
do k=1,n_det
if (psi_bilinear_matrix_columns(k) == -i) then
psi_bilinear_matrix_values(k,1) = 0.d0
endif
enddo
else
na -= 1
do k=1,n_det
if (psi_bilinear_matrix_rows(k) == i) then
psi_bilinear_matrix_values(k,1) = 0.d0
@ -67,27 +70,11 @@ program e_curve
cycle
endif
u_0 = psi_bilinear_matrix_values(1:N_det,1:N_states)
v_t = 0.d0
s_t = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_states)
call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_states, N_det)
double precision, external :: u_dot_u, u_dot_v
do i=1,N_states
e_0(i) = u_dot_v(v_t(1,i),u_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
do k=1,N_states
psi_coef(1:N_det,k) = psi_bilinear_matrix_values(1:N_det,k)
call dset_order(psi_coef(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
TOUCH psi_det psi_coef
m = 0
do k=1,n_det
@ -100,10 +87,18 @@ program e_curve
exit
endif
E = E_0(1) + nuclear_repulsion
norm = u_dot_u(u_0(1,1),N_det)
print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)', thresh, m, &
dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, E
double precision :: cost0, cost
cost0 = elec_alpha_num**3 + elec_beta_num**3
cost = (na-1) * elec_alpha_num**2 + &
(nb-1) * elec_beta_num**2 + &
elec_alpha_num**3 + elec_beta_num**3
cost = cost/cost0
double precision :: u_dot_u
norm = dsqrt(u_dot_u(psi_coef(1,1),N_det))
print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)', thresh, m, &
cost, norm, psi_energy(1)
thresh = thresh * dsqrt(10.d0)
enddo
print *, '=========================================================='

75
devel/qmcchem/qmc_e_curve2.irp.f Normal file
View File

@ -0,0 +1,75 @@
program e_curve
use bitmasks
implicit none
integer :: i,j,k, kk, nab, m, l
double precision :: norm, E, hij, num, ci, cj
double precision :: e_0(N_states)
double precision, allocatable :: psi_save(:)
PROVIDE mo_two_e_integrals_in_map mo_one_e_integrals
if (.not.read_wf) then
stop 'Please set read_wf to true'
endif
PROVIDE psi_bilinear_matrix nuclear_repulsion
PROVIDE psi_coef_sorted psi_det psi_coef
print *, ''
print *, '=============================='
print *, 'Energies at different cut-offs'
print *, '=============================='
print *, ''
print *, '=========================================================='
print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)', 'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
print *, '=========================================================='
double precision :: thresh
integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
thresh = 1.d-10
integer :: n_det_prime
nab = n_det_alpha_unique+n_det_beta_unique
allocate(psi_save(1:N_det))
psi_save(1:N_det) = psi_coef(1:N_det,1)
do while (thresh < 1.d0)
norm = 0.d0
n_det_prime = n_det
do k=1,n_det
psi_coef(k,1) = psi_save(k)
if (dabs(psi_coef(k,1)) < thresh) then
psi_coef(k,1) = 0.d0
n_det_prime -= 1
endif
norm = norm + psi_coef(k,1)**2
enddo
TOUCH psi_coef psi_det
norm = norm/dsqrt(norm)
psi_coef(:,1) = psi_coef(:,1)/norm
integer :: na, nb
na = 0
do k=1,n_det_alpha_unique
if (det_alpha_norm(k) > 0.d0) then
na += 1
endif
enddo
nb = 0
do k=1,n_det_beta_unique
if (det_beta_norm(k) > 0.d0) then
nb += 1
endif
enddo
E = psi_energy(1)
double precision :: cost0, cost
cost0 = elec_alpha_num**3 + elec_beta_num**3
cost = (na-1) * elec_alpha_num**2 + &
(nb-1) * elec_beta_num**2 + &
elec_alpha_num**3 + elec_beta_num**3
cost = cost/cost0
print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)', thresh, n_det_prime, &
cost, norm, E
thresh = thresh * dsqrt(10.d0)
enddo
print *, '=========================================================='
end

166
devel/svdwf/Evar_TruncSVD.irp.f Normal file
View File

@ -0,0 +1,166 @@
program Evar_TruncSVD
implicit none
BEGIN_DOC
! study energy variation with truncated SVD
END_DOC
read_wf = .True.
TOUCH read_wf
! !!!
call run()
! !!!
end
subroutine run
implicit none
include 'constants.include.F'
double precision, allocatable :: A(:,:), U(:,:), V(:,:), D(:)
integer :: r, i, j, k, l, m, n, iter, iter_max
double precision, allocatable :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:), r1(:,:)
! !!!
m = n_det_alpha_unique
n = n_det_beta_unique
r = n
print *, 'matrix:', m,'x',n
print *, 'N det:', N_det
print *, 'rank = ', r
iter_max = 20
! !!!
allocate( Z(m,r) , P(n,r) ) ! Z(m,r) = A(m,n) @ P(n,r)
Z(:,:) = 0.d0
! !!!
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
! first we apply a RSVD for a pre-fixed rank (r)
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,r1)
allocate(r1(N_det,2))
!$OMP DO
do l=1,r
call random_number(r1)
r1(:,1) = dsqrt(-2.d0*dlog(r1(:,1)))
r1(:,1) = r1(:,1) * dcos(dtwo_pi*r1(:,2))
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1)
enddo
enddo
!$OMP END DO
deallocate(r1)
!$OMP END PARALLEL
! !!!
! Power iterations
! !!!
do iter=1,iter_max
! !!!
print *, 'Power iteration ', iter, '/', 20
! !!!
! P(n,r) = At(n,m) @ Z(m,r)
! !!!
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
!$OMP DO
do l=1,r
P(:,l) = 0.d0
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
P(j,l) = P(j,l) + psi_bilinear_matrix_values(k,1) * Z(i,l)
enddo
enddo
!$OMP END DO
! !!!
! Z(m,r) = A(m,n) @ P(n,r)
! !!!
!$OMP BARRIER
!$OMP DO
do l=1,r
Z(:,l) = 0.d0
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * P(j,l)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! !!!
! Compute QR: at return: Q is Z(m,r)
! !!!
call ortho_qr(Z,size(Z,1),m,r)
! !!!
enddo
! !!!
! Y(r,n) = Zt(r,m) @ A(m,n) or Yt(n,r) = At(n,m) @ Z(m,r)
! !!!
allocate(Yt(n,r))
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
!$OMP DO
do l=1,r
do k=1,n
Yt(k,l) = 0.d0
enddo
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Yt(j,l) = Yt(j,l) + Z(i,l) * psi_bilinear_matrix_values(k,1)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! !!!
! Y = UY @ D @ Vt or Yt = V @ Dt @ UYt
! !!!
allocate(D(r),V(n,r), UYt(r,r))
! !!!
call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r)
deallocate(Yt)
! !!!
! U(m,r) = Z(m,r) @ UY(r,r) or U = Z @ (UYt).T
! !!!
allocate(U(m,r))
call dgemm('N','T',m,r,r,1.d0,Z,size(Z,1),UYt,size(UYt,1),0.d0,U,size(U,1))
deallocate(UYt,Z)
! !!!
!do i=1,r
! print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
! if (D(i) < 1.d-15) then
! k = i
! exit
! endif
!enddo
!print *, 'threshold: ', 2.858 * D(k/2)
! !!!
! Build the new determinant: U @ D @ Vt
! !!!
!!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
!!$OMP DO
!!
!print *, 'ok 1'
!N_det = m * n
!print *, 'ok 11'
!TOUCH N_det
!psi_bilinear_matrix_values(:,1) = 0.d0
!TOUCH psi_bilinear_matrix_values
! print *, size(psi_bilinear_matrix_values,1), size(D), size(U,1), size(U,2), size(V,1), size(V,2)
print*, PSI_energy(1) + nuclear_repulsion
psi_bilinear_matrix(:,:,:) = 0.d0
do r = 1, n
call generate_all_alpha_beta_det_products
do i = 1, N_det_beta_unique
do j = 1, N_det_alpha_unique
psi_bilinear_matrix(j,i,1) = 0.d0
do l = 1, r
psi_bilinear_matrix(j,i,1) = psi_bilinear_matrix(j,i,1) + D(l) * U(j,l) * V(i,l)
enddo
enddo
enddo
TOUCH psi_bilinear_matrix
call update_wf_of_psi_bilinear_matrix(.False.)
print*, r, PSI_energy(1) + nuclear_repulsion, s2_values(1) !CI_energy(1)
call save_wavefunction()
enddo
deallocate(U,D,V)
! !!!
end

1
devel/svdwf/NEED
View File

@ -1 +1,2 @@
determinants
davidson_undressed

13
devel/svdwf/ab_products.irp.f Normal file
View File

@ -0,0 +1,13 @@
program abproducts
implicit none
read_wf = .True.
TOUCH read_wf
call run
end
subroutine run
implicit none
call generate_all_alpha_beta_det_products
call diagonalize_ci
call save_wavefunction
end

465
devel/svdwf/buildpsi_diagSVDit_modif_v2.irp.f Normal file
View File

@ -0,0 +1,465 @@
program buildpsi_diagSVDit_modif_v2
implicit none
BEGIN_DOC
! perturbative approach to build psi_postsvd
END_DOC
read_wf = .True.
TOUCH read_wf
PROVIDE N_int
call run()
end
subroutine run
USE OMP_LIB
implicit none
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: degree, i_state
integer :: i, j, k, l, ii, jj, na, nb
double precision :: norm_psi, inv_sqrt_norm_psi
double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:)
double precision :: err0, err_tmp, e_tmp, E0, overlap, E0_old, tol_energy
double precision :: ctmp, htmp, Ept2
double precision :: E0_postsvd, overlap_postsvd, E_prev
double precision :: norm_coeff_psi, inv_sqrt_norm_coeff_psi
double precision :: overlapU, overlapU_mat, overlapV, overlapV_mat, overlap_psi
double precision, allocatable :: H_diag(:,:), Hkl(:,:), H0(:,:), H(:,:,:,:)
double precision, allocatable :: psi_postsvd(:,:), coeff_psi_perturb(:)
integer :: n_FSVD, n_selected, n_toselect, n_tmp, it_svd, it_svd_max
integer :: n_selected2
integer, allocatable :: numalpha_selected(:), numbeta_selected(:)
integer, allocatable :: numalpha_toselect(:), numbeta_toselect(:)
integer, allocatable :: numalpha_tmp(:), numbeta_tmp(:)
integer(kind=8) :: W_tbeg, W_tend, W_tbeg_it, W_tend_it, W_ir
real(kind=8) :: W_tot_time, W_tot_time_it
real(kind=8) :: CPU_tbeg, CPU_tend, CPU_tbeg_it, CPU_tend_it
real(kind=8) :: CPU_tot_time, CPU_tot_time_it
real(kind=8) :: speedup, speedup_it
integer :: nb_taches
!$OMP PARALLEL
nb_taches = OMP_GET_NUM_THREADS()
!$OMP END PARALLEL
call CPU_TIME(CPU_tbeg)
call SYSTEM_CLOCK(COUNT=W_tbeg, COUNT_RATE=W_ir)
i_state = 1
! ---------------------------------------------------------------------------------------
! construct the initial CISD matrix
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
print *, ' CI matrix:', n_det_alpha_unique,'x',n_det_beta_unique
print *, ' N det :', N_det
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
norm_psi = 0.d0
do k = 1, N_det
norm_psi = norm_psi + psi_bilinear_matrix_values(k,i_state) &
* psi_bilinear_matrix_values(k,i_state)
enddo
print *, ' initial norm = ', norm_psi
allocate( Aref(n_det_alpha_unique,n_det_beta_unique) )
Aref(:,:) = 0.d0
do k = 1, N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Aref(i,j) = psi_bilinear_matrix_values(k,i_state)
enddo
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! perform a Full SVD
allocate( Uref(n_det_alpha_unique,n_det_alpha_unique) )
allocate( Dref(max(n_det_beta_unique,n_det_alpha_unique)) )
allocate( Vref(n_det_beta_unique,n_det_beta_unique) )
allocate( Vtref(n_det_beta_unique,n_det_beta_unique) )
call svd_s(Aref, size(Aref,1), Uref, size(Uref,1), Dref, Vtref, size(Vtref,1) &
, n_det_alpha_unique, n_det_beta_unique)
print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
do l = 1, n_det_beta_unique
do i = 1, n_det_beta_unique
Vref(i,l) = Vtref(l,i)
enddo
enddo
deallocate( Vtref )
! Truncated rank
n_selected = 10 ! n_svd
call write_int(6,n_selected, 'Rank of psi')
!________________________________________________________________________________________________________
!
! loop over SVD iterations
!________________________________________________________________________________________________________
tol_energy = 1.d0
it_svd = 0
it_svd_max = 100
E_prev = 0.d0
allocate(H(n_selected,n_selected,n_det_alpha_unique,n_det_beta_unique))
allocate(H_diag(n_det_alpha_unique,n_det_beta_unique))
allocate(psi_postsvd(n_det_alpha_unique,n_det_beta_unique))
do while( ( it_svd .lt. it_svd_max) .and. ( tol_energy .gt. 1d-6 ) )
call CPU_TIME(CPU_tbeg_it)
call SYSTEM_CLOCK(COUNT=W_tbeg_it, COUNT_RATE=W_ir)
it_svd = it_svd + 1
double precision :: norm
norm = 0.d0
do j = 1, n_selected
norm = norm + Dref(j)*Dref(j)
enddo
Dref = Dref / dsqrt(norm)
! print *, '-- Compute H --'
call const_H_uv(Uref, Vref, H, H_diag, n_selected)
! H0(i,j) = < u_i v_j | H | u_i v_j >
! E0 = < psi_0 | H | psi_0 >
E0 = 0.d0
do j = 1, n_selected
do i = 1, n_selected
E0 = E0 + Dref(i) * H(i,i,j,j) * Dref(j)
enddo
enddo
E0 = E0 + nuclear_repulsion
! print *,' E0 =', E0
double precision, allocatable :: eigval0(:)
double precision, allocatable :: eigvec0(:,:,:)
double precision, allocatable :: H_tmp(:,:,:,:)
allocate( H_tmp(n_selected,n_selected,n_selected,n_selected) )
do l=1,n_selected
do k=1,n_selected
do j=1,n_selected
do i=1,n_selected
H_tmp(i,j,k,l) = H(i,j,k,l)
enddo
enddo
enddo
enddo
allocate( eigval0(n_selected**2),eigvec0(n_selected,n_selected,n_selected**2))
eigvec0 = 0.d0
! print *, ' --- Diag post-SVD --- '
call lapack_diag(eigval0, eigvec0, H_tmp, n_selected**2, n_selected**2)
! print *, 'eig =', eigval0(1) + nuclear_repulsion
deallocate(H_tmp, eigval0)
! print *, ' --- SVD --- '
Dref = 0.d0
call perform_newpostSVD(n_selected, eigvec0(1,1,1), Uref, Vref, Dref)
deallocate(eigvec0)
! print *, ' --- Compute H --- '
call const_H_uv(Uref, Vref, H, H_diag, n_selected)
! H0(i,j) = < u_i v_j | H | u_i v_j >
! E0 = < psi_0 | H | psi_0 >
E0 = 0.d0
norm = 0.d0
do j = 1, n_det_beta_unique
do i = 1, n_selected
E0 = E0 + Dref(i) * H(i,i,j,j) * Dref(j)
enddo
norm = norm + Dref(j)*Dref(j)
enddo
E0 = E0 + nuclear_repulsion
! print *,' E0 =', E0
! print *,' norm =', norm
! print *, ' --- Perturbation --- '
psi_postsvd = 0.d0
do i=1,n_selected
psi_postsvd(i,i) = Dref(i)
enddo
Ept2 = 0.d0
do j=1,n_selected
do i=n_selected+1,n_det_alpha_unique
ctmp = 0.d0
do l=1,n_selected
do k=1,n_selected
ctmp = ctmp + H(k,l,i,j) * psi_postsvd(k,l)
enddo
enddo
psi_postsvd(i,j) = ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
Ept2 += ctmp*ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
enddo
enddo
do j=n_selected+1,n_det_beta_unique
do i=1,n_selected
ctmp = 0.d0
do l=1,n_selected
do k=1,n_selected
ctmp = ctmp + H(k,l,i,j) * psi_postsvd(k,l)
enddo
enddo
psi_postsvd(i,j) = ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
Ept2 += ctmp*ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
enddo
enddo
tol_energy = dabs(E_prev - E0)
print '(I5, 3(3X, F20.10))', it_svd, E0, E0 + Ept2, tol_energy
E_prev = E0
! print *, ' --- SVD --- '
call perform_newpostSVD(n_det_beta_unique, psi_postsvd, Uref, Vref, Dref)
end do
end
subroutine perform_newpostSVD(n_selected, psi_postsvd, Uref, Vref, Dref)
USE OMP_LIB
integer, intent(in) :: n_selected
double precision, intent(in) :: psi_postsvd(n_selected,n_selected)
double precision, intent(inout) :: Uref(n_det_alpha_unique,n_det_alpha_unique)
double precision, intent(inout) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(inout) :: Dref(max(n_det_beta_unique,n_det_alpha_unique))
integer :: mm, nn, i, j, ii0, ii, l, jj, na, nb
double precision :: err0, err_norm, err_tmp, norm_tmp
double precision :: overlapU_mat, overlapV_mat, overlapU, overlapV
double precision, allocatable :: S_mat(:,:), SxVt(:,:)
double precision, allocatable :: U_svd(:,:), V_svd(:,:)
double precision, allocatable :: U_newsvd(:,:), V_newsvd(:,:), Vt_newsvd(:,:), D_newsvd(:), A_newsvd(:,:)
mm = n_det_alpha_unique
nn = n_det_beta_unique
allocate( U_svd(mm,n_selected) , V_svd(nn,n_selected) , S_mat(n_selected,n_selected) )
U_svd(1:mm,1:n_selected) = Uref(1:mm,1:n_selected)
V_svd(1:nn,1:n_selected) = Vref(1:nn,1:n_selected)
S_mat(1:n_selected,1:n_selected) = psi_postsvd(1:n_selected,1:n_selected)
! first compute S_mat x transpose(V_svd)
allocate( SxVt(n_selected,nn) )
call dgemm( 'N', 'T', n_selected, nn, n_selected, 1.d0 &
, S_mat , size(S_mat,1) &
, V_svd , size(V_svd,1) &
, 0.d0, SxVt, size(SxVt ,1) )
deallocate(S_mat)
! then compute U_svd x SxVt
allocate( A_newsvd(mm,nn) )
call dgemm( 'N', 'N', mm, nn, n_selected, 1.d0 &
, U_svd , size(U_svd ,1) &
, SxVt , size(SxVt ,1) &
, 0.d0, A_newsvd, size(A_newsvd,1) )
deallocate( SxVt )
! perform new SVD
allocate( U_newsvd(mm,mm), Vt_newsvd(nn,nn), D_newsvd(max(mm,nn)) )
call svd_s( A_newsvd, size(A_newsvd,1), &
U_newsvd, size(U_newsvd,1), &
D_newsvd, &
Vt_newsvd, size(Vt_newsvd,1), &
mm, nn)
deallocate(A_newsvd)
allocate( V_newsvd(nn,nn) )
do l = 1, nn
do j = 1, nn
V_newsvd(j,l) = Vt_newsvd(l,j)
enddo
enddo
deallocate(Vt_newsvd)
do l = 1, n_selected
Dref(l) = D_newsvd(l)
Uref(1:mm,l) = U_newsvd(1:mm,l)
Vref(1:nn,l) = V_newsvd(1:nn,l)
enddo
deallocate(U_newsvd)
deallocate(V_newsvd)
deallocate(D_newsvd)
end subroutine perform_newpostSVD
subroutine const_H_uv(Uref, Vref, H, H_diag, n_selected)
USE OMP_LIB
implicit none
integer, intent(in) :: n_selected
double precision, intent(in) :: Uref(n_det_alpha_unique,n_det_alpha_unique)
double precision, intent(in) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(out) :: H(n_selected,n_selected, n_det_alpha_unique, n_det_beta_unique)
double precision, intent(out) :: H_diag(n_det_alpha_unique,n_det_beta_unique)
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: i, j, k, l, degree
integer :: ii0, jj0, ii, jj, n, m, np, mp
integer :: nn0, mm0, na, nb, mm, ind_gs
integer :: p,q,r,s
double precision :: h12, x
double precision, allocatable :: H0(:,:,:,:)
double precision, allocatable :: H1(:,:,:,:)
double precision, allocatable :: tmp3(:,:,:)
double precision, allocatable :: tmp1(:,:), tmp0(:,:)
double precision :: c_tmp
na = n_det_alpha_unique
nb = n_det_beta_unique
det1(:,1) = psi_det_alpha_unique(:,1)
det2(:,1) = psi_det_alpha_unique(:,1)
det1(:,2) = psi_det_beta_unique(:,1)
det2(:,2) = psi_det_beta_unique(:,1)
call i_H_j(det1, det2, N_int, h12)
call wall_time(t0)
tmp3 = 0.d0
allocate( H0(na,nb,n_selected,n_selected) )
allocate (tmp3(nb,nb,nb))
H0 = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,k,l,m,n,det1,det2,degree,h12,c_tmp,tmp1,tmp0)&
!$OMP SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique,&
!$OMP N_int,tmp3,Uref,Vref,H_diag,H0,n_selected)
allocate(tmp1(na,na), tmp0(na,na))
do i=1,na
do m=1,na
tmp1(m,i) = Uref(i,m)
enddo
enddo
!$OMP DO
do l = 1, nb
det2(:,2) = psi_det_beta_unique(:,l)
do j = 1, nb
det1(:,2) = psi_det_beta_unique(:,j)
call get_excitation_degree_spin(det1(1,2),det2(1,2),degree,N_int)
if (degree > 2) cycle
do k = 1, na
det2(:,1) = psi_det_alpha_unique(:,k)
do i = 1, na
det1(:,1) = psi_det_alpha_unique(:,i)
call get_excitation_degree(det1,det2,degree,N_int)
if ( degree > 2) cycle
call i_H_j(det1, det2, N_int, h12)
if (h12 == 0.d0) cycle
do m=1,nb
tmp3(m,j,l) = tmp3(m,j,l) + h12 * tmp1(m,i) * tmp1(m,k)
enddo
do n=1,n_selected
c_tmp = h12 * Vref(j,n)
if (c_tmp == 0.d0) cycle
do m=1,n_selected
H0(k,l,m,n) = H0(k,l,m,n) + c_tmp * tmp1(m,i)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do m=1,nb
do l=1,nb
do j=1,nb
tmp1(j,l) = tmp3(m,j,l)
enddo
enddo
call DGEMM('N','N',nb,nb,nb,1.d0, &
tmp1, size(tmp1,1), &
Vref, size(Vref,1), &
0.d0, tmp0, size(tmp0,1))
do n=1,nb
H_diag(m,n) = 0.d0
do j=1,nb
H_diag(m,n) = H_diag(m,n) + tmp0(j,n) * Vref(j,n)
enddo
enddo
enddo
!$OMP END DO
deallocate(tmp1, tmp0)
!$OMP END PARALLEL
call wall_time(t1)
allocate( H1(nb,n_selected,n_selected,na) )
call DGEMM('T','N', nb * n_selected * n_selected, na, na, &
1.d0, H0, size(H0,1), Uref, size(Uref,1), 0.d0, H1, size(H1,1)*size(H1,2)*size(H1,3))
deallocate( H0 )
! (l,p,q,r) -> (p,q,r,s)
call DGEMM('T','N', n_selected * n_selected * na, nb, nb, &
1.d0, H1, size(H1,1), Vref, size(Vref,1), 0.d0, H, size(H,1)*size(H,2)*size(H,3))
! do j=1,n_selected
! do i=1,n_selected
! print *, H_diag(i,j), H(i,j,i,j)
! enddo
! enddo
deallocate(H1)
call wall_time(t2)
! print *, 't=', t1-t0, t2-t1
double precision :: t0, t1, t2
! stop
end

18
devel/svdwf/random_svd.irp.f
View File

@ -22,13 +22,17 @@ subroutine run
allocate(Z(m,r))
! Z(m,r) = A(m,n).P(n,r)
Z(:,:) = 0.d0
! Z(m,r) = A(m,n).P(n,r)
do j=1,r
do i=1,m
Z(i,j) = 0.d0
enddo
enddo
allocate(P(n,r))
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,r1)
allocate(r1(N_det,2))
!$OMP DO
!$OMP DO
do l=1,r
call random_number(r1)
r1(:,1) = dsqrt(-2.d0*dlog(r1(:,1)))
@ -36,7 +40,7 @@ subroutine run
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1)
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1(k,1)
enddo
enddo
!$OMP END DO
@ -84,7 +88,9 @@ subroutine run
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
!$OMP DO
do l=1,r
Yt(:,l) = 0.d0
do k=1,n
Yt(k,l) = 0.d0
enddo
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
@ -98,7 +104,7 @@ subroutine run
call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r)
deallocate(Yt)
! U(m,r) = Z(m,r).UY(r,r)
! U(m,r) = Z(m,r).UY(r,r)
allocate(U(m,r))
call dgemm('N','T',m,r,r,1.d0,Z,size(Z,1),UYt,size(UYt,1),0.d0,U,size(U,1))
deallocate(UYt,Z)

41
devel/svdwf/svdwf.irp.f
View File

@ -12,36 +12,57 @@ subroutine run
implicit none
double precision, allocatable :: U(:,:), Vt(:,:), D(:), A(:,:)
integer :: i, j, k, p, q
double precision :: entropy
allocate( A (n_det_alpha_unique, n_det_beta_unique), &
U (n_det_alpha_unique, n_det_alpha_unique), &
Vt(n_det_beta_unique, n_det_beta_unique), &
D(max(n_det_beta_unique,n_det_alpha_unique)) )
A = 0.D0
do j=1,n_det_beta_unique
do i=1,n_det_alpha_unique
A(i,j) = 0.D0
enddo
enddo
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
A(i,j) = psi_bilinear_matrix_values(k,1)
enddo
call randomized_svd(A, size(A,1), &
if (N_det == 1) then
D(1) = 1.d0
U(1,1) = 1.d0
Vt(1,1) = 1.d0
else
call randomized_svd(A, size(A,1), &
U, size(U,1), D, Vt, size(Vt,1), n_det_alpha_unique, n_det_beta_unique, &
6,1000)
6,min(n_det_beta_unique,1000))
endif
entropy = 0.d0
do i=1,n_det_beta_unique
print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
entropy -= D(i) * dlog(D(i))
if (D(i) < 1.d-15) then
k = i
exit
endif
enddo
print *, 'threshold: ', 2.858 * D(k/2)
do i=1,n_det_alpha_unique
print '(I6,4(X,F12.8))', i, U(i,1:4)
enddo
print *, ''
do i=1,n_det_beta_unique
print '(I6,4(X,F12.8))', i, Vt(1:4,i)
enddo
print *, 'Entropy : ', entropy
call ezfio_set_spindeterminants_psi_svd_alpha(U)
call ezfio_set_spindeterminants_psi_svd_beta (Vt)
call ezfio_set_spindeterminants_psi_svd_coefs(D)
! do i=1,n_det_alpha_unique
! print '(I6,4(X,F12.8))', i, U(i,1:4)
! enddo
! print *, ''
! do i=1,n_det_beta_unique
! print '(I6,4(X,F12.8))', i, Vt(1:4,i)
! enddo
end

12
devel/trexio/EZFIO.cfg Normal file
View File

@ -0,0 +1,12 @@
[backend]
type: integer
doc: Back-end used in TREXIO. 0: HDF5, 1:Text
interface: ezfio, ocaml, provider
default: 0
[trexio_file]
type: character*(256)
doc: Name of the exported TREXIO file
interface: ezfio, ocaml, provider
default: None

2
devel/trexio/LIB Normal file
View File

@ -0,0 +1,2 @@
-L/home/scemama/TREX/trexio/_install/lib -ltrexio

6
devel/trexio/NEED Normal file
View File

@ -0,0 +1,6 @@
ezfio_files
determinants
mo_one_e_ints
mo_two_e_ints
ao_two_e_ints
ao_one_e_ints

4
devel/trexio/README.rst Normal file
View File

@ -0,0 +1,4 @@
======
trexio
======

262
devel/trexio/export_trexio.irp.f Normal file
View File

@ -0,0 +1,262 @@
program export_trexio
use trexio
implicit none
BEGIN_DOC
! Exports the wave function in TREXIO format
END_DOC
integer(8) :: f ! TREXIO file handle
integer :: rc
print *, 'TREXIO file : '//trim(trexio_filename)
print *, ''
if (backend == 0) then
f = trexio_open(trexio_filename, 'w', TREXIO_HDF5)
else if (backend == 1) then
f = trexio_open(trexio_filename, 'w', TREXIO_TEXT)
endif
if (f == 0) then
print *, 'Unable to open TREXIO file for writing'
stop -1
endif
! ------------------------------------------------------------------------------
! Electrons
! ---------
rc = trexio_write_electron_up_num(f, elec_alpha_num)
call check_success(rc)
rc = trexio_write_electron_dn_num(f, elec_beta_num)
call check_success(rc)
! Nuclei
! ------
rc = trexio_write_nucleus_num(f, nucl_num)
call check_success(rc)
rc = trexio_write_nucleus_charge(f, nucl_charge)
call check_success(rc)
rc = trexio_write_nucleus_coord(f, nucl_coord_transp)
call check_success(rc)
rc = trexio_write_nucleus_label(f, nucl_label, 32)
call check_success(rc)
! Pseudo-potentials
! -----------------
double precision, allocatable :: tmp_double(:,:)
integer, allocatable :: tmp_int(:,:)
! rc = trexio_write_ecp_lmax_plus_1(f, pseudo_lmax+1)
! call check_success(rc)
!
! rc = trexio_write_ecp_z_core(f, nucl_charge_remove)
! call check_success(rc)
!
! rc = trexio_write_ecp_local_num_n_max(f, pseudo_klocmax)
! call check_success(rc)
!
! rc = trexio_write_ecp_local_power(f, pseudo_n_k_transp)
! call check_success(rc)
!
! rc = trexio_write_ecp_local_exponent(f, pseudo_dz_k_transp)
! call check_success(rc)
!
! rc = trexio_write_ecp_local_coef(f, pseudo_v_k_transp)
! call check_success(rc)
!
! rc = trexio_write_ecp_non_local_num_n_max(f, pseudo_kmax)
! call check_success(rc)
!
! rc = trexio_write_ecp_non_local_power(f, pseudo_n_kl_transp)
! call check_success(rc)
!
! rc = trexio_write_ecp_non_local_exponent(f, pseudo_dz_kl_transp)
! call check_success(rc)
!
! rc = trexio_write_ecp_non_local_coef(f, pseudo_v_kl_transp)
! call check_success(rc)
! Basis
! -----
rc = trexio_write_basis_type(f, 'Gaussian', len('Gaussian'))
call check_success(rc)
rc = trexio_write_basis_num(f, shell_num)
call check_success(rc)
rc = trexio_write_basis_nucleus_shell_num(f, nucleus_shell_num)
call check_success(rc)
rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
call check_success(rc)
rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
call check_success(rc)
rc = trexio_write_basis_prim_num(f, prim_num)
call check_success(rc)
rc = trexio_write_basis_shell_prim_num(f, shell_prim_num)
call check_success(rc)
double precision, allocatable :: factor(:)
allocate(factor(shell_num))
if (ao_normalized) then
factor(1:shell_num) = shell_normalization_factor(1:shell_num)
else
factor(1:shell_num) = 1.d0
endif
rc = trexio_write_basis_shell_factor(f, factor)
call check_success(rc)
deallocate(factor)
rc = trexio_write_basis_shell_prim_index(f, shell_prim_index)
call check_success(rc)
rc = trexio_write_basis_exponent(f, prim_expo)
call check_success(rc)
rc = trexio_write_basis_coefficient(f, prim_coef)
call check_success(rc)
allocate(factor(prim_num))
if (primitives_normalized) then
factor(1:prim_num) = prim_normalization_factor(1:prim_num)
else
factor(1:prim_num) = 1.d0
endif
rc = trexio_write_basis_prim_factor(f, factor)
call check_success(rc)
deallocate(factor)
! Atomic orbitals
! ---------------
rc = trexio_write_ao_num(f, ao_num)
call check_success(rc)
rc = trexio_write_ao_cartesian(f, 1)
call check_success(rc)
rc = trexio_write_ao_shell(f, ao_shell)
call check_success(rc)
integer :: i
allocate(factor(ao_num))
if (ao_normalized) then
do i=1,ao_num
factor(i) = ao_coef_normalization_factor(i) / shell_normalization_factor( ao_shell(i) )
enddo
else
factor(:) = 1.d0
endif
rc = trexio_write_ao_normalization(f, factor)
call check_success(rc)
deallocate(factor)
! One-e AO integrals
! ------------------
rc = trexio_write_ao_1e_int_overlap(f,ao_overlap)
call check_success(rc)
rc = trexio_write_ao_1e_int_kinetic(f,ao_kinetic_integrals)
call check_success(rc)
rc = trexio_write_ao_1e_int_potential_n_e(f,ao_integrals_n_e)
call check_success(rc)
if (do_pseudo) then
rc = trexio_write_ao_1e_int_ecp_local(f,ao_pseudo_integrals_local)
call check_success(rc)
rc = trexio_write_ao_1e_int_ecp_non_local(f,ao_pseudo_integrals_non_local)
call check_success(rc)
endif
rc = trexio_write_ao_1e_int_core_hamiltonian(f,ao_one_e_integrals)
call check_success(rc)
! Molecular orbitals
! ------------------
! rc = trexio_write_mo_type(f, mo_label)
! call check_success(rc)
rc = trexio_write_mo_num(f, mo_num)
call check_success(rc)
rc = trexio_write_mo_coefficient(f, mo_coef)
call check_success(rc)
! One-e MO integrals
! ------------------
rc = trexio_write_mo_1e_int_kinetic(f,mo_kinetic_integrals)
call check_success(rc)
rc = trexio_write_mo_1e_int_potential_n_e(f,mo_integrals_n_e)
call check_success(rc)
if (do_pseudo) then
rc = trexio_write_mo_1e_int_ecp_local(f,mo_pseudo_integrals_local)
call check_success(rc)
rc = trexio_write_mo_1e_int_ecp_non_local(f,mo_pseudo_integrals_non_local)
call check_success(rc)
endif
!
rc = trexio_write_mo_1e_int_core_hamiltonian(f,one_e_dm_mo)
call check_success(rc)
! RDM
! ----
! rc = trexio_write_rdm_one_e(f,one_e_dm_mo)
! call check_success(rc)
!
! rc = trexio_write_rdm_one_e_up(f,one_e_dm_mo_alpha_average)
! call check_success(rc)
!
! rc = trexio_write_rdm_one_e_dn(f,one_e_dm_mo_beta_average)
! call check_success(rc)
! ------------------------------------------------------------------------------
rc = trexio_close(f)
call check_success(rc)
end
subroutine check_success(rc)
use trexio
implicit none
integer, intent(in) :: rc
character*(128) :: str
if (rc /= TREXIO_SUCCESS) then
call trexio_string_of_error(rc,str)
print *, 'TREXIO Error: ' //trim(str)
stop -1
endif
end
! -*- mode: f90 -*-

3232
devel/trexio/trexio_f.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

20
devel/trexio/trexio_file.irp.f Normal file
View File

@ -0,0 +1,20 @@
BEGIN_PROVIDER [ character*(1024), trexio_filename ]
implicit none
BEGIN_DOC
! Name of the TREXIO file
END_DOC
character*(1024) :: prefix
trexio_filename = trexio_file
if (trexio_file == 'None') then
prefix = trim(ezfio_work_dir)//trim(ezfio_filename)
if (backend == 0) then
trexio_filename = trim(prefix)//'.h5'
else if (backend == 1) then
trexio_filename = trim(prefix)
endif
endif
END_PROVIDER

7
stable/amplitudes/EZFIO.cfg.save
View File

@ -1,7 +0,0 @@
[t1_amplitudes]
type: double precision
doc: Amplitudes for the single-excitation operator
interface: ezfio,provider
size: (mo_basis.mo_num,mo_basis.mo_num)

4
stable/utilities/hcore_guess.irp.f Normal file
View File

@ -0,0 +1,4 @@
program hcore_guess_prog
implicit none
call hcore_guess
end

8
stable/utilities/truncate_wf.irp.f
View File

@ -55,16 +55,16 @@ subroutine routine_s2
integer :: i,j,k
double precision :: accu(N_states)
print *, 'Weights of the SOP'
print *, 'Weights of the CFG'
do i=1,N_det
print *, i, real(weight_occ_pattern(det_to_occ_pattern(i),:)), real(sum(weight_occ_pattern(det_to_occ_pattern(i),:)))
print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
enddo
print*, 'Min weight of the occupation pattern ?'
read(5,*) wmin
ndet_max = 0
do i=1,N_det
if (maxval(weight_occ_pattern( det_to_occ_pattern(i),:)) < wmin) cycle
if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle
ndet_max = ndet_max+1
enddo
@ -73,7 +73,7 @@ subroutine routine_s2
accu = 0.d0
k=0
do i = 1, N_det
if (maxval(weight_occ_pattern( det_to_occ_pattern(i),:)) < wmin) cycle
if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle
k = k+1
do j = 1, N_int
psi_det_tmp(j,1,k) = psi_det(j,1,i)