10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-30 00:44:28 +02:00

Merge pull request #79 from scemama/master

qp_install cleaning
This commit is contained in:
Thomas Applencourt 2015-06-22 11:20:10 +02:00
commit 76823bf788
14 changed files with 53 additions and 701 deletions

View File

@ -11,7 +11,7 @@
#
[COMMON]
FC : gfortran -ffree-line-length-none -I .
LAPACK_LIB : -lblas -llapack
LAPACK_LIB : -llapack -lblas
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32

View File

@ -1,226 +1,3 @@
BEGIN_PROVIDER [ double precision, aux_pseudo_integral, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
implicit none
BEGIN_DOC
! Pseudo-potential
END_DOC
if (do_pseudo) then
! aux_pseudo_integral = aux_pseudo_integral_local + aux_pseudo_integral_non_local
! aux_pseudo_integral = aux_pseudo_integral_local
aux_pseudo_integral = aux_pseudo_integral_non_local
else
aux_pseudo_integral = 0.d0
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, aux_pseudo_integral_local, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
implicit none
BEGIN_DOC
! Local pseudo-potential
END_DOC
double precision :: alpha, beta, gama, delta
integer :: 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 :: Vloc, Vpseudo
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
integer :: thread_num
aux_pseudo_integral_local = 0.d0
!! Dump array
integer, allocatable :: n_k_dump(:)
double precision, allocatable :: v_k_dump(:), dz_k_dump(:)
allocate(n_k_dump(1:pseudo_klocmax), v_k_dump(1:pseudo_klocmax), dz_k_dump(1:pseudo_klocmax))
! _
! / _. | _ |
! \_ (_| | (_ |_| |
!
print*, 'Providing the nuclear electron pseudo integrals '
call wall_time(wall_1)
call cpu_time(cpu_1)
!$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 v_k_dump,n_k_dump, dz_k_dump, &
!$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (aux_basis_num_sqrt,aux_basis_prim_num,aux_basis_expo_transp,aux_basis_power,aux_basis_nucl,nucl_coord,aux_basis_coef_transp, &
!$OMP aux_pseudo_integral_local,nucl_num,nucl_charge, &
!$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k,pseudo_n_k, pseudo_dz_k, &
!$OMP wall_1)
!$OMP DO SCHEDULE (guided)
do j = 1, aux_basis_num_sqrt
num_A = aux_basis_nucl(j)
power_A(1:3)= aux_basis_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, aux_basis_num_sqrt
num_B = aux_basis_nucl(i)
power_B(1:3)= aux_basis_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,aux_basis_prim_num(j)
alpha = aux_basis_expo_transp(l,j)
do m=1,aux_basis_prim_num(i)
beta = aux_basis_expo_transp(m,i)
double precision :: c
c = 0.d0
do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
v_k_dump = pseudo_v_k(k,1:pseudo_klocmax)
n_k_dump = pseudo_n_k(k,1:pseudo_klocmax)
dz_k_dump = pseudo_dz_k(k,1:pseudo_klocmax)
c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump, &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
aux_pseudo_integral_local(i,j) = aux_pseudo_integral_local(i,j) + &
aux_basis_coef_transp(l,j)*aux_basis_coef_transp(m,i)*c
enddo
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(j)/float(aux_basis_num_sqrt), '% in ', &
wall_2-wall_1, 's'
endif
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
deallocate(n_k_dump,v_k_dump, dz_k_dump)
END_PROVIDER
BEGIN_PROVIDER [ double precision, aux_pseudo_integral_non_local, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
implicit none
BEGIN_DOC
! Local pseudo-potential
END_DOC
double precision :: alpha, beta, gama, delta
integer :: 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 :: Vloc, Vpseudo
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
integer :: thread_num
aux_pseudo_integral_non_local = 0.d0
!! Dump array
integer, allocatable :: n_kl_dump(:,:)
double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:)
allocate(n_kl_dump(pseudo_kmax,0:pseudo_lmax), v_kl_dump(pseudo_kmax,0:pseudo_lmax), dz_kl_dump(pseudo_kmax,0:pseudo_lmax))
! _
! / _. | _ |
! \_ (_| | (_ |_| |
!
print*, 'Providing the nuclear electron pseudo integrals '
call wall_time(wall_1)
call cpu_time(cpu_1)
!$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 n_kl_dump, v_kl_dump, dz_kl_dump, &
!$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (aux_basis_num_sqrt,aux_basis_prim_num,aux_basis_expo_transp,aux_basis_power,aux_basis_nucl,nucl_coord,aux_basis_coef_transp, &
!$OMP aux_pseudo_integral_non_local,nucl_num,nucl_charge, &
!$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl, &
!$OMP wall_1)
!$OMP DO SCHEDULE (guided)
do j = 1, aux_basis_num_sqrt
num_A = aux_basis_nucl(j)
power_A(1:3)= aux_basis_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, aux_basis_num_sqrt
num_B = aux_basis_nucl(i)
power_B(1:3)= aux_basis_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,aux_basis_prim_num(j)
alpha = aux_basis_expo_transp(l,j)
do m=1,aux_basis_prim_num(i)
beta = aux_basis_expo_transp(m,i)
double precision :: c
c = 0.d0
do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
n_kl_dump = pseudo_n_kl(k,1:pseudo_kmax,0:pseudo_lmax)
v_kl_dump = pseudo_v_kl(k,1:pseudo_kmax,0:pseudo_lmax)
dz_kl_dump = pseudo_dz_kl(k,1:pseudo_kmax,0:pseudo_lmax)
c = c + Vpseudo(pseudo_lmax,pseudo_kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
aux_pseudo_integral_non_local(i,j) = aux_pseudo_integral_non_local(i,j) + &
aux_basis_coef_transp(l,j)*aux_basis_coef_transp(m,i)*c
enddo
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(j)/float(aux_basis_num_sqrt), '% in ', &
wall_2-wall_1, 's'
endif
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ]
implicit none
@ -240,9 +17,11 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l
integer :: n_a(3)
double precision :: a(3), c(3), g_a
integer :: i,j,k,l,m,n,p
double precision :: r(pseudo_grid_size), dr, Ulc
double precision :: dr, Ulc
double precision :: y
double precision, allocatable :: r(:)
allocate (r(pseudo_grid_size))
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
do j=2,pseudo_grid_size
@ -269,6 +48,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l
enddo
enddo
enddo
deallocate(r)
END_PROVIDER
@ -291,8 +71,11 @@ BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l
integer :: n_a(3)
double precision :: a(3), c(3), g_a
integer :: i,j,k,l,m,n,p
double precision :: r(pseudo_grid_size), dr, Ulc
double precision :: dr, Ulc
double precision :: y
double precision, allocatable :: r(:)
allocate (r(pseudo_grid_size))
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
@ -315,6 +98,7 @@ BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l
enddo
enddo
enddo
deallocate(r)
END_PROVIDER

View File

@ -3,169 +3,6 @@ subroutine write_pseudopotential
BEGIN_DOC
! Write the pseudo_potential into the EZFIO file
END_DOC
! call ezfio_set_pseudo_pseudo_matrix(pseudo_matrix)
! call ezfio_set_pseudo_ao_pseudo_grid(ao_pseudo_grid)
call ezfio_set_pseudo_mo_pseudo_grid(mo_pseudo_grid)
end
BEGIN_PROVIDER [ double precision, pseudo_matrix, (aux_basis_num_sqrt,aux_basis_num_sqrt) ]
implicit none
BEGIN_DOC
! Pseudo-potential expressed in the basis of ao products
END_DOC
integer :: i,j,k,l
integer :: info, m,n, lwork, lda, ldu, ldvt
integer, allocatable :: iwork(:)
character :: jobz
double precision, allocatable :: a(:,:),work(:)
double precision,allocatable :: U(:,:)
double precision,allocatable :: Vt(:,:)
double precision,allocatable :: S(:), B(:)
jobz = 'A'
m = aux_basis_num
n = aux_basis_num
lda = size(aux_basis_overlap_matrix,1)
ldu = lda
ldvt = lda
lwork = -1
! allocate (A(lda,n), U(ldu,n), Vt(ldvt,n), S(n), work(1), b(n), iwork(8*n))
allocate (A(lda,n), U(ldu,n), Vt(ldvt,n), S(n), work(1), b(n),iwork(1))
work(1) = 1
do i=1,n
do j=1,n
A(i,j) = aux_basis_overlap_matrix(i,j)
enddo
enddo
! call dgesdd(jobz, m, n, A, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
call dgesvd(jobz, jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
lwork = int(work(1))
deallocate(work)
print *, 'Fitting pseudo-potentials'
allocate(work(lwork))
! call dgesdd(jobz, m, n, A, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
call dgesvd(jobz, jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
deallocate(work)
do i=1,n
print *, i, s(i)
enddo
do k=1,n
if (s(k) < 1.d-1) then
s(k) = 0.d0
else
s(k) = 1.d0/s(k)
endif
do m=1,n
Vt(m,k) = S(m) * Vt(m,k)
enddo
enddo
call dgemm('N','N',n,n,n,1.d0,U,lda,Vt,ldvt,0.d0,A,lda)
! do k=1,n
! do l=1,n
! A(k,l) = 0.d0
! do m=1,n
! A(k,l) = A(k,l) + U(k,m) * Vt(m,l)
! enddo
! enddo
do k=1,n
i = aux_basis_idx(1,k)
j = aux_basis_idx(2,k)
b(k) = aux_pseudo_integral(i,j)
enddo
do k=1,n
S(k) = 0.d0
enddo
do l=1,n
do k=1,n
S(k) = S(k) + A(k,l) * b(l)
enddo
enddo
do k=1,aux_basis_num
i = aux_basis_idx(1,k)
j = aux_basis_idx(2,k)
pseudo_matrix(i,j) = S(k)
pseudo_matrix(j,i) = S(k)
enddo
deallocate(a,b,s,iwork,u,vt)
print *, 'Done'
if (info /= 0) then
print *, info
stop 'pseudo fit failed'
endif
END_PROVIDER
!BEGIN_PROVIDER [ double precision, pseudo_matrix, (ao_num,ao_num) ]
! implicit none
! BEGIN_DOC
! ! Pseudo-potential expressed in the basis of ao products
! END_DOC
!
! integer :: i,j,k
! integer :: info, n, lwork, lda, ldb, nrhs
! character :: uplo
! integer, allocatable :: ipiv(:)
! double precision, allocatable :: a(:,:),work(:), b(:)
!
! uplo = 'L'
! n = aux_basis_num
! nrhs = 1
! lda = size(aux_basis_overlap_matrix,1)
! ldb = n
! lwork = -1
!
! print *, 'Fitting pseudo-potentials'
! allocate(work(1),a(lda,n),ipiv(n),b(n))
! work(1) = 1
! do i=1,n
! do j=1,n
! a(i,j) = aux_basis_overlap_matrix(i,j)
! enddo
! enddo
!
! do k=1,n
! i = aux_basis_idx(1,k)
! j = aux_basis_idx(2,k)
! b(k) = ao_pseudo_integral(i,j)
! enddo
! call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info )
! lwork = int(work(1))
! deallocate(work)
!
! allocate(work(lwork))
! call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info )
! deallocate(work,ipiv)
! do k=1,aux_basis_num
! i = aux_basis_idx(1,k)
! j = aux_basis_idx(2,k)
! pseudo_matrix(i,j) = b(k)
! pseudo_matrix(j,i) = b(k)
! enddo
! deallocate(a,b)
!
!print *, 'Done'
! if (info /= 0) then
! print *, info
! stop 'pseudo fit failed'
! endif
!END_PROVIDER

View File

@ -18,7 +18,9 @@ p = re.compile(ur'-I IRPF90_temp/\S*\s+')
mod = re.compile(ur'module\s+(?P<mod>\S+).+end\s?module\s+(?P=mod)?',
re.MULTILINE | re.IGNORECASE)
TMPDIR = "/tmp/qp_compiler/"
tmpdir_root = os.environ.get("TMPDIR",failobj="/dev/shm")
TMPDIR = os.path.join(tmpdir_root,os.environ["USER"],"qp_compiler")
def return_filename_to_cache(command):
@ -67,7 +69,7 @@ def cache_utility(command):
# Create temp directory
try:
os.mkdir("/tmp/qp_compiler/")
os.makedirs(TMPDIR)
except OSError:
pass

View File

@ -19,7 +19,11 @@ try:
from read_compilation_cfg import get_compilation_option
from docopt import docopt
except ImportError:
print "source .quantum_package.rc"
f = os.path.realpath(os.path.join(os.path.dirname(__file__),"..","..","quantum_package.rc"))
print """
Error:
source %s
"""%f
sys.exit(1)
header = r"""#
@ -253,7 +257,7 @@ def ninja_ezfio_rule():
install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio.a")
l_cmd = ["cd {0}".format(QP_EZFIO)] + l_flag
l_cmd += ["ninja && ln -f {0} {1}".format(install_lib_ezfio, EZFIO_LIB)]
l_cmd += ["ninja && ln -sf {0} {1}".format(install_lib_ezfio, EZFIO_LIB)]
l_string = ["rule build_ezfio",
" command = {0}".format(" ; ".join(l_cmd)),
@ -677,7 +681,8 @@ def ninja_dot_tree_build(path_module):
def create_build_ninja_module(path_module):
l_string = ["rule update_build_ninja_root",
" command = qp_create_ninja.py update", ""]
" command = {0} update".format(__file__),
""]
l_string += ["rule make_local_binaries",
" command = ninja -f {0} module_{1}".format(
@ -709,7 +714,8 @@ def create_build_ninja_module(path_module):
def create_build_ninja_global():
l_string = ["rule update_build_ninja_root",
" command = qp_create_ninja.py update", ""]
" command = {0} update".format(__file__),
""]
l_string += ["rule make_all_binaries",
" command = ninja -f {0}".format(ROOT_BUILD_NINJA),

View File

@ -1,10 +1,12 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Usage: qp_install_module.py list (--installed|--avalaible-local|--avalaible-remote)
qp_install_module.py install <name>...
Usage:
qp_install_module.py create -n <name> [<children_module>...]
qp_install_module.py download -n <name> [<path_folder>...]
qp_install_module.py install <name>...
qp_install_module.py list (--installed|--avalaible-local|--avalaible-remote)
qp_install_module.py uninstall <name>...
Options:
@ -69,7 +71,7 @@ if __name__ == '__main__':
m_instance = ModuleHandler(l_repository)
for module in m_instance.l_module:
for module in sorted(m_instance.l_module):
print "* {0}".format(module)
elif arguments["create"]:
@ -143,6 +145,29 @@ if __name__ == '__main__':
for module_to_cp in l_module_to_cp:
src = os.path.join(qp_root_plugin, module_to_cp)
des = os.path.join(qp_root_src, module_to_cp)
os.symlink(src, des)
try:
os.symlink(src, des)
except OSError:
print "Your src directory is broken. Please remove %s"%des
raise
print "Done"
print "You can now compile as usual"
elif arguments["uninstall"]:
d_local = get_dict_child([qp_root_src])
l_name = arguments["<name>"]
l_failed = [ name for name in l_name if name not in d_local ]
if l_failed:
print "Modules not installed:"
for name in sorted(l_failed):
print "* %s"%name
sys.exit(1)
else:
def unlink(x):
try:
os.unlink(os.path.join(qp_root_src,x))
except OSError:
print "%s is a core module which can not be renmoved"%x
map(unlink,l_name)

View File

@ -1,5 +0,0 @@
IRPF90_temp/
IRPF90_man/
irpf90.make
irpf90_entities
tags

View File

@ -1 +0,0 @@
AOs Pseudo

View File

@ -1,82 +0,0 @@
=================
DensityFit Module
=================
In this module, the basis of all the products of atomic orbitals is built.
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. by the `update_README.py` script.
`aux_basis_coef <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L94>`_
Exponents and coefficients of the auxiliary basis
`aux_basis_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L37>`_
Exponents of the auxiliary basis
`aux_basis_expo <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L93>`_
Exponents and coefficients of the auxiliary basis
`aux_basis_expo_transp <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L36>`_
Exponents of the auxiliary basis
`aux_basis_four_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/overlap.irp.f#L1>`_
\int \chi_i(r) \chi_j(r) \chi_k(r) \chi_l(r) dr
`aux_basis_idx <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L20>`_
aux_basis_idx(k) -> i,j
`aux_basis_nucl <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L40>`_
Exponents of the auxiliary basis
`aux_basis_num <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L2>`_
Number of auxiliary basis functions
`aux_basis_num_8 <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L3>`_
Number of auxiliary basis functions
`aux_basis_num_sqrt <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L1>`_
Number of auxiliary basis functions
`aux_basis_overlap_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L69>`_
Auxiliary basis set
`aux_basis_power <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L39>`_
Exponents of the auxiliary basis
`aux_basis_prim_num <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L38>`_
Exponents of the auxiliary basis
`aux_basis_prim_num_max <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L111>`_
= ao_prim_num_max
`save_aux_basis <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L120>`_
Undocumented
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. by the `update_README.py` script.
.. image:: tree_dependency.png
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `Pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo>`_

View File

@ -1,12 +0,0 @@
aux_basis
aux_basis_num integer
aux_basis_num_sqrt integer
aux_basis_idx integer (2,aux_basis_aux_basis_num)
aux_basis_prim_num integer (aux_basis_aux_basis_num_sqrt)
aux_basis_nucl integer (aux_basis_aux_basis_num_sqrt)
aux_basis_power integer (aux_basis_aux_basis_num_sqrt,3)
aux_basis_prim_num_max integer = maxval(aux_basis_aux_basis_prim_num)
aux_basis_coef double precision (aux_basis_aux_basis_num_sqrt,aux_basis_aux_basis_prim_num_max)
aux_basis_expo double precision (aux_basis_aux_basis_num_sqrt,aux_basis_aux_basis_prim_num_max)

View File

@ -1,130 +0,0 @@
BEGIN_PROVIDER [ integer, aux_basis_num_sqrt ]
&BEGIN_PROVIDER [ integer, aux_basis_num ]
&BEGIN_PROVIDER [ integer, aux_basis_num_8 ]
implicit none
BEGIN_DOC
! Number of auxiliary basis functions
END_DOC
integer :: align_double
if (do_pseudo) then
! aux_basis_num_sqrt = ao_num + ao_pseudo_num
aux_basis_num_sqrt = ao_num
else
endif
aux_basis_num = aux_basis_num_sqrt * (aux_basis_num_sqrt+1)/2
aux_basis_num_8 = align_double(aux_basis_num)
END_PROVIDER
BEGIN_PROVIDER [ integer, aux_basis_idx, (2,aux_basis_num) ]
implicit none
BEGIN_DOC
! aux_basis_idx(k) -> i,j
END_DOC
integer :: i,j,k
k=0
do j=1,aux_basis_num_sqrt
do i=1,j
k = k+1
aux_basis_idx(1,k) = i
aux_basis_idx(2,k) = j
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, aux_basis_expo_transp, (ao_prim_num_max_align,aux_basis_num_sqrt) ]
&BEGIN_PROVIDER [ double precision, aux_basis_coef_transp, (ao_prim_num_max_align,aux_basis_num_sqrt) ]
&BEGIN_PROVIDER [ integer, aux_basis_prim_num, (aux_basis_num_sqrt) ]
&BEGIN_PROVIDER [ integer, aux_basis_power, (aux_basis_num_sqrt,3) ]
&BEGIN_PROVIDER [ integer, aux_basis_nucl, (aux_basis_num_sqrt) ]
implicit none
BEGIN_DOC
! Exponents of the auxiliary basis
END_DOC
integer :: i,j
do j=1,ao_num
do i=1,ao_prim_num_max
aux_basis_expo_transp(i,j) = ao_expo_ordered_transp(i,j)
aux_basis_coef_transp(i,j) = ao_coef_normalized_ordered_transp(i,j)
enddo
enddo
do i=1,ao_num
aux_basis_prim_num(i) = ao_prim_num(i)
aux_basis_nucl(i) = ao_nucl(i)
aux_basis_power(i,1:3) = ao_power(i,1:3)
enddo
! do j=1,ao_pseudo_num
! aux_basis_expo_transp(1,ao_num+j) = 0.5d0*pseudo_ao_expo(j)
! aux_basis_coef_transp(1,ao_num+j) = 1.d0
! aux_basis_power(ao_num+j,1:3) = 0
! aux_basis_prim_num(ao_num+j) = 1
! aux_basis_nucl(ao_num+j) = pseudo_ao_nucl(j)
! enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, aux_basis_overlap_matrix, (aux_basis_num_8,aux_basis_num) ]
implicit none
BEGIN_DOC
! Auxiliary basis set
END_DOC
integer :: m,n,i,j,k,l
double precision :: aux_basis_four_overlap
aux_basis_overlap_matrix(1,1) = aux_basis_four_overlap(1,1,1,1)
!$OMP PARALLEL DO PRIVATE(i,j,k,l,m,n) SCHEDULE(GUIDED)
do m=1,aux_basis_num
i = aux_basis_idx(1,m)
j = aux_basis_idx(2,m)
do n=1,m
k = aux_basis_idx(1,n)
l = aux_basis_idx(2,n)
aux_basis_overlap_matrix(m,n) = aux_basis_four_overlap(i,j,k,l)
aux_basis_overlap_matrix(n,m) = aux_basis_overlap_matrix(m,n)
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, aux_basis_expo, (aux_basis_num_sqrt,aux_basis_prim_num_max) ]
&BEGIN_PROVIDER [ double precision, aux_basis_coef, (aux_basis_num_sqrt,aux_basis_prim_num_max) ]
implicit none
BEGIN_DOC
! Exponents and coefficients of the auxiliary basis
END_DOC
integer :: i,j
aux_basis_expo = 0.d0
aux_basis_coef = 0.d0
do j=1,aux_basis_num_sqrt
do i=1,aux_basis_prim_num(j)
aux_basis_expo(j,i) = aux_basis_expo_transp(i,j)
aux_basis_coef(j,i) = aux_basis_coef_transp(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, aux_basis_prim_num_max ]
implicit none
BEGIN_DOC
! = ao_prim_num_max
END_DOC
aux_basis_prim_num_max = ao_prim_num_max
END_PROVIDER
subroutine save_aux_basis
implicit none
call ezfio_set_aux_basis_aux_basis_num(aux_basis_num)
call ezfio_set_aux_basis_aux_basis_num_sqrt(aux_basis_num_sqrt)
call ezfio_set_aux_basis_aux_basis_idx(aux_basis_idx)
call ezfio_set_aux_basis_aux_basis_prim_num(aux_basis_prim_num)
call ezfio_set_aux_basis_aux_basis_nucl(aux_basis_nucl)
call ezfio_set_aux_basis_aux_basis_power(aux_basis_power)
call ezfio_set_aux_basis_aux_basis_coef(aux_basis_coef)
call ezfio_set_aux_basis_aux_basis_expo(aux_basis_expo)
end

View File

@ -1,66 +0,0 @@
double precision function aux_basis_four_overlap(i,j,k,l)
implicit none
BEGIN_DOC
! \int \chi_i(r) \chi_j(r) \chi_k(r) \chi_l(r) dr
END_DOC
integer,intent(in) :: i,j,k,l
integer :: p,q,r,s
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
double precision :: overlap_x,overlap_y,overlap_z, overlap
include 'Utils/constants.include.F'
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
integer :: iorder_p(3), iorder_q(3)
dim1 = n_pt_max_integrals
num_i = aux_basis_nucl(i)
num_j = aux_basis_nucl(j)
num_k = aux_basis_nucl(k)
num_l = aux_basis_nucl(l)
aux_basis_four_overlap = 0.d0
do p = 1, 3
I_power(p) = aux_basis_power(i,p)
J_power(p) = aux_basis_power(j,p)
K_power(p) = aux_basis_power(k,p)
L_power(p) = aux_basis_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = nucl_coord(num_j,p)
K_center(p) = nucl_coord(num_k,p)
L_center(p) = nucl_coord(num_l,p)
enddo
do p = 1, aux_basis_prim_num(i)
double precision :: coef1
coef1 = aux_basis_coef_transp(p,i)
do q = 1, aux_basis_prim_num(j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
aux_basis_expo_transp(p,i),aux_basis_expo_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
double precision :: coef2
coef2 = coef1*aux_basis_coef_transp(q,j)*fact_p
do r = 1, aux_basis_prim_num(k)
double precision :: coef3
coef3 = coef2*aux_basis_coef_transp(r,k)
do s = 1, aux_basis_prim_num(l)
double precision :: general_primitive_integral
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q, &
aux_basis_expo_transp(r,k),aux_basis_expo_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
double precision :: coef4
coef4 = coef3*aux_basis_coef_transp(s,l)*fact_q
call overlap_gaussian_xyz(P_center,Q_center,pp,qq,iorder_p,iorder_q,overlap_x,overlap_y,overlap_z,overlap,dim1)
aux_basis_four_overlap += coef4 * overlap
enddo ! s
enddo ! r
enddo ! q
enddo ! p
end
! TODO : Schwartz acceleration

Binary file not shown.

Before

Width:  |  Height:  |  Size: 24 KiB

View File

@ -79,10 +79,4 @@ doc: QMC grid
interface: ezfio
size: (ao_basis.ao_num,-pseudo.pseudo_lmax:pseudo.pseudo_lmax,0:pseudo.pseudo_lmax,nuclei.nucl_num,pseudo.pseudo_grid_size)
[pseudo_matrix]
type: double precision
doc: QMC grid
interface: ezfio
size: (aux_basis.aux_basis_num_sqrt,aux_basis.aux_basis_num_sqrt)