From 6c29d22aa2f8b5b9208ae47bc6f9f97fee7c693a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Apr 2014 18:37:27 +0200 Subject: [PATCH] Added Utils, Electrons, Ezfio, Nuclei, include, AOs. --- README.md | 2 +- src/AOs/Makefile | 9 + src/AOs/aos.ezfio_config | 10 + src/AOs/aos.irp.f | 192 ++++++++ src/Electrons/electrons.ezfio_config | 5 + src/Electrons/electrons.irp.f | 21 + src/Ezfio_files/Makefile | 8 + src/Ezfio_files/ezfio.irp.f | 15 + src/Ezfio_files/get_unit_and_open.irp.f | 43 ++ src/Makefile | 21 + src/Makefile.common | 39 ++ src/Makefile.config.example | 29 ++ src/Nuclei/Makefile | 9 + src/Nuclei/nuclei.ezfio_config | 6 + src/Nuclei/nuclei.irp.f | 134 ++++++ src/Nuclei/test.irp.f | 3 + src/Utils/LinearAlgebra.irp.f | 195 ++++++++ src/Utils/integration.irp.f | 579 ++++++++++++++++++++++++ src/Utils/one_e_integration.irp.f | 135 ++++++ src/Utils/sort.irp.f | 408 +++++++++++++++++ src/Utils/util.irp.f | 178 ++++++++ src/include/constants.F | 10 + src/selected_ci.ezfio_config | 6 + 23 files changed, 2056 insertions(+), 1 deletion(-) create mode 100644 src/AOs/Makefile create mode 100644 src/AOs/aos.ezfio_config create mode 100644 src/AOs/aos.irp.f create mode 100644 src/Electrons/electrons.ezfio_config create mode 100644 src/Electrons/electrons.irp.f create mode 100644 src/Ezfio_files/Makefile create mode 100644 src/Ezfio_files/ezfio.irp.f create mode 100644 src/Ezfio_files/get_unit_and_open.irp.f create mode 100644 src/Makefile create mode 100644 src/Makefile.common create mode 100644 src/Makefile.config.example create mode 100644 src/Nuclei/Makefile create mode 100644 src/Nuclei/nuclei.ezfio_config create mode 100644 src/Nuclei/nuclei.irp.f create mode 100644 src/Nuclei/test.irp.f create mode 100644 src/Utils/LinearAlgebra.irp.f create mode 100644 src/Utils/integration.irp.f create mode 100644 src/Utils/one_e_integration.irp.f create mode 100644 src/Utils/sort.irp.f create mode 100644 src/Utils/util.irp.f create mode 100644 src/include/constants.F create mode 100644 src/selected_ci.ezfio_config diff --git a/README.md b/README.md index 8de94373..1192bd72 100644 --- a/README.md +++ b/README.md @@ -4,5 +4,5 @@ Quantum package Set of quantum chemistry programs and libraries. For more information, you can visit the -`wiki of the project `_ +wiki of the project : http://github.com/LCPQ/quantum_package/wiki diff --git a/src/AOs/Makefile b/src/AOs/Makefile new file mode 100644 index 00000000..9f436863 --- /dev/null +++ b/src/AOs/Makefile @@ -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) diff --git a/src/AOs/aos.ezfio_config b/src/AOs/aos.ezfio_config new file mode 100644 index 00000000..53ac8d85 --- /dev/null +++ b/src/AOs/aos.ezfio_config @@ -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) + + diff --git a/src/AOs/aos.irp.f b/src/AOs/aos.irp.f new file mode 100644 index 00000000..ae5b67c7 --- /dev/null +++ b/src/AOs/aos.irp.f @@ -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 + diff --git a/src/Electrons/electrons.ezfio_config b/src/Electrons/electrons.ezfio_config new file mode 100644 index 00000000..d89823b7 --- /dev/null +++ b/src/Electrons/electrons.ezfio_config @@ -0,0 +1,5 @@ +electrons + elec_alpha_num integer + elec_beta_num integer + elec_num integer = electrons_elec_alpha_num + electrons_elec_beta_num + diff --git a/src/Electrons/electrons.irp.f b/src/Electrons/electrons.irp.f new file mode 100644 index 00000000..9fab8d76 --- /dev/null +++ b/src/Electrons/electrons.irp.f @@ -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 + + + diff --git a/src/Ezfio_files/Makefile b/src/Ezfio_files/Makefile new file mode 100644 index 00000000..ac9b7c52 --- /dev/null +++ b/src/Ezfio_files/Makefile @@ -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) diff --git a/src/Ezfio_files/ezfio.irp.f b/src/Ezfio_files/ezfio.irp.f new file mode 100644 index 00000000..90c71a5b --- /dev/null +++ b/src/Ezfio_files/ezfio.irp.f @@ -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, ' ' + stop 1 + endif + call getarg(1,ezfio_filename) + call ezfio_set_file(ezfio_filename) +END_PROVIDER + diff --git a/src/Ezfio_files/get_unit_and_open.irp.f b/src/Ezfio_files/get_unit_and_open.irp.f new file mode 100644 index 00000000..259b2eec --- /dev/null +++ b/src/Ezfio_files/get_unit_and_open.irp.f @@ -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 + diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 00000000..f9c33910 --- /dev/null +++ b/src/Makefile @@ -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 + diff --git a/src/Makefile.common b/src/Makefile.common new file mode 100644 index 00000000..ee853f2e --- /dev/null +++ b/src/Makefile.common @@ -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)) + diff --git a/src/Makefile.config.example b/src/Makefile.config.example new file mode 100644 index 00000000..82e0ca27 --- /dev/null +++ b/src/Makefile.config.example @@ -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 diff --git a/src/Nuclei/Makefile b/src/Nuclei/Makefile new file mode 100644 index 00000000..0179932e --- /dev/null +++ b/src/Nuclei/Makefile @@ -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) diff --git a/src/Nuclei/nuclei.ezfio_config b/src/Nuclei/nuclei.ezfio_config new file mode 100644 index 00000000..07be4d29 --- /dev/null +++ b/src/Nuclei/nuclei.ezfio_config @@ -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) + diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f new file mode 100644 index 00000000..7cf4ce11 --- /dev/null +++ b/src/Nuclei/nuclei.irp.f @@ -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 diff --git a/src/Nuclei/test.irp.f b/src/Nuclei/test.irp.f new file mode 100644 index 00000000..46f28cce --- /dev/null +++ b/src/Nuclei/test.irp.f @@ -0,0 +1,3 @@ +program test + print *, nuclear_repulsion +end diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f new file mode 100644 index 00000000..b6023179 --- /dev/null +++ b/src/Utils/LinearAlgebra.irp.f @@ -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) = 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 diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f new file mode 100644 index 00000000..c3ebaf5b --- /dev/null +++ b/src/Utils/integration.irp.f @@ -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 diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f new file mode 100644 index 00000000..ab246ff5 --- /dev/null +++ b/src/Utils/one_e_integration.irp.f @@ -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 + + + + diff --git a/src/Utils/sort.irp.f b/src/Utils/sort.irp.f new file mode 100644 index 00000000..e821895c --- /dev/null +++ b/src/Utils/sort.irp.f @@ -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 (j1) 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 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 + + diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f new file mode 100644 index 00000000..90d4603f --- /dev/null +++ b/src/Utils/util.irp.f @@ -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 + diff --git a/src/include/constants.F b/src/include/constants.F new file mode 100644 index 00000000..632dc50b --- /dev/null +++ b/src/include/constants.F @@ -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)) diff --git a/src/selected_ci.ezfio_config b/src/selected_ci.ezfio_config new file mode 100644 index 00000000..8c64c790 --- /dev/null +++ b/src/selected_ci.ezfio_config @@ -0,0 +1,6 @@ +work + empty logical + +save + empty logical +