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

Added Utils, Electrons, Ezfio, Nuclei, include, AOs.

Browse Source
This commit is contained in:
Anthony Scemama 2014-04-01 18:37:27 +02:00
parent 0c2db61521
commit 6c29d22aa2
23 changed files with 2056 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
README.md
View File

@ -4,5 +4,5 @@ Quantum package
Set of quantum chemistry programs and libraries.
For more information, you can visit the
`wiki of the project <http://github.com/LCPQ/quantum_package/wiki>`_
wiki of the project : http://github.com/LCPQ/quantum_package/wiki

9
src/AOs/Makefile Normal file
View File

@ -0,0 +1,9 @@
default: all
INCLUDE_DIRS = Ezfio_files Utils Nuclei
include ../Makefile.config
include ../Makefile.common
include irpf90.make
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO)
$(IRPF90)

10
src/AOs/aos.ezfio_config Normal file
View File

@ -0,0 +1,10 @@
ao_basis
ao_num integer
ao_prim_num integer (ao_basis_ao_num)
ao_nucl integer (ao_basis_ao_num)
ao_power integer (ao_basis_ao_num,3)
ao_prim_num_max integer = maxval(ao_basis_ao_prim_num)
ao_coef double precision (ao_basis_ao_num,ao_basis_ao_prim_num_max)
ao_expo double precision (ao_basis_ao_num,ao_basis_ao_prim_num_max)

192
src/AOs/aos.irp.f Normal file
View File

@ -0,0 +1,192 @@
BEGIN_PROVIDER [ integer, ao_num ]
&BEGIN_PROVIDER [ integer, ao_num_align ]
implicit none
BEGIN_DOC
! Number of atomic orbitals
END_DOC
ao_num = -1
PROVIDE ezfio_filename
call ezfio_get_ao_basis_ao_num(ao_num)
if (ao_num <= 0) then
stop 'Number of contracted gaussians should be > 0'
endif
integer :: align_double
ao_num_align = align_double(ao_num)
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_power, (ao_num_align,3) ]
&BEGIN_PROVIDER [ double precision, ao_expo, (ao_num_align,ao_prim_num_max) ]
&BEGIN_PROVIDER [ double precision, ao_coef, (ao_num_align,ao_prim_num_max) ]
implicit none
double precision, allocatable :: buffer(:,:)
allocate ( buffer(ao_num,ao_prim_num_max) )
integer :: ibuffer(ao_num,3)
integer :: i,j,k
ibuffer = 0
BEGIN_DOC
! coefficient of the primitives on the AO basis
! ao_coef(i,j) = coefficient of the jth primitive on the ith ao
END_DOC
PROVIDE ezfio_filename
call ezfio_get_ao_basis_ao_power(ibuffer)
ao_power = 0
do j = 1, 3
do i = 1, ao_num
ao_power(i,j) = ibuffer(i,j)
enddo
enddo
ao_expo = 0.d0
buffer = 0.d0
call ezfio_get_ao_basis_ao_expo(buffer)
do j = 1, ao_prim_num_max
do i = 1, ao_num
ao_expo(i,j) = buffer(i,j)
enddo
enddo
ao_coef = 0.d0
buffer = 0.d0
call ezfio_get_ao_basis_ao_coef(buffer)
do j = 1, ao_prim_num_max
do i = 1, ao_num
ao_coef(i,j) = buffer(i,j)
enddo
enddo
deallocate(buffer)
! Normalization of the AO coefficients
! ------------------------------------
double precision :: norm, norm2,overlap_x,overlap_y,overlap_z,C_A(3)
integer :: l, powA(3), nz
nz=100
C_A(1) = 0.d0
C_A(2) = 0.d0
C_A(3) = 0.d0
do i=1,ao_num
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
do j=1,ao_prim_num(i)
call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
ao_coef(i,j) = ao_coef(i,j)/sqrt(norm)
enddo
enddo
! Sorting of the exponents for efficient integral calculations
integer :: iorder(ao_prim_num_max)
double precision :: d(ao_prim_num_max,2)
do i=1,ao_num
do j=1,ao_prim_num(i)
iorder(j) = j
d(j,1) = ao_expo(i,j)
d(j,2) = ao_coef(i,j)
enddo
call dsort(d(1,1),iorder,ao_prim_num(i))
call dset_order(d(1,2),iorder,ao_prim_num(i))
do j=1,ao_prim_num(i)
ao_expo(i,j) = d(j,1)
ao_coef(i,j) = d(j,2)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap, (ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! matrix of the overlap for tha AOs
! S(i,j) = overlap between the ith and the jth atomic basis function
END_DOC
integer :: i,j,k,l,nz,num_i,num_j,powA(3),powB(3)
double precision :: accu,overlap_x,overlap_y,overlap_z,A_center(3),B_center(3),norm
nz=100
do i = 1, ao_num
num_i = ao_nucl(i)
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
A_center(1)=nucl_coord(num_i,1)
A_center(2)=nucl_coord(num_i,2)
A_center(3)=nucl_coord(num_i,3)
do j = i, ao_num
num_j = ao_nucl(j)
powB(1) = ao_power(j,1)
powB(2) = ao_power(j,2)
powB(3) = ao_power(j,3)
B_center(1)=nucl_coord(num_j,1)
B_center(2)=nucl_coord(num_j,2)
B_center(3)=nucl_coord(num_j,3)
accu = 0.d0
do k = 1, ao_prim_num(i)
do l = 1, ao_prim_num(j)
call overlap_gaussian_xyz(A_center,B_center,ao_expo(i,k),ao_expo(j,l),powA,powB,overlap_x,overlap_y,overlap_z,norm,nz)
accu = accu + norm * ao_coef(i,k) * ao_coef(j,l)
enddo
enddo
ao_overlap(i,j) = accu
ao_overlap(j,i) = accu
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ]
implicit none
BEGIN_DOC
! Transposed ao_coef and ao_expo
END_DOC
integer :: i,j
do j=1, ao_num
do i=1, ao_prim_num_max
ao_coef_transp(i,j) = ao_coef(j,i)
ao_expo_transp(i,j) = ao_expo(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_prim_num, (ao_num_align) ]
implicit none
BEGIN_DOC
! Number of primitives per atomic orbital
END_DOC
ao_prim_num = 0
PROVIDE ezfio_filename
call ezfio_get_ao_basis_ao_prim_num(ao_prim_num)
integer :: i
character*(80) :: message
do i=1,ao_num
if (ao_prim_num(i) <= 0) then
write(message,'(A,I6,A)') 'Number of primitives of contraction ',i,' should be > 0'
print *, message
stop
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_prim_num_max ]
&BEGIN_PROVIDER [ integer, ao_prim_num_max_align ]
implicit none
ao_prim_num_max = 0
PROVIDE ezfio_filename
call ezfio_get_ao_basis_ao_prim_num_max(ao_prim_num_max)
integer :: align_double
ao_prim_num_max_align = align_double(ao_prim_num_max)
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_nucl, (ao_num)]
BEGIN_DOC
! Index of the nuclei on which the ao is centered
END_DOC
implicit none
PROVIDE ezfio_filename
call ezfio_get_ao_basis_ao_nucl(ao_nucl)
END_PROVIDER

5
src/Electrons/electrons.ezfio_config Normal file
View File

@ -0,0 +1,5 @@
electrons
elec_alpha_num integer
elec_beta_num integer
elec_num integer = electrons_elec_alpha_num + electrons_elec_beta_num

21
src/Electrons/electrons.irp.f Normal file
View File

@ -0,0 +1,21 @@
BEGIN_PROVIDER [ integer, elec_alpha_num ]
&BEGIN_PROVIDER [ integer, elec_beta_num ]
&BEGIN_PROVIDER [ integer, elec_num ]
&BEGIN_PROVIDER [ integer, elec_num_tab, (2) ]
implicit none
BEGIN_DOC
! Numbers of alpha ("up") , beta ("down) and total electrons
END_DOC
PROVIDE ezfio_filename
call ezfio_get_electrons_elec_alpha_num(elec_alpha_num)
call ezfio_get_electrons_elec_beta_num(elec_beta_num)
call ezfio_get_electrons_elec_num(elec_num)
elec_num_tab(1) = elec_alpha_num
elec_num_tab(2) = elec_beta_num
ASSERT (elec_alpha_num > 0)
ASSERT (elec_beta_num >= 0)
END_PROVIDER

8
src/Ezfio_files/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
include ../Makefile.config
include ../Makefile.common
include irpf90.make
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO)
$(IRPF90)

15
src/Ezfio_files/ezfio.irp.f Normal file
View File

@ -0,0 +1,15 @@
BEGIN_PROVIDER [ character*(128), ezfio_filename ]
implicit none
BEGIN_DOC
! Name of EZFIO file
END_DOC
integer :: iargc
call getarg(0,ezfio_filename)
if (iargc() /= 1) then
print *, ezfio_filename, ' <ezfio_file>'
stop 1
endif
call getarg(1,ezfio_filename)
call ezfio_set_file(ezfio_filename)
END_PROVIDER

43
src/Ezfio_files/get_unit_and_open.irp.f Normal file
View File

@ -0,0 +1,43 @@
integer function getUnitAndOpen(f,mode)
implicit none
character*(*) :: f
character*(128) :: new_f
integer :: iunit
logical :: is_open, exists
character :: mode
is_open = .True.
iunit = 10
new_f = f
do while (is_open)
inquire(unit=iunit,opened=is_open)
if (.not.is_open) then
getUnitAndOpen = iunit
endif
iunit = iunit+1
enddo
if (mode.eq.'r') then
inquire(file=f,exist=exists)
if (.not.exists) then
open(unit=getUnitAndOpen,file=f,status='NEW',action='WRITE',form='FORMATTED')
close(unit=getUnitAndOpen)
endif
open(unit=getUnitAndOpen,file=f,status='OLD',action='READ',form='FORMATTED')
else if (mode.eq.'R') then
inquire(file=f,exist=exists)
if (.not.exists) then
open(unit=getUnitAndOpen,file=f,status='NEW',action='WRITE',form='UNFORMATTED')
close(unit=getUnitAndOpen)
endif
open(unit=getUnitAndOpen,file=f,status='OLD',action='READ',form='UNFORMATTED')
else if (mode.eq.'W') then
open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',form='UNFORMATTED')
else if (mode.eq.'w') then
open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',form='FORMATTED')
else if (mode.eq.'a') then
open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',position='APPEND',form='FORMATTED')
else if (mode.eq.'x') then
open(unit=getUnitAndOpen,file=new_f,form='FORMATTED')
endif
end function getUnitAndOpen

21
src/Makefile Normal file
View File

@ -0,0 +1,21 @@
default: all
INCLUDE_DIRS = Ezfio_files Nuclei Utils AOs Electrons # MonoInts BiInts MOs Output Bitmask
include Makefile.common
include Makefile.config
LIB+=$(MKL)
SRC=#BiInts/map_module.f90 Bitmask/bitmasks_module.f90
OBJ=#IRPF90_temp/BiInts/map_module.o IRPF90_temp/Bitmask/bitmasks_module.o
PYTHON_SCRIPTS=
include irpf90.make
all:$(ALL) $(PYTHON_SCRIPTS)
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO)
$(IRPF90)
all_clean:
for i in */ ; do $(MAKE) -C $$i veryclean clean_links ; done

39
src/Makefile.common Normal file
View File

@ -0,0 +1,39 @@
ifndef SCI_ROOT
$(error SCI_ROOT undefined. Run the setup_environment.sh script)
endif
IRP_VERSION_OK=$(shell $(IRPF90) -v | python -c "import sys ; print float(sys.stdin.readline().rsplit('.',1)[0]) >= 1.3")
ifeq ($(IRP_VERSION_OK),False)
$(error 'IRPF90 version >= 1.3 is required')
endif
MAKEFILE_OK=$(shell ls $(SCI_ROOT)/src/Makefile.config 2> /dev/null && echo True || echo False)
ifeq ($(MAKEFILE_OK),False)
$(error 'Makefile.config is not present. Please modify Makefile.config.example to create Makefile.config')
endif
EZFIO_DIR=$(SCI_ROOT)/EZFIO
EZFIO=$(EZFIO_DIR)/lib/libezfio_irp.a
$(EZFIO): $(wildcard $(SCI_ROOT)/src/*.ezfio_config) $(wildcard $(SCI_ROOT)/src/*/*.ezfio_config)
@echo Building EZFIO library
@cp $(wildcard $(SCI_ROOT)/src/*.ezfio_config) $(wildcard $(SCI_ROOT)/src/*/*.ezfio_config) $(EZFIO_DIR)/config
@cd $(EZFIO_DIR) ; export FC="$(FC)" ; export FCFLAGS="$(FCFLAGS)" ; export IRPF90="$(IRPF90)" ; $(MAKE) ; $(MAKE) Python
INCLUDE_DIRS+=include
ifneq ($(PWD),$(SCI_ROOT)/src)
$(shell $(SCI_ROOT)/scripts/prepare_module.sh $(INCLUDE_DIRS))
clean_links:
rm $(INCLUDE_DIRS) $$(basename $$PWD)
else
clean_links:
endif
LIB+=$(EZFIO) $(MKL)
IRPF90+=$(patsubst %, -I %, $(INCLUDE_DIRS))

29
src/Makefile.config.example Normal file
View File

@ -0,0 +1,29 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90_FLAGS+= --align=32
FC = ifort -g
FCFLAGS=
FCFLAGS+= -xHost
#FCFLAGS+= -xAVX
FCFLAGS+= -O2
FCFLAGS+= -ip
FCFLAGS+= -opt-prefetch
FCFLAGS+= -ftz
MKL=-mkl=parallel
ifeq ($(PROFILE),1)
FC += -p -g
CXX += -pg
endif
ifeq ($(OPENMP),1)
FC += -openmp
CXX += -fopenmp
endif
ifeq ($(DEBUG),1)
FC += -C -traceback -fpe0
#FCFLAGS =-O0
endif

9
src/Nuclei/Makefile Normal file
View File

@ -0,0 +1,9 @@
default: all
INCLUDE_DIRS = Ezfio_files Utils
include ../Makefile.config
include ../Makefile.common
include irpf90.make
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO)
$(IRPF90)

6
src/Nuclei/nuclei.ezfio_config Normal file
View File

@ -0,0 +1,6 @@
nuclei
nucl_num integer
nucl_label character*(32) (nuclei_nucl_num)
nucl_charge double precision (nuclei_nucl_num)
nucl_coord double precision (nuclei_nucl_num,3)

134
src/Nuclei/nuclei.irp.f Normal file
View File

@ -0,0 +1,134 @@
BEGIN_PROVIDER [ integer, nucl_num ]
&BEGIN_PROVIDER [ integer, nucl_num_aligned ]
implicit none
BEGIN_DOC
! Number of nuclei
END_DOC
PROVIDE ezfio_filename
nucl_num = 0
call ezfio_get_nuclei_nucl_num(nucl_num)
ASSERT (nucl_num > 0)
integer :: align_double
nucl_num_aligned = align_double(nucl_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_charge, (nucl_num) ]
implicit none
BEGIN_DOC
! Nuclear charges
END_DOC
PROVIDE ezfio_filename
nucl_charge = -1.d0
call ezfio_get_nuclei_nucl_charge(nucl_charge)
ASSERT (nucl_charge(:) >= 0.d0)
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num_aligned,3) ]
implicit none
BEGIN_DOC
! Nuclear coordinates in the format (:, {x,y,z})
END_DOC
PROVIDE ezfio_filename
double precision, allocatable :: buffer(:,:)
nucl_coord = 0.d0
allocate (buffer(nucl_num,3))
buffer = 0.d0
call ezfio_get_nuclei_nucl_coord(buffer)
integer :: i,j
do i=1,3
do j=1,nucl_num
nucl_coord(j,i) = buffer(j,i)
enddo
enddo
deallocate(buffer)
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_coord_transp, (3,nucl_num) ]
implicit none
BEGIN_DOC
! Transposed array of nucl_coord
END_DOC
integer :: i, k
nucl_coord_transp = 0.
do i=1,nucl_num
nucl_coord_transp(1,i) = nucl_coord(i,1)
nucl_coord_transp(2,i) = nucl_coord(i,2)
nucl_coord_transp(3,i) = nucl_coord(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_dist_2, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_x, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_y, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_z, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist, (nucl_num_aligned,nucl_num) ]
implicit none
BEGIN_DOC
! nucl_dist : Nucleus-nucleus distances
! nucl_dist_2 : Nucleus-nucleus distances squared
! nucl_dist_vec : Nucleus-nucleus distances vectors
END_DOC
integer :: ie1, ie2, l
integer,save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
nucl_dist = 0.d0
nucl_dist_2 = 0.d0
nucl_dist_vec_x = 0.d0
nucl_dist_vec_y = 0.d0
nucl_dist_vec_z = 0.d0
endif
do ie2 = 1,nucl_num
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do ie1 = 1,nucl_num_aligned
nucl_dist_vec_x(ie1,ie2) = nucl_coord(ie1,1) - nucl_coord(ie2,1)
nucl_dist_vec_y(ie1,ie2) = nucl_coord(ie1,2) - nucl_coord(ie2,2)
nucl_dist_vec_z(ie1,ie2) = nucl_coord(ie1,3) - nucl_coord(ie2,3)
enddo
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do ie1 = 1,nucl_num_aligned
nucl_dist_2(ie1,ie2) = nucl_dist_vec_x(ie1,ie2)*nucl_dist_vec_x(ie1,ie2) + &
nucl_dist_vec_y(ie1,ie2)*nucl_dist_vec_y(ie1,ie2) + &
nucl_dist_vec_z(ie1,ie2)*nucl_dist_vec_z(ie1,ie2)
nucl_dist(ie1,ie2) = sqrt(nucl_dist_2(ie1,ie2))
ASSERT (nucl_dist(ie1,ie2) > 0.d0)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
implicit none
BEGIN_DOC
! Nuclear repulsion energy
END_DOC
integer :: k,l
double precision :: Z12, r2, x(3)
nuclear_repulsion = 0.d0
do l = 1, nucl_num
do k = 1, nucl_num
if(k /= l) then
Z12 = nucl_charge(k)*nucl_charge(l)
x(1) = nucl_coord(k,1) - nucl_coord(l,1)
x(2) = nucl_coord(k,2) - nucl_coord(l,2)
x(3) = nucl_coord(k,3) - nucl_coord(l,3)
r2 = x(1)*x(1) + x(2)*x(2) + x(3)*x(3)
nuclear_repulsion += Z12/dsqrt(r2)
endif
enddo
enddo
nuclear_repulsion *= 0.5d0
END_PROVIDER

3
src/Nuclei/test.irp.f Normal file
View File

@ -0,0 +1,3 @@
program test
print *, nuclear_repulsion
end

195
src/Utils/LinearAlgebra.irp.f Normal file
View File

@ -0,0 +1,195 @@
subroutine ortho_lowdin(overlap,lda,n,C,ldc,m)
implicit none
! Compute U.S^-1/2 canonical orthogonalization
integer, intent(in) :: lda, ldc, n, m
double precision, intent(in) :: overlap(lda,n)
double precision, intent(inout) :: C(ldc,n)
double precision :: U(ldc,n)
double precision :: Vt(lda,n)
double precision :: D(n)
double precision :: S_half(lda,n)
double precision,allocatable :: work(:)
!DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D, work
integer :: info, lwork, i, j, k
double precision,allocatable :: overlap_tmp(:,:)
allocate (overlap_tmp(lda,n))
overlap_tmp = overlap
allocate(work(1))
lwork = -1
call dgesvd('A','A', n, n, overlap_tmp, lda, &
D, U, ldc, Vt, lda, work, lwork, info)
lwork = work(1)
deallocate(work)
allocate(work(lwork))
call dgesvd('A','A', n, n, overlap_tmp, lda, &
D, U, ldc, Vt, lda, work, lwork, info)
deallocate(work,overlap_tmp)
if (info /= 0) then
print *, info, ': SVD failed'
stop
endif
do i=1,n
if ( D(i) < 1.d-6 ) then
D(i) = 0.d0
else
D(i) = 1.d0/dsqrt(D(i))
endif
enddo
S_half = 0.d0
do k=1,n
do j=1,n
do i=1,n
S_half(i,j) += U(i,k)*D(k)*Vt(k,j)
enddo
enddo
enddo
do j=1,n
do i=1,m
U(i,j) = C(i,j)
enddo
enddo
C = 0.d0
do j=1,n
do i=1,m
do k=1,n
C(i,j) += U(i,k)*S_half(k,j)
enddo
enddo
enddo
end
subroutine get_pseudo_inverse(A,m,n,C,LDA)
! Find C = A^-1
implicit none
integer, intent(in) :: m,n, LDA
double precision, intent(in) :: A(LDA,n)
double precision, intent(out) :: C(n,m)
double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:)
integer :: info, lwork
integer :: i,j,k
allocate (D(n),U(m,n),Vt(n,n),work(1),A_tmp(m,n))
do j=1,n
do i=1,m
A_tmp(i,j) = A(i,j)
enddo
enddo
lwork = -1
call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info)
if (info /= 0) then
print *, info, ': SVD failed'
stop
endif
lwork = work(1)
deallocate(work)
allocate(work(lwork))
call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info)
if (info /= 0) then
print *, info, ': SVD failed'
stop
endif
do i=1,n
if (abs(D(i)) > 1.d-6) then
D(i) = 1.d0/D(i)
else
D(i) = 0.d0
endif
enddo
C = 0.d0
do i=1,m
do j=1,n
do k=1,n
C(j,i) += U(i,k) * D(k) * Vt(k,j)
enddo
enddo
enddo
deallocate(U,D,Vt,work,A_tmp)
end
subroutine find_rotation(A,LDA,B,m,C,n)
! Find A.C = B
implicit none
integer, intent(in) :: m,n, LDA
double precision, intent(in) :: A(LDA,n), B(LDA,n)
double precision, intent(out) :: C(n,n)
double precision, allocatable :: A_inv(:,:)
allocate(A_inv(LDA,n))
call get_pseudo_inverse(A,m,n,A_inv,LDA)
integer :: i,j,k
call dgemm('N','N',n,n,m,1.d0,A_inv,n,B,LDA,0.d0,C,n)
deallocate(A_inv)
end
subroutine apply_rotation(A,LDA,R,LDR,B,LDB,m,n)
implicit none
double precision, intent(in) :: R(LDR,n)
double precision, intent(in) :: A(LDA,n)
double precision, intent(out) :: B(LDB,n)
integer, intent(in) :: m,n, LDA, LDB, LDR
call dgemm('N','N',m,n,n,1.d0,A,LDA,R,LDR,0.d0,B,LDB)
end
subroutine jacobi_lapack(eigvalues,eigvectors,H,nmax,n)
implicit none
integer, intent(in) :: n,nmax
double precision, intent(out) :: eigvectors(nmax,n)
double precision, intent(out) :: eigvalues(n)
double precision, intent(in) :: H(nmax,n)
double precision,allocatable :: eigenvalues(:)
double precision,allocatable :: work(:)
double precision,allocatable :: A(:,:)
!eigvectors(i,j) = <d_i|psi_j> where d_i is the basis function and psi_j is the j th eigenvector
print*,nmax,n
allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax))
integer :: LWORK, info, i,j,l,k
double precision :: cpu_2, cpu_1
A=H
call cpu_time (cpu_1)
LWORK = 4*nmax
call dsyev( 'V', &
'U', &
n, &
A, &
nmax, &
eigenvalues, &
work, &
LWORK, &
info )
if (info < 0) then
print *, irp_here, ': the ',-info,'-th argument had an illegal value'
stop
else if (info > 0) then
print *, irp_here, ': the algorithm failed to converge; ',info,' off-diagonal'
print *, 'elements of an intermediate tridiagonal form did not converge to zero.'
stop
endif
call cpu_time (cpu_2)
eigvectors = 0.d0
eigvalues = 0.d0
do j = 1, n
eigvalues(j) = eigenvalues(j)
do i = 1, n
eigvectors(i,j) = A(i,j)
enddo
enddo
end

579
src/Utils/integration.irp.f Normal file
View File

@ -0,0 +1,579 @@
subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim)
! subroutine that transform the product of
! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
! into
! fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2)
implicit none
integer, intent(in) :: dim
integer, intent(in) :: a,b ! powers : (x-xa)**a_x = (x-A(1))**a(1)
double precision, intent(in) :: alpha, beta ! exponents
double precision, intent(in) :: A_center ! A center
double precision, intent(in) :: B_center ! B center
double precision, intent(out) :: P_center ! new center
double precision, intent(out) :: p ! new exponent
double precision, intent(out) :: fact_k ! constant factor
include 'constants.F'
double precision, intent(out) :: P_new(0:max_dim) ! polynom
integer, intent(out) :: iorder ! order of the polynoms
double precision :: P_a(0:max_dim), P_b(0:max_dim)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b
integer :: n_new,i,j
P_new = 0.d0
! you do the gaussiant product to get the new center and the new exponent
double precision :: p_inv,ab,d_AB
p = alpha+beta
p_inv = 1.d0/p
ab = alpha * beta
d_AB = (A_center - B_center) * (A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
fact_k = exp(-ab*p_inv * d_AB)
! you recenter the polynomw P_a and P_b on x
!call recentered_poly(P_a(0),A_center,P_center,a,dim)
!call recentered_poly(P_b(0),B_center,P_center,b,dim)
!DIR$ FORCEINLINE
call recentered_poly2(P_a(0),A_center,P_center,a,P_b(0),B_center,P_center,b)
n_new = 0
!DEC$ FORCEINLINE
call multiply_poly(P_a(0),a,P_b(0),b,P_new(0),n_new)
iorder = a + b
end
subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim)
! subroutine that transform the product of
! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
! into
! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
implicit none
include 'constants.F'
integer, intent(in) :: dim
integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
double precision, intent(in) :: alpha, beta ! exponents
double precision, intent(in) :: A_center(3) ! A center
double precision, intent(in) :: B_center (3) ! B center
double precision, intent(out) :: P_center(3) ! new center
double precision, intent(out) :: p ! new exponent
double precision, intent(out) :: fact_k ! constant factor
double precision, intent(out) :: P_new(0:max_dim,3) ! polynom
integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynoms
double precision :: P_a(0:max_dim,3), P_b(0:max_dim,3)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b
integer :: n_new,i,j
iorder(1) = 0
iorder(2) = 0
iorder(3) = 0
P_new(0,1) = 0.d0
P_new(0,2) = 0.d0
P_new(0,3) = 0.d0
!DEC$ FORCEINLINE
call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
if (fact_k < 1.d-8) then
fact_k = 0.d0
return
endif
!DEC$ FORCEINLINE
call recentered_poly2(P_a(0,1),A_center(1),P_center(1),a(1),P_b(0,1),B_center(1),P_center(1),b(1))
iorder(1) = a(1) + b(1)
!DIR$ VECTOR ALIGNED
do i=0,iorder(1)
P_new(i,1) = 0.d0
enddo
n_new=0
!DEC$ FORCEINLINE
call multiply_poly(P_a(0,1),a(1),P_b(0,1),b(1),P_new(0,1),n_new)
!DEC$ FORCEINLINE
call recentered_poly2(P_a(0,2),A_center(2),P_center(2),a(2),P_b(0,2),B_center(2),P_center(2),b(2))
iorder(2) = a(2) + b(2)
!DIR$ VECTOR ALIGNED
do i=0,iorder(2)
P_new(i,2) = 0.d0
enddo
n_new=0
!DEC$ FORCEINLINE
call multiply_poly(P_a(0,2),a(2),P_b(0,2),b(2),P_new(0,2),n_new)
!DEC$ FORCEINLINE
call recentered_poly2(P_a(0,3),A_center(3),P_center(3),a(3),P_b(0,3),B_center(3),P_center(3),b(3))
iorder(3) = a(3) + b(3)
!DIR$ VECTOR ALIGNED
do i=0,iorder(3)
P_new(i,3) = 0.d0
enddo
n_new=0
!DEC$ FORCEINLINE
call multiply_poly(P_a(0,3),a(3),P_b(0,3),b(3),P_new(0,3),n_new)
end
subroutine gaussian_product(a,xa,b,xb,k,p,xp)
implicit none
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
double precision , intent(in) :: a,b ! Exponents
double precision , intent(in) :: xa(3),xb(3) ! Centers
double precision , intent(out) :: p ! New exponent
double precision , intent(out) :: xp(3) ! New center
double precision , intent(out) :: k ! Constant
double precision :: p_inv
ASSERT (a>0.)
ASSERT (b>0.)
double precision :: xab(3), ab
!DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab
p = a+b
p_inv = 1.d0/(a+b)
ab = a*b
xab(1) = xa(1)-xb(1)
xab(2) = xa(2)-xb(2)
xab(3) = xa(3)-xb(3)
ab = ab*p_inv
k = ab*(xab(1)*xab(1)+xab(2)*xab(2)+xab(3)*xab(3))
if (k > 20.d0) then
k=0.d0
return
endif
k = exp(-k)
xp(1) = (a*xa(1)+b*xb(1))*p_inv
xp(2) = (a*xa(2)+b*xb(2))*p_inv
xp(3) = (a*xa(3)+b*xb(3))*p_inv
end subroutine
subroutine gaussian_product_x(a,xa,b,xb,k,p,xp)
implicit none
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
double precision , intent(in) :: a,b ! Exponents
double precision , intent(in) :: xa,xb ! Centers
double precision , intent(out) :: p ! New exponent
double precision , intent(out) :: xp ! New center
double precision , intent(out) :: k ! Constant
double precision :: p_inv
ASSERT (a>0.)
ASSERT (b>0.)
double precision :: xab, ab
p = a+b
p_inv = 1.d0/(a+b)
ab = a*b
xab = xa-xb
ab = ab*p_inv
k = ab*xab*xab
if (k > 20.d0) then
k=0.d0
return
endif
k = exp(-k)
xp = (a*xa+b*xb)*p_inv
end subroutine
subroutine multiply_poly(b,nb,c,nc,d,nd)
implicit none
! D(t) += B(t)*C(t)
! 4251722 + 3928617
! 4481076
! 185461
! 418740
integer, intent(in) :: nb, nc
integer, intent(out) :: nd
double precision, intent(in) :: b(0:nb), c(0:nc)
double precision, intent(inout) :: d(0:nb+nc)
integer :: ndtmp
integer :: ib, ic, id, k
if(ior(nc,nb) >= 0) then ! True if nc>=0 and nb>=0
continue
else
return
endif
ndtmp = nb+nc
!DIR$ VECTOR ALIGNED
do ic = 0,nc
d(ic) = d(ic) + c(ic) * b(0)
enddo
do ib=1,nb
d(ib) = d(ib) + c(0) * b(ib)
do ic = 1,nc
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
enddo
enddo
do nd = ndtmp,0,-1
if (d(nd) == 0.d0) then
cycle
endif
exit
enddo
end
subroutine add_poly(b,nb,c,nc,d,nd)
implicit none
! D(t) += B(t)+C(t)
integer, intent(inout) :: nb, nc
integer, intent(out) :: nd
double precision, intent(in) :: b(0:nb), c(0:nc)
double precision, intent(out) :: d(0:nb+nc)
nd = nb+nc
integer :: ib, ic, id
do ib=0,max(nb,nc)
d(ib) = d(ib) + c(ib) + b(ib)
enddo
do while ( (d(nd) == 0.d0).and.(nd>=0) )
nd -= 1
if (nd < 0) then
exit
endif
enddo
end
subroutine add_poly_multiply(b,nb,cst,d,nd)
implicit none
! D(t) += cst * B(t)
integer, intent(in) :: nb
integer, intent(inout) :: nd
double precision, intent(in) :: b(0:nb),cst
double precision, intent(inout) :: d(0:max(nb,nd))
nd = max(nd,nb)
if (nd /= -1) then
integer :: ib, ic, id
do ib=0,nb
d(ib) = d(ib) + cst*b(ib)
enddo
do while ( d(nd) == 0.d0 )
nd -= 1
if (nd < 0) then
exit
endif
enddo
endif
end
subroutine recentered_poly2(P_new,x_A,x_P,a,P_new2,x_B,x_Q,b)
! you enter with (x-x_A)^a
! you leave with sum_i=0,a c_i * (x-x_P)^i ==== P_new(i) * (x-x_P)^i
implicit none
integer, intent(in) :: a,b
double precision, intent(in) :: x_A,x_P,x_B,x_Q
double precision, intent(out) :: P_new(0:a),P_new2(0:b)
double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
double precision :: binom_func
integer :: i,j,k,l, minab, maxab
if ((a<0).or.(b<0) ) return
maxab = max(a,b)
minab = max(min(a,b),0)
pows_a(0) = 1.d0
pows_a(1) = (x_P - x_A)
pows_b(0) = 1.d0
pows_b(1) = (x_Q - x_B)
do i = 2,maxab
pows_a(i) = pows_a(i-1)*pows_a(1)
pows_b(i) = pows_b(i-1)*pows_b(1)
enddo
P_new (0) = pows_a(a)
P_new2(0) = pows_b(b)
do i = 1,min(minab,20)
P_new (i) = binom_transp(a-i,a) * pows_a(a-i)
P_new2(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = minab+1,min(a,20)
P_new (i) = binom_transp(a-i,a) * pows_a(a-i)
enddo
do i = minab+1,min(b,20)
P_new2(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = 101,a
P_new(i) = binom_func(a,a-i) * pows_a(a-i)
enddo
do i = 101,b
P_new2(i) = binom_func(b,b-i) * pows_b(b-i)
enddo
end
double precision function F_integral(n,p)
! function that calculates the following integral
! sum (x) between [-infty;+infty] of x^n exp(-p*x^2)
implicit none
integer :: n
double precision :: p
integer :: i,j
double precision :: accu,sqrt_p,fact_ratio,tmp,fact
include 'include/constants.F'
if(n < 0)then
F_integral = 0.d0
endif
if(iand(n,1).ne.0)then
F_integral = 0.d0
return
endif
sqrt_p = 1.d0/dsqrt(p)
if(n==0)then
F_integral = sqpi * sqrt_p
return
endif
F_integral = sqpi * 0.5d0**n * sqrt_p**(n+1) * fact(n)/fact(ishft(n,-1))
end
double precision function rint(n,rho)
implicit none
double precision :: rho,u,rint1,v,val0,rint_large_n,u_inv
integer :: n,k
double precision, parameter :: pi=3.14159265359d0
double precision, parameter :: dsqpi=1.77245385091d0
double precision :: two_rho_inv
! double precision :: rint_brut
! rint = rint_brut(n,rho,10000)
! return
if(n.eq.0)then
if(rho == 0.d0)then
rint=1.d0
else
u_inv=1.d0/dsqrt(rho)
u=rho*u_inv
rint=0.5d0*u_inv*dsqpi*erf(u)
endif
return
endif
if(rho.lt.1.d0)then
rint=rint1(n,rho)
else
if(n.le.20)then
u_inv=1.d0/dsqrt(rho)
v=dexp(-rho)
u=rho*u_inv
two_rho_inv = 0.5d0*u_inv*u_inv
val0=0.5d0*u_inv*dsqpi*erf(u)
rint=(val0-v)*two_rho_inv
do k=2,n
rint=(rint*dfloat(k+k-1)-v)*two_rho_inv
enddo
else
rint=rint_large_n(n,rho)
endif
endif
end
double precision function rint_sum(n_pt_out,rho,d1)
implicit none
integer, intent(in) :: n_pt_out
double precision, intent(in) :: rho,d1(0:n_pt_out)
double precision :: u,rint1,v,val0,rint_large_n,u_inv
integer :: n,k,i
double precision, parameter :: pi=3.14159265359d0
double precision, parameter :: dsqpi=1.77245385091d0
double precision :: two_rho_inv, rint_tmp, di
if(rho < 1.d0)then
if(rho == 0.d0)then
rint_sum=d1(0)
else
u_inv=1.d0/dsqrt(rho)
u=rho*u_inv
rint_sum=0.5d0*u_inv*dsqpi*erf(u) *d1(0)
endif
do i=2,n_pt_out,2
n = ishft(i,-1)
rint_sum = rint_sum + d1(i)*rint1(n,rho)
enddo
else
v=dexp(-rho)
u_inv=1.d0/dsqrt(rho)
u=rho*u_inv
two_rho_inv = 0.5d0*u_inv*u_inv
val0=0.5d0*u_inv*dsqpi*erf(u)
rint_sum=val0*d1(0)
rint_tmp=(val0-v)*two_rho_inv
di = 3.d0
do i=2,min(n_pt_out,40),2
rint_sum = rint_sum + d1(i)*rint_tmp
rint_tmp = (rint_tmp*di-v)*two_rho_inv
di = di+2.d0
enddo
do i=42,n_pt_out,2
n = ishft(i,-1)
rint_sum = rint_sum + d1(i)*rint_large_n(n,rho)
enddo
endif
end
double precision function hermite(n,x)
implicit none
integer :: n,k
double precision :: h0,x,h1,h2
h0=1.d0
if(n.eq.0)then
hermite=h0
return
endif
h1=x+x
if(n.eq.1)then
hermite=h1
return
endif
do k=1,n-1
h2=(x+x)*h1-dfloat(k+k)*h0
h0=h1
h1=h2
enddo
hermite=h2
end
double precision function rint_large_n(n,rho)
implicit none
integer :: n,k,l
double precision :: rho,u,accu,eps,t1,t2,fact,alpha_k,rajout,hermite
u=dsqrt(rho)
accu=0.d0
k=0
eps=1.d0
do while (eps.gt.1.d-15)
t1=1.d0
do l=0,k
t1=t1*(n+n+l+1.d0)
enddo
t2=0.d0
do l=0,k
t2=t2+(-1.d0)**l/fact(l+1)/fact(k-l)
enddo
alpha_k=t2*fact(k+1)*fact(k)*(-1.d0)**k
alpha_k= alpha_k/t1
rajout=(-1.d0)**k*u**k*hermite(k,u)*alpha_k/fact(k)
accu=accu+rajout
eps=dabs(rajout)/accu
k=k+1
enddo
rint_large_n=dexp(-rho)*accu
end
double precision function rint_brut(n,rho,npts)
implicit double precision(a-h,o-z)
double precision :: fi(4), t2(4), accu(4), dt(4), rho4(4), four(4)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fi, t2, accu, dt, rho4, four
accu(1:4)=0.d0
dt(1)=1.d0/dfloat(npts)
dt(2)=dt(1)
dt(3)=dt(1)
dt(4)=dt(1)
rho4(1:4) = (/ rho, rho, rho, rho /)
fi(1:4)= (/ 0.5d0, 1.5d0, 2.5d0, 3.5d0 /)
four(1:4) = (/ 4.d0, 4.d0, 4.d0, 4.d0 /)
select case(n)
case (1)
do i=1,npts,4
!DIR$ VECTOR ALIGNED
t2(1:4)=fi(1:4)*dt(1:4)
!DIR$ VECTOR ALIGNED
t2(1:4) = t2(1:4)*t2(1:4)
!DIR$ VECTOR ALIGNED
accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*t2(1:4)
!DIR$ VECTOR ALIGNED
fi(1:4) = fi(1:4)+four(1:4)
enddo
case (2)
do i=1,npts,4
!DIR$ VECTOR ALIGNED
t2(1:4)=fi(1:4)*dt(1:4)
!DIR$ VECTOR ALIGNED
t2(1:4) = t2(1:4)*t2(1:4)
!DIR$ VECTOR ALIGNED
accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*t2(1:4)*t2(1:4)
!DIR$ VECTOR ALIGNED
fi(1:4) = fi(1:4)+four(1:4)
enddo
case (3)
do i=1,npts,4
!DIR$ VECTOR ALIGNED
t2(1:4)=fi(1:4)*dt(1:4)
!DIR$ VECTOR ALIGNED
t2(1:4) = t2(1:4)*t2(1:4)
!DIR$ VECTOR ALIGNED
accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*t2(1:4)*t2(1:4)*t2(1:4)
!DIR$ VECTOR ALIGNED
fi(1:4) = fi(1:4)+four(1:4)
enddo
case default
do i=1,npts,4
!DIR$ VECTOR ALIGNED
t2(1:4)=fi(1:4)*dt(1:4)
!DIR$ VECTOR ALIGNED
t2(1:4) = t2(1:4)*t2(1:4)
!DIR$ VECTOR ALIGNED
accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*(t2(1:4)**(n))
!DIR$ VECTOR ALIGNED
fi(1:4) = fi(1:4)+four(1:4)
enddo
end select
rint_brut=sum(accu)*dt(1)
end
double precision function rint1(n,rho)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: rho
double precision, parameter :: eps=1.d-15
double precision :: rho_tmp, diff
integer :: k
rint1=inv_int(n+n+1)
rho_tmp = 1.d0
do k=1,20
rho_tmp = -rho_tmp*rho
diff=rho_tmp*fact_inv(k)*inv_int(ishft(k+n,1)+1)
rint1=rint1+diff
if (dabs(diff) > eps) then
cycle
endif
return
enddo
write(*,*)'pb in rint1 k too large!'
stop
end

135
src/Utils/one_e_integration.irp.f Normal file
View File

@ -0,0 +1,135 @@
double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_A,power_B,dim)
implicit none
! calculates the following overlap :
! sum (x) between [-infty;+infty] of (x-A_x)^ax (x-B_x)^bx exp(-alpha(x-A_x)^2) exp(-beta(x-B_X)^2)
include 'include/constants.F'
integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynoms
double precision,intent(in) :: A_center,B_center ! center of the x1 functions
integer,intent(in) :: power_A, power_B ! power of the x1 functions
double precision :: P_new(0:max_dim),P_center,fact_p,p,alpha,beta
integer :: iorder_p
call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
if(fact_p.lt.0.000001d0)then
overlap_gaussian_x = 0.d0
return
endif
overlap_gaussian_x = 0.d0
integer :: i
double precision :: F_integral
do i = 0,iorder_p
overlap_gaussian_x += P_new(i) * F_integral(i,p)
enddo
overlap_gaussian_x*= fact_p
end
subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim)
implicit none
! .. math::
!
! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx \\
! S = S_x S_y S_z
!
!
include 'include/constants.F'
integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynoms
double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions
double precision, intent(in) :: alpha,beta
integer,intent(in) :: power_A(3), power_B(3) ! power of the x1 functions
double precision, intent(out) :: overlap_x,overlap_y,overlap_z,overlap
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p
double precision :: F_integral_tab(0:max_dim)
integer :: iorder_p(3)
call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
if(fact_p.lt.0.000001d0)then
overlap_x = 0.d0
overlap_y = 0.d0
overlap_z = 0.d0
overlap = 0.d0
return
endif
integer :: nmax
double precision :: F_integral
nmax = maxval(iorder_p)
do i = 0,nmax
F_integral_tab(i) = F_integral(i,p)
enddo
overlap_x = P_new(0,1) * F_integral_tab(0)
overlap_y = P_new(0,2) * F_integral_tab(0)
overlap_z = P_new(0,3) * F_integral_tab(0)
integer :: i
do i = 1,iorder_p(1)
overlap_x += P_new(i,1) * F_integral_tab(i)
enddo
call gaussian_product_x(alpha,A_center(1),beta,B_center(1),fact_p,p,P_center(1))
overlap_x *= fact_p
do i = 1,iorder_p(2)
overlap_y += P_new(i,2) * F_integral_tab(i)
enddo
call gaussian_product_x(alpha,A_center(2),beta,B_center(2),fact_p,p,P_center(2))
overlap_y *= fact_p
do i = 1,iorder_p(3)
overlap_z += P_new(i,3) * F_integral_tab(i)
enddo
call gaussian_product_x(alpha,A_center(3),beta,B_center(3),fact_p,p,P_center(3))
overlap_z *= fact_p
overlap = overlap_x * overlap_y * overlap_z
end
subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
implicit none
! compute the following integral :
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ]
integer :: i,j,k,l
integer,intent(in) :: power_A,power_B
double precision, intent(in) :: lower_exp_val
double precision,intent(in) :: A_center, B_center,alpha,beta
double precision, intent(out) :: overlap_x,dx
integer, intent(in) :: nx
double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho
double precision :: P_center
if(power_A.lt.0.or.power_B.lt.0)then
overlap_x = 0.d0
dx = 0.d0
return
endif
p = alpha + beta
p_inv= 1.d0/p
rho = alpha * beta * p_inv
dist = (A_center - B_center)*(A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
factor = dexp(-rho * dist)
double precision :: tmp
tmp = dsqrt(lower_exp_val/p)
x_min = P_center - tmp
x_max = P_center + tmp
domain = x_max-x_min
dx = domain/dble(nx)
overlap_x = 0.d0
x = x_min
do i = 1, nx
x += dx
overlap_x += abs((x-A_center)**power_A * (x-B_center)**power_B) * dexp(-p * (x-P_center)*(x-P_center))
enddo
overlap_x *= factor * dx
end

408
src/Utils/sort.irp.f Normal file
View File

@ -0,0 +1,408 @@
BEGIN_TEMPLATE
subroutine insertion_$Xsort (x,iorder,isize)
implicit none
$type,intent(inout) :: x(isize)
integer,intent(inout) :: iorder(isize)
integer,intent(in) :: isize
$type :: xtmp
integer :: i, i0, j, jmax
do i=1,isize
xtmp = x(i)
i0 = iorder(i)
j = i-1
do j=i-1,1,-1
if ( x(j) > xtmp ) then
x(j+1) = x(j)
iorder(j+1) = iorder(j)
else
exit
endif
enddo
x(j+1) = xtmp
iorder(j+1) = i0
enddo
end subroutine insertion_$Xsort
subroutine heap_$Xsort(x,iorder,isize)
implicit none
$type,intent(inout) :: x(isize)
integer,intent(inout) :: iorder(isize)
integer,intent(in) :: isize
integer :: i, k, j, l, i0
$type :: xtemp
l = isize/2+1
k = isize
do while (.True.)
if (l>1) then
l=l-1
xtemp = x(l)
i0 = iorder(l)
else
xtemp = x(k)
i0 = iorder(k)
x(k) = x(1)
iorder(k) = iorder(1)
k = k-1
if (k == 1) then
x(1) = xtemp
iorder(1) = i0
exit
endif
endif
i=l
j = ishft(l,1)
do while (j<k)
if ( x(j) < x(j+1) ) then
j=j+1
endif
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = ishft(j,1)
else
j = k+1
endif
enddo
if (j==k) then
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = ishft(j,1)
else
j = k+1
endif
endif
x(i) = xtemp
iorder(i) = i0
enddo
end subroutine heap_$Xsort
subroutine heap_$Xsort_big(x,iorder,isize)
implicit none
$type,intent(inout) :: x(isize)
integer*8,intent(inout) :: iorder(isize)
integer*8,intent(in) :: isize
integer*8 :: i, k, j, l, i0
$type :: xtemp
l = isize/2+1
k = isize
do while (.True.)
if (l>1) then
l=l-1
xtemp = x(l)
i0 = iorder(l)
else
xtemp = x(k)
i0 = iorder(k)
x(k) = x(1)
iorder(k) = iorder(1)
k = k-1
if (k == 1) then
x(1) = xtemp
iorder(1) = i0
exit
endif
endif
i=l
j = ishft(l,1)
do while (j<k)
if ( x(j) < x(j+1) ) then
j=j+1
endif
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = ishft(j,1)
else
j = k+1
endif
enddo
if (j==k) then
if (xtemp < x(j)) then
x(i) = x(j)
iorder(i) = iorder(j)
i = j
j = ishft(j,1)
else
j = k+1
endif
endif
x(i) = xtemp
iorder(i) = i0
enddo
end subroutine heap_$Xsort$big
subroutine $Xsort(x,iorder,isize)
implicit none
$type,intent(inout) :: x(isize)
integer,intent(inout) :: iorder(isize)
integer,intent(in) :: isize
if (isize < 32) then
call insertion_$Xsort(x,iorder,isize)
else
call heap_$Xsort(x,iorder,isize)
endif
end subroutine $Xsort
SUBST [ X, type ]
; real ;;
d ; double precision ;;
i ; integer ;;
i8 ; integer*8 ;;
i2 ; integer*2 ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine $Xset_order(x,iorder,isize)
implicit none
integer :: isize
$type :: x(*)
$type,allocatable :: xtmp(:)
integer :: iorder(*)
integer :: i
allocate(xtmp(isize))
do i=1,isize
xtmp(i) = x(iorder(i))
enddo
do i=1,isize
x(i) = xtmp(i)
enddo
deallocate(xtmp)
end
SUBST [ X, type ]
; real ;;
d ; double precision ;;
i ; integer ;;
i8; integer*8 ;;
i2; integer*2 ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine insertion_$Xsort_big (x,iorder,isize)
implicit none
$type,intent(inout) :: x(isize)
integer*8,intent(inout) :: iorder(isize)
integer*8,intent(in) :: isize
$type :: xtmp
integer*8 :: i, i0, j, jmax
do i=1_8,isize
xtmp = x(i)
i0 = iorder(i)
j = i-1_8
do j=i-1_8,1_8,-1_8
if ( x(j) > xtmp ) then
x(j+1_8) = x(j)
iorder(j+1_8) = iorder(j)
else
exit
endif
enddo
x(j+1_8) = xtmp
iorder(j+1_8) = i0
enddo
end subroutine insertion_$Xsort
subroutine $Xset_order_big(x,iorder,isize)
implicit none
integer*8 :: isize
$type :: x(*)
$type, allocatable :: xtmp(:)
integer*8 :: iorder(*)
integer*8 :: i
allocate(xtmp(isize))
do i=1_8,isize
xtmp(i) = x(iorder(i))
enddo
do i=1_8,isize
x(i) = xtmp(i)
enddo
deallocate(xtmp)
end
SUBST [ X, type ]
; real ;;
d ; double precision ;;
i ; integer ;;
i8; integer*8 ;;
i2; integer*2 ;;
END_TEMPLATE
BEGIN_TEMPLATE
recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix)
implicit none
$int_type, intent(in) :: isize
$int_type, intent(inout) :: iorder(isize)
$type, intent(inout) :: x(isize)
integer, intent(in) :: iradix
integer :: iradix_new
$type, allocatable :: x2(:), x1(:)
$int_type, allocatable :: iorder1(:),iorder2(:)
$int_type :: i0, i1, i2, i3, i
integer, parameter :: integer_size=$octets
$type, parameter :: zero=$zero
$type :: mask
integer :: nthreads, omp_get_num_threads
!DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1
if (iradix == -1) then
! Find most significant bit
i0 = 0_8
i3 = -1_8
do i=1,isize
i3 = max(i3,x(i))
enddo
iradix_new = integer_size-1-leadz(i3)
mask = ibset(zero,iradix_new)
nthreads = 1
! nthreads = 1+ishft(omp_get_num_threads(),-1)
integer :: err
allocate(x1(isize/nthreads+1),iorder1(isize/nthreads+1),x2(isize/nthreads+1),iorder2(isize/nthreads+1),stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to allocate arrays'
stop
endif
i1=1_8
i2=1_8
do i=1,isize
if (iand(mask,x(i)) == zero) then
iorder1(i1) = iorder(i)
x1(i1) = x(i)
i1 = i1+1_8
else
iorder2(i2) = iorder(i)
x2(i2) = x(i)
i2 = i2+1_8
endif
enddo
i1=i1-1_8
i2=i2-1_8
do i=1,i1
iorder(i0+i) = iorder1(i)
x(i0+i) = x1(i)
enddo
i0 = i0+i1
i3 = i0
deallocate(x1,iorder1,stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to deallocate arrays x1, iorder1'
stop
endif
do i=1,i2
iorder(i0+i) = iorder2(i)
x(i0+i) = x2(i)
enddo
i0 = i0+i2
deallocate(x2,iorder2,stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to deallocate arrays x2, iorder2'
stop
endif
if (i3>1) then
call $Xradix_sort$big(x,iorder,i3,iradix_new-1)
endif
if (isize-i3>1) then
call $Xradix_sort$big(x(i3+1),iorder(i3+1),isize-i3,iradix_new-1)
endif
return
endif
ASSERT (iradix > 0)
if (isize < 48) then
call insertion_$Xsort$big(x,iorder,isize)
return
endif
allocate(x2(isize),iorder2(isize),stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to allocate arrays x1, iorder1'
stop
endif
mask = ibset(zero,iradix)
i0=1
i1=1
do i=1,isize
if (iand(mask,x(i)) == zero) then
iorder(i0) = iorder(i)
x(i0) = x(i)
i0 = i0+1
else
iorder2(i1) = iorder(i)
x2(i1) = x(i)
i1 = i1+1
endif
enddo
i0=i0-1
i1=i1-1
do i=1,i1
iorder(i0+i) = iorder2(i)
x(i0+i) = x2(i)
enddo
deallocate(x2,iorder2,stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to allocate arrays x2, iorder2'
stop
endif
if (iradix == 0) then
return
endif
if (i1>1) then
call $Xradix_sort$big(x(i0+1),iorder(i0+1),i1,iradix-1)
endif
if (i0>1) then
call $Xradix_sort$big(x,iorder,i0,iradix-1)
endif
end
SUBST [ X, type, octets, is_big, big, int_type, zero ]
i ; integer ; 32 ; .False. ; ; integer ; 0;;
i8 ; integer*8 ; 32 ; .False. ; ; integer ; 0_8;;
i2 ; integer*2 ; 32 ; .False. ; ; integer ; 0;;
i ; integer ; 64 ; .True. ; _big ; integer*8 ; 0 ;;
i8 ; integer*8 ; 64 ; .True. ; _big ; integer*8 ; 0_8 ;;
END_TEMPLATE

178
src/Utils/util.irp.f Normal file
View File

@ -0,0 +1,178 @@
double precision function binom_func(i,j)
implicit none
integer,intent(in) :: i,j
double precision :: fact, f
integer, save :: ifirst
double precision, save :: memo(0:15,0:15)
!DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: memo
integer :: k,l
if (ifirst == 0) then
ifirst = 1
do k=0,15
f = fact(k)
do l=0,15
memo(k,l) = f/(fact(l)*fact(k-l))
enddo
enddo
endif
if ( (i<=15).and.(j<=15) ) then
binom_func = memo(i,j)
else
binom_func = fact(i)/(fact(j)*fact(i-j))
endif
end
BEGIN_PROVIDER [ double precision, binom, (0:20,0:20) ]
&BEGIN_PROVIDER [ double precision, binom_transp, (0:20,0:20) ]
implicit none
BEGIN_DOC
! Binomial coefficients
END_DOC
integer :: k,l
double precision :: fact, f
do k=0,20
f = fact(k)
do l=0,20
binom(k,l) = f/(fact(l)*fact(k-l))
binom_transp(l,k) = binom(k,l)
enddo
enddo
END_PROVIDER
integer function align_double(n)
implicit none
integer :: n
include 'include/constants.F'
if (mod(n,SIMD_vector/4) /= 0) then
align_double= n + SIMD_vector/4 - mod(n,SIMD_vector/4)
else
align_double= n
endif
end
double precision function fact(n)
implicit none
integer :: n
double precision, save :: memo(1:100)
integer, save :: memomax = 1
if (n<=memomax) then
if (n<2) then
fact = 1.d0
else
fact = memo(n)
endif
return
endif
integer :: i
memo(1) = 1.d0
do i=memomax+1,min(n,100)
memo(i) = memo(i-1)*float(i)
enddo
memomax = min(n,100)
fact = memo(memomax)
do i=101,n
fact = fact*float(i)
enddo
end function
BEGIN_PROVIDER [ double precision, fact_inv, (128) ]
implicit none
BEGIN_DOC
! 1.d0/fact(k)
END_DOC
integer :: i
double precision :: fact
do i=1,size(fact_inv)
fact_inv(i) = 1.d0/fact(i)
enddo
END_PROVIDER
double precision function dble_fact(n) result(fact2)
implicit none
integer :: n
double precision, save :: memo(1:100)
integer, save :: memomax = 1
ASSERT (iand(n,1) /= 0)
if (n<=memomax) then
if (n<3) then
fact2 = 1.d0
else
fact2 = memo(n)
endif
return
endif
integer :: i
memo(1) = 1.d0
do i=memomax+2,min(n,99),2
memo(i) = memo(i-2)* float(i)
enddo
memomax = min(n,99)
fact2 = memo(memomax)
do i=101,n,2
fact2 = fact2*float(i)
enddo
end function
subroutine write_git_log(iunit)
implicit none
integer, intent(in) :: iunit
write(iunit,*) '----------------'
write(iunit,*) 'Last git commit:'
BEGIN_SHELL [ /bin/bash ]
git log -1 | sed "s/'//g"| sed "s/^/ write(iunit,*) '/g" | sed "s/$/'/g"
END_SHELL
write(iunit,*) '----------------'
end
BEGIN_PROVIDER [ double precision, inv_int, (128) ]
implicit none
BEGIN_DOC
! 1/i
END_DOC
integer :: i
do i=1,size(inv_int)
inv_int(i) = 1.d0/dble(i)
enddo
END_PROVIDER
subroutine wall_time(t)
implicit none
double precision, intent(out) :: t
integer :: c
integer, save :: rate = 0
if (rate == 0) then
CALL SYSTEM_CLOCK(count_rate=rate)
endif
CALL SYSTEM_CLOCK(count=c)
t = dble(c)/dble(rate)
end
BEGIN_PROVIDER [ integer, nproc ]
implicit none
BEGIN_DOC
! Number of current openmp threads
END_DOC
integer :: omp_get_num_threads
nproc = 1
!$OMP PARALLEL
!$OMP MASTER
!$ nproc = omp_get_num_threads()
!$OMP END MASTER
!$OMP END PARALLEL
END_PROVIDER

10
src/include/constants.F Normal file
View File

@ -0,0 +1,10 @@
integer, parameter :: max_dim = 255
integer, parameter :: SIMD_vector = 32
double precision, parameter :: pi = dacos(-1.d0)
double precision, parameter :: sqpi = dsqrt(dacos(-1.d0))
double precision, parameter :: pi_5_2 = 34.9868366552d0
double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0)
double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0)
double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0))
double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))

6
src/selected_ci.ezfio_config Normal file
View File

@ -0,0 +1,6 @@
work
empty logical
save
empty logical