LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity

Added Mono Integrals. Works in openMP

Browse Source
This commit is contained in:
Anthony Scemama 2014-04-16 22:16:32 +02:00
parent e1ed65a1a6
commit 0764aa259d
34 changed files with 255561 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
scripts/create_test_ref.sh
View File

@ -40,7 +40,14 @@ printf "Running tests...."
for dir in ${QPACKAGE_ROOT}/data/inputs/*.ezfio
do
printf " '%s' : {\n " $(basename ${dir})
./${TEST_EXE} ${dir} | sed "s/\([^ ]*\) *\(:\) *\([^ ]*\)/'\1' : \3,/g"
OMP_NUM_THREADS=1 python << EOF
import subprocess
lines = subprocess.check_output("./${TEST_EXE} ${dir}", shell=True)
for line in lines.splitlines():
buffer = line.split(':')
print "'%s' : %s,"%(buffer[0].strip(), buffer[1].strip())
EOF
# sed "s/ \(.*\) *\(:\) *\([^ ]*\)/'\1' : \3,/g"
printf " },\n"
done >> ${REF_FILE}
printf "}\n" >> ${REF_FILE}

21
scripts/run_tests.py
View File

@ -1,16 +1,22 @@
#!/usr/bin/python
import sys,os
import subprocess
DEBUG=False
def run_test(test_name,inp):
command = './'+test_name+" ${QPACKAGE_ROOT}/data/inputs/"+inp
result = subprocess.check_output(command, shell=True)
try:
result = subprocess.check_output(command, shell=True)
except:
result = sys.exc_info()[0]
return result
if __name__ == '__main__':
import unittest
import subprocess
from math import *
from multiprocessing import Pool
@ -33,7 +39,7 @@ if __name__ == '__main__':
template = """
class $test(unittest.TestCase):
default_precision = 1.e-10
default_precision = 10
execfile('$test.ref')
@ -41,12 +47,16 @@ class $test(unittest.TestCase):
tasks = {}
def setUp(self):
if DEBUG: return
for d in self.data.keys():
if d not in self.tasks:
self.tasks[d] = pool.apply_async(run_test, [self.name, d])
def _test_input(self,inp):
output = self.tasks[inp].get()
if DEBUG:
output = run_test(self.name,inp)
else:
output = self.tasks[inp].get()
for line in output.splitlines():
buffer = line.split(':')
if len(buffer) == 1:
@ -59,7 +69,8 @@ class $test(unittest.TestCase):
precision = self.default_precision
if type(r) == float:
self.assertAlmostEqual(self.data[inp][l], r,
places=abs(int(log10(precision*max(abs(self.data[inp][l]),1.e-12)))), msg=None)
places=precision, msg=None)
#places=abs(int(log10(precision*max(abs(self.data[inp][l]),1.e-12)))), msg=None)
else:
self.assertEqual(self.data[inp][l], r, msg=None)

17
src/Hartree_Fock/README.rst Normal file
View File

@ -0,0 +1,17 @@
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_

2
src/MOs/utils.irp.f
View File

@ -44,5 +44,5 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label)
deallocate(mo_coef_new,R,eigvalues)
mo_label = label
SOFT_TOUCH mo_coef
SOFT_TOUCH mo_coef mo_label
end

2
src/Makefile.common
View File

@ -124,7 +124,7 @@ Makefile.depend: Makefile
include irpf90.make
test:
make -j 1 -C tests
DEBUG=1 make -j 1 -C tests
# Dummy rule to enable to force recompilation
FORCE:

0
src/MonoInts/ASSUMPTIONS.rst Normal file
View File

8
src/MonoInts/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/MonoInts/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs Ezfio_files MOs Nuclei Output Utils

13
src/MonoInts/README.rst Normal file
View File

@ -0,0 +1,13 @@
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_

121
src/MonoInts/ao_mono_ints.irp.f Normal file
View File

@ -0,0 +1,121 @@
BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Overlap between atomic basis functions:
! :math:`\int \chi_i(r) \chi_j(r) dr)`
END_DOC
integer :: i,j,n,l
double precision :: f
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
dim1=100
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,c) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_transp,ao_nucl, &
!$OMP ao_expo_transp,dim1)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
ao_overlap(i,j)= 0.d0
ao_overlap_x(i,j)= 0.d0
ao_overlap_y(i,j)= 0.d0
ao_overlap_z(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
c = ao_coef_transp(n,j) * ao_coef_transp(l,i)
ao_overlap(i,j) += c * overlap
ao_overlap_x(i,j) += c * overlap_x
ao_overlap_y(i,j) += c * overlap_y
ao_overlap_z(i,j) += c * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Overlap between absolute value of atomic basis functions:
! :math:`\int |\chi_i(r)| |\chi_j(r)| dr)`
END_DOC
integer :: i,j,n,l
double precision :: f
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx
dim1=100
lower_exp_val = 40.d0
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,dx) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_abs,ao_num,ao_coef_transp,ao_nucl, &
!$OMP ao_expo_transp,dim1,lower_exp_val)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
ao_overlap_abs(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_transp(l,i)
call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)
call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1)
call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1)
ao_overlap_abs(i,j) += abs(ao_coef_transp(n,j) * ao_coef_transp(l,i)) * overlap_x * overlap_y * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

31
src/MonoInts/check_orthonormality.irp.f Normal file
View File

@ -0,0 +1,31 @@
program check_ortho
implicit none
call do_print
print *, '--'
print *, '--'
call orthonormalize_mos
call do_print
end
subroutine do_print
implicit none
integer :: i,j
real :: off_diag, diag
off_diag = 0.
diag = 0.
do j=1,mo_tot_num
do i=1,mo_tot_num
off_diag += abs(mo_overlap(i,j))
enddo
diag += abs(mo_overlap(j,j))
off_diag -= abs(mo_overlap(j,j))
enddo
print *, 'Diag = ', abs(1.d0-diag/dble(mo_tot_num))
print *, 'Off-Diag = ', off_diag/(dble(mo_tot_num)**2-dble(mo_tot_num))
print *, '--'
end

16
src/MonoInts/dimensions.irp.f Normal file
View File

@ -0,0 +1,16 @@
BEGIN_PROVIDER [ integer, n_pt_max_integrals ]
&BEGIN_PROVIDER [ integer, n_pt_max_i_x]
implicit none
integer :: n_pt_sup
integer :: prim_power_l_max
include 'constants.F'
prim_power_l_max = maxval(ao_power)
n_pt_max_integrals = 24 * prim_power_l_max + 4
n_pt_max_i_x = 8 * prim_power_l_max
ASSERT (n_pt_max_i_x-1 <= max_dim)
if (n_pt_max_i_x-1 > max_dim) then
print *, 'Increase max_dim in constants.F to ', n_pt_max_i_x-1
stop 1
endif
END_PROVIDER

148
src/MonoInts/kin_ao_ints.irp.f Normal file
View File

@ -0,0 +1,148 @@
BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num_align,ao_num) ]
implicit none
integer :: i,j,n,l
double precision :: f
integer :: dim1
double precision :: overlap, overlap_y, overlap_z
double precision :: overlap_x0, overlap_y0, overlap_z0
double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: d_a_2,d_2
dim1=100
BEGIN_DOC
! second derivatives matrix elements in the ao basis
! .. math::
!
! {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
END_DOC
! -- Dummy call to provide everything
A_center(:) = 0.d0
B_center(:) = 1.d0
alpha = 1.d0
beta = .1d0
power_A = 1
power_B = 0
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1)
! --
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,c,d_a_2,d_2,deriv_tmp, &
!$OMP overlap_x0,overlap_y0,overlap_z0) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_transp,ao_nucl, &
!$OMP ao_expo_transp,dim1)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
ao_deriv2_x(i,j)= 0.d0
ao_deriv2_y(i,j)= 0.d0
ao_deriv2_z(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
c = ao_coef_transp(n,j) * ao_coef_transp(l,i)
! if (abs(c) < 1.d-8) then
! cycle
! endif
power_A(1) = power_A(1)-2
if (power_A(1)>-1) then
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_a_2,overlap_y,overlap_z,overlap,dim1)
else
d_a_2 = 0.d0
endif
power_A(1) = power_A(1)+4
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_2,overlap_y,overlap_z,overlap,dim1)
power_A(1) = power_A(1)-2
double precision :: deriv_tmp
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(1) +1.d0) * overlap_x0 &
+power_A(1) * (power_A(1)-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0
ao_deriv2_x(i,j) += c*deriv_tmp
power_A(2) = power_A(2)-2
if (power_A(2)>-1) then
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1)
else
d_a_2 = 0.d0
endif
power_A(2) = power_A(2)+4
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1)
power_A(2) = power_A(2)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(2) +1.d0 ) * overlap_y0 &
+power_A(2) * (power_A(2)-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0
ao_deriv2_y(i,j) += c*deriv_tmp
power_A(3) = power_A(3)-2
if (power_A(3)>-1) then
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_a_2,overlap,dim1)
else
d_a_2 = 0.d0
endif
power_A(3) = power_A(3)+4
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1)
power_A(3) = power_A(3)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(3) +1.d0 ) * overlap_z0 &
+power_A(3) * (power_A(3)-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0
ao_deriv2_z(i,j) += c*deriv_tmp
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
implicit none
BEGIN_DOC
! array of the priminitve basis kinetic integrals
! \langle \chi_i |\hat{T}| \chi_j \rangle
END_DOC
integer :: i,j,k,l
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i,j) &
!$OMP SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
do j = 1, ao_num
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do i = 1, ao_num
ao_kinetic_integral(i,j) = -0.5d0 * (ao_deriv2_x(i,j) + ao_deriv2_y(i,j) + ao_deriv2_z(i,j) )
enddo
do i = ao_num +1,ao_num_align
ao_kinetic_integral(i,j) = 0.d0
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

25
src/MonoInts/kin_mo_ints.irp.f Normal file
View File

@ -0,0 +1,25 @@
BEGIN_PROVIDER [double precision, mo_kinetic_integral, (mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i1,j1,i,j
double precision :: c_i1
mo_kinetic_integral = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j,i1,j1,c_i1) &
!$OMP SHARED(mo_tot_num,ao_num,mo_coef,ao_Kinetic_integral, &
!$OMP mo_kinetic_integral)
do i = 1,mo_tot_num
do j = 1,mo_tot_num
do i1 = 1,ao_num
c_i1 = mo_coef(i1,i)
!DIR$ VECTOR ALIGNED
do j1 = 1,ao_num
mo_kinetic_integral(j,i) = mo_kinetic_integral(j,i) + c_i1*mo_coef(j1,j) *&
ao_Kinetic_integral(j1,i1)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

47
src/MonoInts/mo_mono_ints.irp.f Normal file
View File

@ -0,0 +1,47 @@
BEGIN_PROVIDER [ double precision, mo_overlap,(mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i,j,n,l
double precision :: f
integer :: lmax
lmax = iand(ao_num,-4)
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) &
!$OMP PRIVATE(i,j,n,l) &
!$OMP SHARED(mo_overlap,mo_coef,ao_overlap, &
!$OMP mo_tot_num,ao_num,lmax)
do j=1,mo_tot_num
do i= 1,mo_tot_num
mo_overlap(i,j) = 0.d0
do n = 1, lmax,4
!DIR$ VECTOR ALIGNED
do l = 1, ao_num
mo_overlap(i,j) = mo_overlap(i,j) + mo_coef(l,i) * &
( mo_coef(n ,j) * ao_overlap(l,n ) &
+ mo_coef(n+1,j) * ao_overlap(l,n+1) &
+ mo_coef(n+2,j) * ao_overlap(l,n+2) &
+ mo_coef(n+3,j) * ao_overlap(l,n+3) )
enddo
enddo
do n = lmax+1, ao_num
!DIR$ VECTOR ALIGNED
do l = 1, ao_num
mo_overlap(i,j) = mo_overlap(i,j) + mo_coef(n,j) * mo_coef(l,i) * ao_overlap(l,n)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i,j,n,l
BEGIN_DOC
! array of the mono electronic hamiltonian on the MOs basis
! : sum of the kinetic and nuclear electronic potential
END_DOC
do i = 1, mo_tot_num
do j = 1, mo_tot_num
mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j)
enddo
enddo
END_PROVIDER

11
src/MonoInts/orthonormalize.irp.f Normal file
View File

@ -0,0 +1,11 @@
subroutine orthonormalize_mos
implicit none
integer :: m,p,s
m = size(mo_coef,1)
p = size(mo_overlap,1)
call ortho_lowdin(mo_overlap,p,mo_tot_num,mo_coef,m,ao_num)
mo_label = 'Orthonormalized'
SOFT_TOUCH mo_coef mo_label
end

520
src/MonoInts/pot_ao_ints.irp.f Normal file
View File

@ -0,0 +1,520 @@
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)]
BEGIN_DOC
! interaction nuclear electron
END_DOC
implicit none
double precision :: alpha, beta, gama, delta
integer :: i_c,num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt_in,m
double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
integer :: nucl_numC
! Important for OpenMP
PROVIDE all_utils
ao_nucl_elec_integral = 0.d0
! -- Dummy run to provide everything called by NAI_pol_mult
!power_A(:)= ao_power(1,:)
!A_center(:) = nucl_coord(1,:)
!B_center(:) = nucl_coord(1,:)+0.5d0
!C_center(:) = 0.5d0*(B_center(:)-A_center(:))
!alpha = 1.d0
!beta = 0.1d0
!Z = 1.d0
!n_pt_in = n_pt_max_integrals
!c = NAI_pol_mult(A_center,B_center,power_A,power_A,alpha,beta,C_center,n_pt_in)
!c = NAI_pol_mult(A_center,A_center,power_A,power_A,alpha,beta,C_center,n_pt_in)
!c = NAI_pol_mult(A_center,A_center,power_A,power_A,alpha,beta,A_center,n_pt_in)
! --
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
!$OMP num_A,num_B,Z,c,n_pt_in) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
!$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge)
n_pt_in = n_pt_max_integrals
!$OMP DO SCHEDULE (guided)
do j = 1, ao_num
power_A(1)= ao_power(j,1)
power_A(2)= ao_power(j,2)
power_A(3)= ao_power(j,3)
num_A = ao_nucl(j)
A_center(1) = nucl_coord(num_A,1)
A_center(2) = nucl_coord(num_A,2)
A_center(3) = nucl_coord(num_A,3)
do i = 1, ao_num
power_B(1)= ao_power(i,1)
power_B(2)= ao_power(i,2)
power_B(3)= ao_power(i,3)
num_B = ao_nucl(i)
B_center(1) = nucl_coord(num_B,1)
B_center(2) = nucl_coord(num_B,2)
B_center(3) = nucl_coord(num_B,3)
do l=1,ao_prim_num(j)
alpha = ao_expo_transp(l,j)
do m=1,ao_prim_num(i)
beta = ao_expo_transp(m,i)
c = 0.d0
do k = 1, nucl_num
double precision :: Z,c
Z = nucl_charge(k)
C_center(1) = nucl_coord(k,1)
C_center(2) = nucl_coord(k,2)
C_center(3) = nucl_coord(k,3)
c = c+Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)
enddo
ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - &
ao_coef_transp(l,j)*ao_coef_transp(m,i)*c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)
! function that calculate the folowing integral :
! int{dr} of (x-A_x)^ax (x-B_X)^bx exp(-alpha (x-A_x)^2 - beta (x-B_x)^2 ) 1/(r-R_c)
implicit none
double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt
double precision :: P_center(3)
double precision :: d(0:n_pt_in),pouet,coeff,rho,dist,const,pouet_2,p,p_inv,factor
double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi
double precision :: V_e_n,const_factor,dist_integral,tmp
include 'constants.F'
if ( (A_center(1)/=B_center(1)).or. &
(A_center(2)/=B_center(2)).or. &
(A_center(3)/=B_center(3)).or. &
(A_center(1)/=C_center(1)).or. &
(A_center(2)/=C_center(2)).or. &
(A_center(3)/=C_center(3))) then
continue
else
NAI_pol_mult = V_e_n(power_A(1),power_A(2),power_A(3),power_B(1),power_B(2),power_B(3),alpha,beta)
return
endif
p = alpha + beta
! print*, "a"
p_inv = 1.d0/p
rho = alpha * beta * p_inv
dist = 0.d0
dist_integral = 0.d0
do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
dist += (A_center(i) - B_center(i))*(A_center(i) - B_center(i))
dist_integral += (P_center(i) - C_center(i))*(P_center(i) - C_center(i))
enddo
const_factor = dist*rho
const = p * dist_integral
factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv
lmax = 20
! print*, "b"
do i = 0, n_pt_in
d(i) = 0.d0
enddo
n_pt = 2 * ( (power_A(1) + power_B(1)) +(power_A(2) + power_B(2)) +(power_A(3) + power_B(3)) )
if (n_pt == 0) then
epsilo = 1.d0
pouet = rint(0,const)
NAI_pol_mult = coeff * pouet
return
endif
call give_polynom_mult_center_mono_elec(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out)
if(n_pt_out<0)then
NAI_pol_mult = 0.d0
return
endif
double precision :: accu,epsilo,rint
integer :: n_pt_in,n_pt_out,lmax
accu = 0.d0
! 1/r1 standard attraction integral
epsilo = 1.d0
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
do i =0 ,n_pt_out,2
accu += d(i) * rint(i/2,const)
enddo
NAI_pol_mult = accu * coeff
end
subroutine give_polynom_mult_center_mono_elec(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out)
!!!! subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw ::
!!!! I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q)
!!!! it is for the nuclear electron atraction
implicit none
integer, intent(in) :: n_pt_in
integer,intent(out) :: n_pt_out
double precision, intent(in) :: A_center(3), B_center(3),C_center(3)
double precision, intent(in) :: alpha,beta
integer, intent(in) :: power_A(3), power_B(3)
integer :: a_x,b_x,a_y,b_y,a_z,b_z
double precision :: d(0:n_pt_in)
double precision :: d1(0:n_pt_in)
double precision :: d2(0:n_pt_in)
double precision :: d3(0:n_pt_in)
double precision :: accu, pq_inv, p10_1, p10_2, p01_1, p01_2
double precision :: p,P_center(3),rho,p_inv,p_inv_2
!print*,'n_pt_in = ',n_pt_in
accu = 0.d0
!COMPTEUR irp_rdtsc1 = irp_rdtsc()
ASSERT (n_pt_in > 1)
p = alpha+beta
p_inv = 1.d0/p
p_inv_2 = 0.5d0/p
do i =1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
enddo
! print*,'passed the P_center'
double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
R1x(0) = (P_center(1) - A_center(1))
R1x(1) = 0.d0
R1x(2) = -(P_center(1) - C_center(1))
! R1x = (P_x - A_x) - (P_x - C_x) t^2
R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(1) - C_center(1))
!R1xp = (P_x - B_x) - (P_x - C_x) t^2
R2x(0) = p_inv_2
R2x(1) = 0.d0
R2x(2) = -p_inv_2
!R2x = 0.5 / p - 0.5/p t^2
do i = 0,n_pt_in
d(i) = 0.d0
enddo
do i = 0,n_pt_in
d1(i) = 0.d0
enddo
do i = 0,n_pt_in
d2(i) = 0.d0
enddo
do i = 0,n_pt_in
d3(i) = 0.d0
enddo
integer :: n_pt1,n_pt2,n_pt3,dim,i
n_pt1 = n_pt_in
n_pt2 = n_pt_in
n_pt3 = n_pt_in
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_mono_elec(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in)
! print*,'passed the first I_x1'
if(n_pt1<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(2) - A_center(2))
R1x(1) = 0.d0
R1x(2) = -(P_center(2) - C_center(2))
! R1x = (P_x - A_x) - (P_x - C_x) t^2
R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(2) - C_center(2))
!R1xp = (P_x - B_x) - (P_x - C_x) t^2
a_y = power_A(2)
b_y = power_B(2)
call I_x1_pol_mult_mono_elec(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in)
! print*,'passed the second I_x1'
if(n_pt2<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(3) - A_center(3))
R1x(1) = 0.d0
R1x(2) = -(P_center(3) - C_center(3))
! R1x = (P_x - A_x) - (P_x - C_x) t^2
R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(3) - C_center(3))
!R2x = 0.5 / p - 0.5/p t^2
a_z = power_A(3)
b_z = power_B(3)
! print*,'a_z = ',a_z
! print*,'b_z = ',b_z
call I_x1_pol_mult_mono_elec(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in)
! print*,'passed the third I_x1'
if(n_pt3<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
integer :: n_pt_tmp
n_pt_tmp = 0
call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp)
do i = 0,n_pt_tmp
d1(i) = 0.d0
enddo
n_pt_out = 0
call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out)
do i = 0, n_pt_out
d(i) = d1(i)
enddo
end
recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
!!!! recursive function involved in the electron nucleus potential
implicit none
integer , intent(in) :: n_pt_in
double precision,intent(inout) :: d(0:n_pt_in)
integer,intent(inout) :: nd
integer, intent(in):: a,c
double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2)
include 'constants.F'
double precision :: X(0:max_dim)
double precision :: Y(0:max_dim)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
integer :: nx, ix,dim,iy,ny
dim = n_pt_in
! print*,'a,c = ',a,c
! print*,'nd_in = ',nd
if( (a==0) .and. (c==0))then
! print*,'coucou !'
nd = 0
d(0) = 1.d0
return
elseif( (c<0).or.(nd<0) )then
nd = -1
return
else if ((a==0).and.(c.ne.0)) then
call I_x2_pol_mult_mono_elec(c,R1x,R1xp,R2x,d,nd,n_pt_in)
! print*,'nd 0,c',nd
else if (a==1) then
nx = nd
do ix=0,n_pt_in
X(ix) = 0.d0
Y(ix) = 0.d0
enddo
call I_x2_pol_mult_mono_elec(c-1,R1x,R1xp,R2x,X,nx,n_pt_in)
do ix=0,nx
X(ix) *= c
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
ny=0
call I_x2_pol_mult_mono_elec(c,R1x,R1xp,R2x,Y,ny,n_pt_in)
call multiply_poly(Y,ny,R1x,2,d,nd)
else
do ix=0,n_pt_in
X(ix) = 0.d0
Y(ix) = 0.d0
enddo
nx = 0
call I_x1_pol_mult_mono_elec(a-2,c,R1x,R1xp,R2x,X,nx,n_pt_in)
! print*,'nx a-2,c= ',nx
do ix=0,nx
X(ix) *= a-1
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
! print*,'nd out = ',nd
nx = nd
do ix=0,n_pt_in
X(ix) = 0.d0
enddo
call I_x1_pol_mult_mono_elec(a-1,c-1,R1x,R1xp,R2x,X,nx,n_pt_in)
! print*,'nx a-1,c-1 = ',nx
do ix=0,nx
X(ix) *= c
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
ny=0
call I_x1_pol_mult_mono_elec(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in)
call multiply_poly(Y,ny,R1x,2,d,nd)
endif
end
recursive subroutine I_x2_pol_mult_mono_elec(c,R1x,R1xp,R2x,d,nd,dim)
implicit none
integer , intent(in) :: dim
include 'constants.F'
double precision :: d(0:max_dim)
integer,intent(inout) :: nd
integer, intent(in):: c
double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2)
integer :: i
!print*,'X2,c = ',c
!print*,'nd_in = ',nd
if(c==0) then
nd = 0
d(0) = 1.d0
! print*,'nd IX2 = ',nd
return
elseif ((nd<0).or.(c<0))then
nd = -1
return
else
integer :: nx, ix,ny
double precision :: X(0:max_dim),Y(0:max_dim)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
do ix=0,dim
X(ix) = 0.d0
Y(ix) = 0.d0
enddo
nx = 0
call I_x1_pol_mult_mono_elec(0,c-2,R1x,R1xp,R2x,X,nx,dim)
! print*,'nx 0,c-2 = ',nx
do ix=0,nx
X(ix) *= c-1
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
! print*,'nd = ',nd
ny = 0
do ix=0,dim
Y(ix) = 0.d0
enddo
call I_x1_pol_mult_mono_elec(0,c-1,R1x,R1xp,R2x,Y,ny,dim)
! print*,'ny = ',ny
! do ix=0,ny
! print*,'Y(ix) = ',Y(ix)
! enddo
if(ny.ge.0)then
call multiply_poly(Y,ny,R1xp,2,d,nd)
endif
endif
end
double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
implicit none
!!! primitve nuclear attraction between the two primitves centered on the same atom ::
!!!! primitive_1 = x**(a_x) y**(a_y) z**(a_z) exp(-alpha * r**2)
!!!! primitive_2 = x**(b_x) y**(b_y) z**(b_z) exp(- beta * r**2)
integer :: a_x,a_y,a_z,b_x,b_y,b_z
double precision :: alpha,beta
double precision :: V_r, V_phi, V_theta
if(iand((a_x+b_x),1)==1.or.iand(a_y+b_y,1)==1.or.iand((a_z+b_z),1)==1)then
V_e_n = 0.d0
else
V_e_n = V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta) &
& * V_phi(a_x+b_x,a_y+b_y) &
& * V_theta(a_z+b_z,a_x+b_x+a_y+b_y+1)
endif
end
double precision function int_gaus_pol(alpha,n)
!!!! calculate the integral of
!! integral on "x" with boundaries (- infinity; + infinity) of [ x**n exp(-alpha * x**2) ]
implicit none
double precision :: alpha
integer :: n
double precision :: dble_fact
include 'constants.F'
!if(iand(n,1).eq.1)then
! int_gaus_pol= 0.d0
!else
! int_gaus_pol = dsqrt(pi/alpha) * dble_fact(n -1)/(alpha+alpha)**(n/2)
!endif
int_gaus_pol = 0.d0
if(iand(n,1).eq.0)then
int_gaus_pol = dsqrt(alpha/pi)
double precision :: two_alpha
two_alpha = alpha+alpha
integer :: i
do i=1,n,2
int_gaus_pol = int_gaus_pol * two_alpha
enddo
int_gaus_pol = dble_fact(n -1) / int_gaus_pol
endif
end
double precision function V_r(n,alpha)
!!!! calculate the radial part of the nuclear attraction integral which is the following integral :
!! integral on "r" with boundaries ( 0 ; + infinity) of [ r**n exp(-alpha * r**2) ]
!!! CAUTION :: this function requires the constant sqpi = dsqrt(pi)
implicit none
double precision :: alpha, fact
integer :: n
include 'constants.F'
if(iand(n,1).eq.1)then
V_r = 0.5d0 * fact(ishft(n,-1)) / (alpha ** (ishft(n,-1) + 1))
else
V_r = sqpi * fact(n) / fact(ishft(n,-1)) * (0.5d0/sqrt(alpha)) ** (n+1)
endif
end
double precision function V_phi(n,m)
implicit none
!!!! calculate the angular "phi" part of the nuclear attraction integral wich is the following integral :
!! integral on "phi" with boundaries ( 0 ; 2 pi) of [ cos(phi) **n sin(phi) **m ]
integer :: n,m, i
double precision :: prod, Wallis
prod = 1.d0
do i = 0,ishft(n,-1)-1
prod = prod/ (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
enddo
V_phi = 4.d0 * prod * Wallis(m)
end
double precision function V_theta(n,m)
implicit none
!!!! calculate the angular "theta" part of the nuclear attraction integral wich is the following integral :
!! integral on "theta" with boundaries ( 0 ; pi) of [ cos(theta) **n sin(theta) **m ]
integer :: n,m,i
double precision :: Wallis, prod
include 'constants.F'
V_theta = 0.d0
prod = 1.d0
do i = 0,ishft(n,-1)-1
prod = prod / (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
enddo
V_theta = (prod+prod) * Wallis(m)
end
double precision function Wallis(n)
!!!! calculate the Wallis integral :
!! integral on "theta" with boundaries ( 0 ; pi/2) of [ cos(theta) **n ]
implicit none
double precision :: fact
integer :: n,p
include 'constants.F'
if(iand(n,1).eq.0)then
Wallis = fact(ishft(n,-1))
Wallis = pi * fact(n) / (dble(ibset(0_8,n)) * (Wallis+Wallis)*Wallis)
else
p = ishft(n,-1)
Wallis = fact(p)
Wallis = dble(ibset(0_8,p+p)) * Wallis*Wallis / fact(p+p+1)
endif
end

25
src/MonoInts/pot_mo_ints.irp.f Normal file
View File

@ -0,0 +1,25 @@
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i1,j1,i,j
double precision :: c_i1,c_j1
mo_nucl_elec_integral = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
!$OMP SHARED(mo_tot_num,ao_num,mo_coef, &
!$OMP mo_nucl_elec_integral, ao_nucl_elec_integral)
do i = 1, mo_tot_num
do j = 1, mo_tot_num
do i1 = 1,ao_num
c_i1 = mo_coef(i1,i)
do j1 = 1,ao_num
c_j1 = c_i1*mo_coef(j1,j)
mo_nucl_elec_integral(j,i) = mo_nucl_elec_integral(j,i) + &
c_j1 * ao_nucl_elec_integral(j1,i1)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

4
src/MonoInts/save_ortho_mos.irp.f Normal file
View File

@ -0,0 +1,4 @@
program save_ortho_mos
call orthonormalize_mos
call save_mos
end

29
src/MonoInts/tests/Makefile Normal file
View File

@ -0,0 +1,29 @@
IRPF90+= -I tests
REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f))
.PHONY: clean executables serial_tests parallel_tests
all: clean executables serial_tests parallel_tests
parallel_tests: $(REF_FILES)
@echo ; echo " ---- Running parallel tests ----" ; echo
@OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py
serial_tests: $(REF_FILES)
@echo ; echo " ---- Running serial tests ----" ; echo
@OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py
executables: $(wildcard *.irp.f) veryclean
$(MAKE) -C ..
%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables
$(QPACKAGE_ROOT)/scripts/create_test_ref.sh $*
clean:
$(MAKE) -C .. clean
veryclean:
$(MAKE) -C .. veryclean

12
src/MonoInts/tests/ao_kin.irp.f Normal file
View File

@ -0,0 +1,12 @@
program test ao_kin
implicit none
integer :: i, j
do i=1,ao_num
do j=1,i
if (abs(ao_kinetic_integral(i,j)) > 1.d-9) then
print '(I4,A,I4,2x,A, 2x,e16.8)' , i,',',j, ' : ', ao_kinetic_integral(i,j)
endif
enddo
enddo
end

50554
src/MonoInts/tests/ao_kin.ref Normal file
View File

File diff suppressed because it is too large Load Diff

12
src/MonoInts/tests/ao_mono.irp.f Normal file
View File

@ -0,0 +1,12 @@
program test ao_mono
implicit none
integer :: i, j
do i=1,ao_num
do j=1,i
if (abs(ao_overlap(i,j)) > 1.d-9) then
print '(I4,A,I4,2x,A, 2x,e16.8)' , i,',',j, ' : ', ao_overlap(i,j)
endif
enddo
enddo
end

49142
src/MonoInts/tests/ao_mono.ref Normal file
View File

File diff suppressed because it is too large Load Diff

12
src/MonoInts/tests/ao_pot.irp.f Normal file
View File

@ -0,0 +1,12 @@
program test ao_pot
implicit none
integer :: i, j
do i=1,ao_num
do j=1,i
if (abs(ao_nucl_elec_integral(i,j)) > 1.d-9) then
print '(I4,A,I4,2x,A, 2x,e16.8)' , i,',',j, ' : ', ao_nucl_elec_integral(i,j)
endif
enddo
enddo
end

64451
src/MonoInts/tests/ao_pot.ref Normal file
View File

File diff suppressed because it is too large Load Diff

12
src/MonoInts/tests/mo_kin.irp.f Normal file
View File

@ -0,0 +1,12 @@
program test mo_kin
implicit none
integer :: i, j
do i=1,mo_tot_num
do j=1,i
if (abs(mo_kinetic_integral(i,j)) > 1.d-9) then
print '(I4,A,I4,2x,A, 2x,e16.8)' , i,',',j, ' : ', mo_kinetic_integral(i,j)
endif
enddo
enddo
end

44717
src/MonoInts/tests/mo_kin.ref Normal file
View File

File diff suppressed because it is too large Load Diff

20
src/MonoInts/tests/mo_mono.irp.f Normal file
View File

@ -0,0 +1,20 @@
program test mo_mono
implicit none
integer :: i, j
integer :: ortho
ortho = 1
do i=1,mo_tot_num
do j=1,i-1
if (abs(mo_overlap(i,j)) > 1.d-9) then
print '(I4,A,I4,2x,A, 2x,e16.8)' , i,',',j, ' : ', ao_overlap(i,j)
ortho = 0
endif
enddo
if (abs(1.d0-mo_overlap(i,j)) > 1.d-9) then
print '(I4,A,I4,2x,A, 2x,e16.8)' , i,',',j, ' : ', ao_overlap(i,j)
ortho = 0
endif
enddo
print *, 'Orthonormal : ', ortho
end

635
src/MonoInts/tests/mo_mono.ref Normal file
View File

@ -0,0 +1,635 @@
data = {
'AlCl.ezfio' : {
'Orthonormal' : 1,
},
'Al.ezfio' : {
'Orthonormal' : 1,
},
'Al+.ezfio' : {
'Orthonormal' : 1,
},
'AlH2.ezfio' : {
'Orthonormal' : 1,
},
'AlH3.ezfio' : {
'Orthonormal' : 1,
},
'AlH.ezfio' : {
'Orthonormal' : 1,
},
'BCl.ezfio' : {
'Orthonormal' : 1,
},
'BeCl.ezfio' : {
'Orthonormal' : 1,
},
'Be.ezfio' : {
'Orthonormal' : 1,
},
'Be+.ezfio' : {
'Orthonormal' : 1,
},
'BeF.ezfio' : {
'Orthonormal' : 1,
},
'BeH2.ezfio' : {
'Orthonormal' : 1,
},
'BeH.ezfio' : {
'Orthonormal' : 1,
},
'BeO.ezfio' : {
'Orthonormal' : 1,
},
'BeOH.ezfio' : {
'Orthonormal' : 1,
},
'BeS.ezfio' : {
'Orthonormal' : 1,
},
'B.ezfio' : {
'Orthonormal' : 1,
},
'B+.ezfio' : {
'Orthonormal' : 1,
},
'BH2.ezfio' : {
'Orthonormal' : 1,
},
'BH3.ezfio' : {
'Orthonormal' : 1,
},
'BH.ezfio' : {
'Orthonormal' : 1,
},
'BO.ezfio' : {
'Orthonormal' : 1,
},
'BS.ezfio' : {
'Orthonormal' : 1,
},
'C2.ezfio' : {
'Orthonormal' : 1,
},
'C2H2.ezfio' : {
'Orthonormal' : 1,
},
'C2H2+.ezfio' : {
'Orthonormal' : 1,
},
'C2H3.ezfio' : {
'Orthonormal' : 1,
},
'C2H3+.ezfio' : {
'Orthonormal' : 1,
},
'C2H4.ezfio' : {
'Orthonormal' : 1,
},
'C2H4+.ezfio' : {
'Orthonormal' : 1,
},
'C2H5.ezfio' : {
'Orthonormal' : 1,
},
'C2H6.ezfio' : {
'Orthonormal' : 1,
},
'C2H.ezfio' : {
'Orthonormal' : 1,
},
'CCl.ezfio' : {
'Orthonormal' : 1,
},
'C-.ezfio' : {
'Orthonormal' : 1,
},
'C.ezfio' : {
'Orthonormal' : 1,
},
'C+.ezfio' : {
'Orthonormal' : 1,
},
'CF.ezfio' : {
'Orthonormal' : 1,
},
'CH2_1A1.ezfio' : {
'Orthonormal' : 1,
},
'CH2_3B1.ezfio' : {
'Orthonormal' : 1,
},
'CH2-.ezfio' : {
'Orthonormal' : 1,
},
'CH3-.ezfio' : {
'Orthonormal' : 1,
},
'CH3.ezfio' : {
'Orthonormal' : 1,
},
'CH4.ezfio' : {
'Orthonormal' : 1,
},
'CH4+.ezfio' : {
'Orthonormal' : 1,
},
'CH-.ezfio' : {
'Orthonormal' : 1,
},
'CH.ezfio' : {
'Orthonormal' : 1,
},
'Cl2-.ezfio' : {
'Orthonormal' : 1,
},
'Cl2.ezfio' : {
'Orthonormal' : 1,
},
'Cl2+.ezfio' : {
'Orthonormal' : 1,
},
'Cl-.ezfio' : {
'Orthonormal' : 1,
},
'Cl.ezfio' : {
'Orthonormal' : 1,
},
'Cl+.ezfio' : {
'Orthonormal' : 1,
},
'ClH2+.ezfio' : {
'Orthonormal' : 1,
},
'ClH.ezfio' : {
'Orthonormal' : 1,
},
'ClH+.ezfio' : {
'Orthonormal' : 1,
},
'ClS.ezfio' : {
'Orthonormal' : 1,
},
'ClSiH3.ezfio' : {
'Orthonormal' : 1,
},
'CN-.ezfio' : {
'Orthonormal' : 1,
},
'CN.ezfio' : {
'Orthonormal' : 1,
},
'CO2.ezfio' : {
'Orthonormal' : 1,
},
'CO.ezfio' : {
'Orthonormal' : 1,
},
'CO+.ezfio' : {
'Orthonormal' : 1,
},
'COS.ezfio' : {
'Orthonormal' : 1,
},
'CP.ezfio' : {
'Orthonormal' : 1,
},
'CS2.ezfio' : {
'Orthonormal' : 1,
},
'CS.ezfio' : {
'Orthonormal' : 1,
},
'CS+.ezfio' : {
'Orthonormal' : 1,
},
'CSi.ezfio' : {
'Orthonormal' : 1,
},
'F2.ezfio' : {
'Orthonormal' : 1,
},
'FAl.ezfio' : {
'Orthonormal' : 1,
},
'FCl.ezfio' : {
'Orthonormal' : 1,
},
'FCl+.ezfio' : {
'Orthonormal' : 1,
},
'F-.ezfio' : {
'Orthonormal' : 1,
},
'F.ezfio' : {
'Orthonormal' : 1,
},
'F+.ezfio' : {
'Orthonormal' : 1,
},
'FH.ezfio' : {
'Orthonormal' : 1,
},
'FH+.ezfio' : {
'Orthonormal' : 1,
},
'FMg.ezfio' : {
'Orthonormal' : 1,
},
'FNa.ezfio' : {
'Orthonormal' : 1,
},
'FP.ezfio' : {
'Orthonormal' : 1,
},
'FS.ezfio' : {
'Orthonormal' : 1,
},
'FSi.ezfio' : {
'Orthonormal' : 1,
},
'FSiH3.ezfio' : {
'Orthonormal' : 1,
},
'H2CNH.ezfio' : {
'Orthonormal' : 1,
},
'H2CO.ezfio' : {
'Orthonormal' : 1,
},
'H2CPH.ezfio' : {
'Orthonormal' : 1,
},
'H2CS.ezfio' : {
'Orthonormal' : 1,
},
'H2.ezfio' : {
'Orthonormal' : 1,
},
'H2NNH2.ezfio' : {
'Orthonormal' : 1,
},
'H2PPH2.ezfio' : {
'Orthonormal' : 1,
},
'H3CCl.ezfio' : {
'Orthonormal' : 1,
},
'H3CF.ezfio' : {
'Orthonormal' : 1,
},
'H3CNH2.ezfio' : {
'Orthonormal' : 1,
},
'H3COH.ezfio' : {
'Orthonormal' : 1,
},
'H3CSH.ezfio' : {
'Orthonormal' : 1,
},
'H3SiSiH3.ezfio' : {
'Orthonormal' : 1,
},
'HBO.ezfio' : {
'Orthonormal' : 1,
},
'HBS.ezfio' : {
'Orthonormal' : 1,
},
'HCF.ezfio' : {
'Orthonormal' : 1,
},
'HCN.ezfio' : {
'Orthonormal' : 1,
},
'HCO.ezfio' : {
'Orthonormal' : 1,
},
'HCP.ezfio' : {
'Orthonormal' : 1,
},
'H.ezfio' : {
'Orthonormal' : 1,
},
'HNO.ezfio' : {
'Orthonormal' : 1,
},
'HOCl.ezfio' : {
'Orthonormal' : 1,
},
'HOF.ezfio' : {
'Orthonormal' : 1,
},
'HOMg.ezfio' : {
'Orthonormal' : 1,
},
'HONa.ezfio' : {
'Orthonormal' : 1,
},
'HOO.ezfio' : {
'Orthonormal' : 1,
},
'HOOH.ezfio' : {
'Orthonormal' : 1,
},
'HSSH.ezfio' : {
'Orthonormal' : 1,
},
'Li2.ezfio' : {
'Orthonormal' : 1,
},
'LiCl.ezfio' : {
'Orthonormal' : 1,
},
'Li.ezfio' : {
'Orthonormal' : 1,
},
'Li+.ezfio' : {
'Orthonormal' : 1,
},
'LiF.ezfio' : {
'Orthonormal' : 1,
},
'LiH.ezfio' : {
'Orthonormal' : 1,
},
'LiN.ezfio' : {
'Orthonormal' : 1,
},
'LiO.ezfio' : {
'Orthonormal' : 1,
},
'LiOH.ezfio' : {
'Orthonormal' : 1,
},
'MgCl.ezfio' : {
'Orthonormal' : 1,
},
'Mg.ezfio' : {
'Orthonormal' : 1,
},
'Mg+.ezfio' : {
'Orthonormal' : 1,
},
'MgH.ezfio' : {
'Orthonormal' : 1,
},
'MgS.ezfio' : {
'Orthonormal' : 1,
},
'N2.ezfio' : {
'Orthonormal' : 1,
},
'N2+.ezfio' : {
'Orthonormal' : 1,
},
'Na2.ezfio' : {
'Orthonormal' : 1,
},
'NaCl.ezfio' : {
'Orthonormal' : 1,
},
'Na.ezfio' : {
'Orthonormal' : 1,
},
'Na+.ezfio' : {
'Orthonormal' : 1,
},
'NaH.ezfio' : {
'Orthonormal' : 1,
},
'N.ezfio' : {
'Orthonormal' : 1,
},
'N+.ezfio' : {
'Orthonormal' : 1,
},
'NF.ezfio' : {
'Orthonormal' : 1,
},
'NH2-.ezfio' : {
'Orthonormal' : 1,
},
'NH2.ezfio' : {
'Orthonormal' : 1,
},
'NH3.ezfio' : {
'Orthonormal' : 1,
},
'NH3+.ezfio' : {
'Orthonormal' : 1,
},
'NH4+.ezfio' : {
'Orthonormal' : 1,
},
'NH-.ezfio' : {
'Orthonormal' : 1,
},
'NH.ezfio' : {
'Orthonormal' : 1,
},
'NO-.ezfio' : {
'Orthonormal' : 1,
},
'NO.ezfio' : {
'Orthonormal' : 1,
},
'NP.ezfio' : {
'Orthonormal' : 1,
},
'NS.ezfio' : {
'Orthonormal' : 1,
},
'NSi.ezfio' : {
'Orthonormal' : 1,
},
'O2Cl.ezfio' : {
'Orthonormal' : 1,
},
'O2-.ezfio' : {
'Orthonormal' : 1,
},
'O2.ezfio' : {
'Orthonormal' : 1,
},
'O2+.ezfio' : {
'Orthonormal' : 1,
},
'O2S.ezfio' : {
'Orthonormal' : 1,
},
'O2Si.ezfio' : {
'Orthonormal' : 1,
},
'O3.ezfio' : {
'Orthonormal' : 1,
},
'OCl.ezfio' : {
'Orthonormal' : 1,
},
'O-.ezfio' : {
'Orthonormal' : 1,
},
'O.ezfio' : {
'Orthonormal' : 1,
},
'O+.ezfio' : {
'Orthonormal' : 1,
},
'OH2.ezfio' : {
'Orthonormal' : 1,
},
'OH2+.ezfio' : {
'Orthonormal' : 1,
},
'OH3+.ezfio' : {
'Orthonormal' : 1,
},
'OH-.ezfio' : {
'Orthonormal' : 1,
},
'OH.ezfio' : {
'Orthonormal' : 1,
},
'OH+.ezfio' : {
'Orthonormal' : 1,
},
'OMg.ezfio' : {
'Orthonormal' : 1,
},
'ONa.ezfio' : {
'Orthonormal' : 1,
},
'OP-.ezfio' : {
'Orthonormal' : 1,
},
'OP.ezfio' : {
'Orthonormal' : 1,
},
'OPH.ezfio' : {
'Orthonormal' : 1,
},
'OS.ezfio' : {
'Orthonormal' : 1,
},
'OSi.ezfio' : {
'Orthonormal' : 1,
},
'P2.ezfio' : {
'Orthonormal' : 1,
},
'P2+.ezfio' : {
'Orthonormal' : 1,
},
'PCl.ezfio' : {
'Orthonormal' : 1,
},
'P-.ezfio' : {
'Orthonormal' : 1,
},
'P.ezfio' : {
'Orthonormal' : 1,
},
'PH2-.ezfio' : {
'Orthonormal' : 1,
},
'PH2.ezfio' : {
'Orthonormal' : 1,
},
'PH2+.ezfio' : {
'Orthonormal' : 1,
},
'PH3.ezfio' : {
'Orthonormal' : 1,
},
'PH3+.ezfio' : {
'Orthonormal' : 1,
},
'PH4+.ezfio' : {
'Orthonormal' : 1,
},
'PH-.ezfio' : {
'Orthonormal' : 1,
},
'PH.ezfio' : {
'Orthonormal' : 1,
},
'PS.ezfio' : {
'Orthonormal' : 1,
},
'S2-.ezfio' : {
'Orthonormal' : 1,
},
'S2.ezfio' : {
'Orthonormal' : 1,
},
'S-.ezfio' : {
'Orthonormal' : 1,
},
'S.ezfio' : {
'Orthonormal' : 1,
},
'S+.ezfio' : {
'Orthonormal' : 1,
},
'SH2.ezfio' : {
'Orthonormal' : 1,
},
'SH2+.ezfio' : {
'Orthonormal' : 1,
},
'SH3+.ezfio' : {
'Orthonormal' : 1,
},
'SH-.ezfio' : {
'Orthonormal' : 1,
},
'SH.ezfio' : {
'Orthonormal' : 1,
},
'SH+.ezfio' : {
'Orthonormal' : 1,
},
'Si2.ezfio' : {
'Orthonormal' : 1,
},
'SiCl.ezfio' : {
'Orthonormal' : 1,
},
'Si-.ezfio' : {
'Orthonormal' : 1,
},
'Si.ezfio' : {
'Orthonormal' : 1,
},
'SiH2_1A1.ezfio' : {
'Orthonormal' : 1,
},
'SiH2_3B1.ezfio' : {
'Orthonormal' : 1,
},
'SiH2-.ezfio' : {
'Orthonormal' : 1,
},
'SiH3-.ezfio' : {
'Orthonormal' : 1,
},
'SiH3.ezfio' : {
'Orthonormal' : 1,
},
'SiH4.ezfio' : {
'Orthonormal' : 1,
},
'SiH4+.ezfio' : {
'Orthonormal' : 1,
},
'SiH-.ezfio' : {
'Orthonormal' : 1,
},
'SiH.ezfio' : {
'Orthonormal' : 1,
},
'SiS.ezfio' : {
'Orthonormal' : 1,
},
}

12
src/MonoInts/tests/mo_pot.irp.f Normal file
View File

@ -0,0 +1,12 @@
program test mo_pot
implicit none
integer :: i, j
do i=1,mo_tot_num
do j=1,i
if (abs(mo_nucl_elec_integral(i,j)) > 1.d-9) then
print '(I4,A,I4,2x,A, 2x,e16.8)' , i,',',j, ' : ', mo_nucl_elec_integral(i,j)
endif
enddo
enddo
end

44923
src/MonoInts/tests/mo_pot.ref Normal file
View File

File diff suppressed because it is too large Load Diff

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix

16
src/Utils/util.irp.f
View File

@ -1,3 +1,18 @@
BEGIN_PROVIDER [ logical, all_utils ]
implicit none
BEGIN_DOC
! Dummy provider to provide all utils
END_DOC
! Do not move this : it greps itself
BEGIN_SHELL [ /bin/bash ]
for i in $(grep "BEGIN_PROVIDER" $QPACKAGE_ROOT/src/Utils/*.irp.f | cut -d ',' -f 2 | cut -d ']' -f 1 | tail --lines=+3 )
do
echo PROVIDE $i
done
END_SHELL
END_PROVIDER
double precision function binom_func(i,j)
implicit none
BEGIN_DOC
@ -162,7 +177,6 @@ BEGIN_PROVIDER [ double precision, inv_int, (128) ]
do i=1,size(inv_int)
inv_int(i) = 1.d0/dble(i)
enddo
END_PROVIDER
subroutine wall_time(t)