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

Try to merge...

Browse Source
This commit is contained in:
Thomas Applencourt 2015-06-03 14:03:11 +02:00
commit e99034d3f9
47 changed files with 1713 additions and 127 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ocaml/qp_create_ezfio_from_xyz.ml
View File

@ -14,7 +14,7 @@ let spec =
+> flag "m" (optional_with_default 1 int) +> flag "m" (optional_with_default 1 int)
~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1." ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1."
+> flag "p" no_arg +> flag "p" no_arg
~doc:"Using pseudo." ~doc:"Using pseudopotentials"
+> anon ("xyz_file" %: string) +> anon ("xyz_file" %: string)
;; ;;

62
scripts/cache_compile.py Executable file
View File

@ -0,0 +1,62 @@
#!/usr/bin/env python
import os
import sys
import shelve
import hashlib
import re
r = re.compile(ur'-c\s+(\S+\.[fF]90)\s+-o\s+(\S+\.o)')
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/"
def main():
# Create temp directory
if "qp_compiler" not in os.listdir("/tmp"):
os.mkdir("/tmp/qp_compiler/")
line = sys.argv[1:]
command = " ".join(line)
command_clean = p.sub('',command)
try:
match = r.search(command_clean)
input = match.group(1)
output = match.group(2)
except:
os.system(command)
return
m = hashlib.md5()
# Fread : read input
with open(input,'r') as file:
fread = file.read()
m.update( " ".join( [ command, fread ] ))
# Md5 Key containing command + content of Fread
key = TMPDIR+m.hexdigest()
try:
# Try to return the content of the .o file
with open(key,'r') as file:
result = file.read()
except IOError:
# Compile the file -> .o
os.system(command)
# Read the .o
with open(output,'r') as file:
result = file.read()
# Copy the .o in database
if not mod.search(fread.replace('\n',' ')):
with open(key,'w') as file:
file.write(result)
else:
print input+' -> module'
else:
# Write the .o file
with open(output,'w') as file:
file.write(result)
if __name__ == '__main__':
main()

5
scripts/ezfio_interface/ei_handler.py
View File

@ -365,7 +365,10 @@ def create_ezfio_stuff(dict_ezfio_cfg, config_or_default="config"):
except ValueError: except ValueError:
a_size_raw.append(dim) a_size_raw.append(dim)
else: else:
a_size_raw.append("{0}+{1}+1".format(begin, end)) if begin[0] == '-':
a_size_raw.append("{0}+{1}+1".format(end, begin[1:]))
else:
a_size_raw.append("{0}-{1}+1".format(end, begin))
size_raw = ",".join(a_size_raw) size_raw = ",".join(a_size_raw)

2
scripts/get_basis.sh
View File

@ -51,5 +51,7 @@ then
${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" ${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}"
else else
${EMSL_API_ROOT}/EMSL_api.py get_basis_data --save --path="${tmpfile}" --basis="${basis}" --db_path="${EMSL_API_ROOT}/db/Pseudo.db" ${EMSL_API_ROOT}/EMSL_api.py get_basis_data --save --path="${tmpfile}" --basis="${basis}" --db_path="${EMSL_API_ROOT}/db/Pseudo.db"
# echo ${EMSL_API_ROOT}/EMSL_api.py get_basis_data --save --path="${tmpfile}" --basis="${basis}" --db_path="${EMSL_API_ROOT}/db/Pseudo.db" 1>&2
# echo $PWD/BASIS
fi fi

2
scripts/install/install_docopt.sh
View File

@ -18,4 +18,4 @@ cd ${QPACKAGE_ROOT}
rm -f -- scripts/${DOCOPT}{,c} rm -f -- scripts/${DOCOPT}{,c}
${QPACKAGE_ROOT}/scripts/install/fetch_from_web.py ${DOCOPT_URL} ${DOCOPT} ${QPACKAGE_ROOT}/scripts/install/fetch_from_web.py ${DOCOPT_URL} ${DOCOPT}
mv ${DOCOPT} scripts/utility/${DOCOPT} mv ${DOCOPT} scripts/utility/${DOCOPT}

1
scripts/install/install_m4.sh
View File

@ -16,6 +16,7 @@ fi
cd ${QPACKAGE_ROOT} cd ${QPACKAGE_ROOT}
rm -f l${QPACKAGE_ROOT}/bin/m4
if [[ -z ${M4} ]] if [[ -z ${M4} ]]
then then
rm -f -- bin/m4 rm -f -- bin/m4

6
scripts/install/install_ninja.sh
View File

@ -1,7 +1,7 @@
#!/bin/bash #!/bin/bash
# #
# Installs m4 for ocaml # Installs the ninja build system
# Thu Oct 23 22:02:08 CEST 2014 # Thu May 28 13:21:16 CEST 2015
BASE="ninja" BASE="ninja"
URL="https://github.com/martine/ninja/archive/master.tar.gz" URL="https://github.com/martine/ninja/archive/master.tar.gz"
@ -21,5 +21,3 @@ rm -rf ${BASE}
mv ${BASE}-master ${BASE} mv ${BASE}-master ${BASE}
cd ${BASE} cd ${BASE}
./configure.py --bootstrap ./configure.py --bootstrap

46
scripts/module/create_Makefile_depend.sh Executable file
View File

@ -0,0 +1,46 @@
#!/bin/bash
#
# This script is automatically invoked by Makefiles and should not be invoked
# by users.
# Creates the Makefile.depend file. This file contains all the external source
# files included by including other modules.
# Thu Apr 3 01:44:09 CEST 2014
if [[ -z ${QPACKAGE_ROOT} ]]
then
print "The QPACKAGE_ROOT environment variable is not set."
print "Please reload the quantum_package.rc file."
exit -1
fi
source ${QPACKAGE_ROOT}/scripts/qp_include.sh
check_current_dir_is_module
SRC=""
OBJ=""
DEPS="$NEEDED_MODULES"
for M in ${DEPS}
do
# X is the list of external source files
X=$(grep '^SRC=' "${QPACKAGE_ROOT}/src/${M}/Makefile" 2>/dev/null |cut -d '=' -f 2)
for f in ${X}
do
SRC+=" ${M}/${f}"
done
X=$(grep '^OBJ=' "${QPACKAGE_ROOT}/src/${M}/Makefile" 2>/dev/null |cut -d '=' -f 2)
for f in ${X}
do
OBJ+=" IRPF90_temp/${M}/${f/IRPF90_temp//}"
done
done
# Create the Makefile.depend
cat << EOF > Makefile.depend
# This file was created by the $0 script. Do not modify it by hand.
SRC+=${SRC}
OBJ+=${OBJ}
EOF

2
scripts/pseudo/put_pseudo_in_ezfio.py
View File

@ -186,7 +186,7 @@ def add_zero(array, size, type):
def make_it_square(matrix, dim, type=float): def make_it_square(matrix, dim, type=float):
""" """
matix the matrix to squate matix the matrix to square
dim array [lmax, kmax] dim array [lmax, kmax]
type the null value you want type the null value you want
[[[28.59107316], [19.37583724]], [[50.25646328]]] [[[28.59107316], [19.37583724]], [[50.25646328]]]

2
setup_environment.sh
View File

@ -29,7 +29,7 @@ EOF
source quantum_package.rc source quantum_package.rc
mkdir -p install_logs
echo "${BLUE}===== Installing IRPF90 ===== ${BLACK}" echo "${BLUE}===== Installing IRPF90 ===== ${BLACK}"
${QPACKAGE_ROOT}/scripts/install/install_irpf90.sh | tee ${QPACKAGE_ROOT}/install_logs/install_irpf90.log ${QPACKAGE_ROOT}/scripts/install/install_irpf90.sh | tee ${QPACKAGE_ROOT}/install_logs/install_irpf90.log
if [[ ! -d ${QPACKAGE_ROOT}/irpf90 ]] || [[ ! -x ${QPACKAGE_ROOT}/bin/irpf90 ]] || [[ ! -x ${QPACKAGE_ROOT}/bin/irpman ]] if [[ ! -d ${QPACKAGE_ROOT}/irpf90 ]] || [[ ! -x ${QPACKAGE_ROOT}/bin/irpf90 ]] || [[ ! -x ${QPACKAGE_ROOT}/bin/irpman ]]

2
src/Bitmask/bitmasks_module.f90
View File

@ -8,4 +8,4 @@ module bitmasks
integer, parameter :: d_part2 = 4 integer, parameter :: d_part2 = 4
integer, parameter :: s_hole = 5 integer, parameter :: s_hole = 5
integer, parameter :: s_part = 6 integer, parameter :: s_part = 6
end module end module bitmasks

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

6
src/DensityFit/Makefile Normal file
View File

@ -0,0 +1,6 @@
# 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/DensityFit/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
AOs Pseudo

70
src/DensityFit/README.rst Normal file
View File

@ -0,0 +1,70 @@
=================
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
.. NEEDED_MODULES file.
`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_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
`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
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
.. image:: tree_dependancy.png
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `Pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo>`_

12
src/DensityFit/aux_basis.ezfio_config Normal file
View File

@ -0,0 +1,12 @@
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)

130
src/DensityFit/aux_basis.irp.f Normal file
View File

@ -0,0 +1,130 @@
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

66
src/DensityFit/overlap.irp.f Normal file
View File

@ -0,0 +1,66 @@
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 'include/constants.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

6
src/Determinants/README.rst
View File

@ -532,12 +532,6 @@ Documentation
`save_casino <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_casino.irp.f#L1>`_ `save_casino <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_casino.irp.f#L1>`_
Undocumented Undocumented
`save_dets_qmcchem <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L1>`_
Undocumented
`save_for_qmc <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L46>`_
Undocumented
`save_natorb <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_natorb.irp.f#L1>`_ `save_natorb <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_natorb.irp.f#L1>`_
Undocumented Undocumented

51
src/Determinants/save_for_qmcchem.irp.f
View File

@ -1,51 +0,0 @@
subroutine save_dets_qmcchem
use bitmasks
implicit none
character :: c(mo_tot_num)
integer :: i,k
integer, allocatable :: occ(:,:,:), occ_tmp(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: occ, occ_tmp
call ezfio_set_determinants_det_num(N_det)
call ezfio_set_determinants_det_coef(psi_coef_sorted(1,1))
allocate (occ(elec_alpha_num,N_det,2))
! OMP PARALLEL DEFAULT(NONE) &
! OMP PRIVATE(occ_tmp,i,k)&
! OMP SHARED(N_det,psi_det_sorted,elec_alpha_num, &
! OMP occ,elec_beta_num,N_int)
allocate (occ_tmp(N_int*bit_kind_size,2))
occ_tmp = 0
! OMP DO
do i=1,N_det
call bitstring_to_list(psi_det_sorted(1,1,i), occ_tmp(1,1), elec_alpha_num, N_int )
call bitstring_to_list(psi_det_sorted(1,2,i), occ_tmp(1,2), elec_beta_num, N_int )
do k=1,elec_alpha_num
occ(k,i,1) = occ_tmp(k,1)
occ(k,i,2) = occ_tmp(k,2)
enddo
enddo
! OMP END DO
deallocate(occ_tmp)
! OMP END PARALLEL
call ezfio_set_determinants_det_occ(occ)
call write_int(output_determinants,N_det,'Determinants saved for QMC')
deallocate(occ)
open(unit=31,file=trim(ezfio_filename)//'/mo_basis/mo_classif')
write(31,'(I1)') 1
write(31,*) mo_tot_num
do i=1,mo_tot_num
write(31,'(A)') 'a'
enddo
close(31)
call system('gzip -f '//trim(ezfio_filename)//'/mo_basis/mo_classif')
end
program save_for_qmc
read_wf = .True.
TOUCH read_wf
! call save_dets_qmcchem
call write_spindeterminants
end

3
src/FCIdump/README.rst
View File

@ -10,9 +10,6 @@ Documentation
.. Do not edit this section. It was auto-generated from the .. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES_CHILDREN file by the `update_README.py` script. .. NEEDED_MODULES_CHILDREN file by the `update_README.py` script.
`fcidump <http://github.com/LCPQ/quantum_package/tree/master/src/FCIdump/fcidump.irp.f#L1>`_
Undocumented
Needed Modules Needed Modules

4
src/Integrals_Monoelec/README.rst
View File

@ -97,10 +97,10 @@ Documentation
`ao_pseudo_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L1>`_ `ao_pseudo_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L1>`_
Pseudo-potential Pseudo-potential
`ao_pseudo_integral_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L13>`_ `ao_pseudo_integral_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L15>`_
Local pseudo-potential Local pseudo-potential
`ao_pseudo_integral_non_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L119>`_ `ao_pseudo_integral_non_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L121>`_
Local pseudo-potential Local pseudo-potential
`mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_mo_ints.irp.f#L1>`_ `mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_mo_ints.irp.f#L1>`_

3
src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
View File

@ -5,6 +5,8 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
END_DOC END_DOC
if (do_pseudo) then if (do_pseudo) then
ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local
! ao_pseudo_integral = ao_pseudo_integral_local
! ao_pseudo_integral = ao_pseudo_integral_non_local
else else
ao_pseudo_integral = 0.d0 ao_pseudo_integral = 0.d0
endif endif
@ -221,3 +223,4 @@ END_PROVIDER
END_PROVIDER END_PROVIDER

1
src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f
View File

@ -30,3 +30,4 @@ BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_n
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
END_PROVIDER END_PROVIDER

147
src/Integrals_Monoelec/pseudopot.f90
View File

@ -201,7 +201,7 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
! |_ (_) (_ (_| | (/_ ! |_ (_) (_ (_| | (/_
! !
double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm double precision :: fourpi,f,prod,prodp,binom_func,accu,bigR,bigI,ylm
double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big
double precision :: areal,freal,breal,t1,t2,int_prod_bessel double precision :: areal,freal,breal,t1,t2,int_prod_bessel
double precision :: arg double precision :: arg
@ -327,7 +327,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
do k1=0,n_a(1) do k1=0,n_a(1)
do k2=0,n_a(2) do k2=0,n_a(2)
do k3=0,n_a(3) do k3=0,n_a(3)
array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3) &
*(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)
enddo enddo
enddo enddo
@ -336,7 +336,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
do k1p=0,n_b(1) do k1p=0,n_b(1)
do k2p=0,n_b(2) do k2p=0,n_b(2)
do k3p=0,n_b(3) do k3p=0,n_b(3)
array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p) &
*(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
enddo enddo
enddo enddo
@ -448,7 +448,7 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
do k2p=0,n_b(2) do k2p=0,n_b(2)
do k3p=0,n_b(3) do k3p=0,n_b(3)
array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p) &
*(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
enddo enddo
enddo enddo
@ -534,7 +534,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
do k2=0,n_a(2) do k2=0,n_a(2)
do k3=0,n_a(3) do k3=0,n_a(3)
array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3) &
*(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)
enddo enddo
@ -705,7 +705,7 @@ integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p
double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0 double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0
double precision,allocatable :: array_R_loc(:,:,:) double precision,allocatable :: array_R_loc(:,:,:)
double precision,allocatable :: array_coefs(:,:,:,:,:,:) double precision,allocatable :: array_coefs(:,:,:,:,:,:)
double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
fourpi=4.d0*dacos(-1.d0) fourpi=4.d0*dacos(-1.d0)
f=fourpi**1.5d0 f=fourpi**1.5d0
@ -762,9 +762,9 @@ allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:n
do k1p=0,n_b(1) do k1p=0,n_b(1)
do k2p=0,n_b(2) do k2p=0,n_b(2)
do k3p=0,n_b(3) do k3p=0,n_b(3)
array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3) &
*(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) & *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) &
*binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p) &
*(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
enddo enddo
enddo enddo
@ -823,7 +823,7 @@ double precision function bigI(lambda,mu,l,m,k1,k2,k3)
implicit none implicit none
integer lambda,mu,l,m,k1,k2,k3 integer lambda,mu,l,m,k1,k2,k3
integer k,i,kp,ip integer k,i,kp,ip
double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom,fact,coef_pm double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm
pi=dacos(-1.d0) pi=dacos(-1.d0)
if(mu.gt.0.and.m.gt.0)then if(mu.gt.0.and.m.gt.0)then
@ -834,8 +834,8 @@ do k=0,mu/2
do i=0,lambda-mu do i=0,lambda-mu
do kp=0,m/2 do kp=0,m/2
do ip=0,l-m do ip=0,l-m
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3) sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3)
enddo enddo
enddo enddo
@ -868,7 +868,7 @@ do i=0,lambda
do kp=0,m/2 do kp=0,m/2
do ip=0,l-m do ip=0,l-m
cylm=factor1*coef_pm(lambda,i) cylm=factor1*coef_pm(lambda,i)
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3) sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3)
enddo enddo
enddo enddo
@ -884,7 +884,7 @@ factor2=dsqrt((2*l+1)/(4.d0*pi))
do k=0,mu/2 do k=0,mu/2
do i=0,lambda-mu do i=0,lambda-mu
do ip=0,l do ip=0,l
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
cylmp=factor2*coef_pm(l,ip) cylmp=factor2*coef_pm(l,ip)
sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3) sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3)
enddo enddo
@ -904,8 +904,8 @@ do k=0,(mu-1)/2
do i=0,lambda-mu do i=0,lambda-mu
do kp=0,(m-1)/2 do kp=0,(m-1)/2
do ip=0,l-m do ip=0,l-m
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3) sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3)
enddo enddo
enddo enddo
@ -926,7 +926,7 @@ do i=0,lambda
do kp=0,(m-1)/2 do kp=0,(m-1)/2
do ip=0,l-m do ip=0,l-m
cylm=factor1*coef_pm(lambda,i) cylm=factor1*coef_pm(lambda,i)
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3) sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3)
enddo enddo
enddo enddo
@ -944,7 +944,7 @@ factor2=dsqrt((2*l+1)/(4.d0*pi))
do k=0,(mu-1)/2 do k=0,(mu-1)/2
do i=0,lambda-mu do i=0,lambda-mu
do ip=0,l do ip=0,l
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
cylmp=factor2*coef_pm(l,ip) cylmp=factor2*coef_pm(l,ip)
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3) sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3)
enddo enddo
@ -964,8 +964,8 @@ do k=0,mu/2
do i=0,lambda-mu do i=0,lambda-mu
do kp=0,(m-1)/2 do kp=0,(m-1)/2
do ip=0,l-m do ip=0,l-m
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3) sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3)
enddo enddo
enddo enddo
@ -985,8 +985,8 @@ do k=0,(mu-1)/2
do i=0,lambda-mu do i=0,lambda-mu
do kp=0,m/2 do kp=0,m/2
do ip=0,l-m do ip=0,l-m
cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3) sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3)
enddo enddo
enddo enddo
@ -1699,14 +1699,14 @@ end
double precision function coef_pm(n,k) double precision function coef_pm(n,k)
implicit none implicit none
integer n,k integer n,k
double precision arg,binom,binom_gen double precision arg,binom_func,binom_gen
if(n.eq.0.and.k.ne.0)stop 'coef_pm not defined' if(n.eq.0.and.k.ne.0)stop 'coef_pm not defined'
if(n.eq.0.and.k.eq.0)then if(n.eq.0.and.k.eq.0)then
coef_pm=1.d0 coef_pm=1.d0
return return
endif endif
arg=0.5d0*dfloat(n+k-1) arg=0.5d0*dfloat(n+k-1)
coef_pm=2.d0**n*binom(n,k)*binom_gen(arg,n) coef_pm=2.d0**n*binom_func(n,k)*binom_gen(arg,n)
end end
!! Ylm_bis uses the series expansion of Ylm in xchap^i ychap^j zchap^k !! Ylm_bis uses the series expansion of Ylm in xchap^i ychap^j zchap^k
@ -1735,7 +1735,7 @@ end
double precision function ylm_bis(l,m,theta,phi) double precision function ylm_bis(l,m,theta,phi)
implicit none implicit none
integer l,m,k,i integer l,m,k,i
double precision x,y,z,theta,phi,sum,factor,pi,binom,fact,coef_pm,cylm double precision x,y,z,theta,phi,sum,factor,pi,binom_func,fact,coef_pm,cylm
pi=dacos(-1.d0) pi=dacos(-1.d0)
x=dsin(theta)*dcos(phi) x=dsin(theta)*dcos(phi)
y=dsin(theta)*dsin(phi) y=dsin(theta)*dsin(phi)
@ -1745,7 +1745,7 @@ if(m.gt.0)then
sum=0.d0 sum=0.d0
do k=0,m/2 do k=0,m/2
do i=0,l-m do i=0,l-m
cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom(m,2*k)*fact(m+i)/fact(i)*coef_pm(l,i+m) cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k)*fact(m+i)/fact(i)*coef_pm(l,i+m)
sum=sum+cylm*x**(m-2*k)*y**(2*k)*z**i sum=sum+cylm*x**(m-2*k)*y**(2*k)*z**i
enddo enddo
enddo enddo
@ -1765,7 +1765,7 @@ m=-m
sum=0.d0 sum=0.d0
do k=0,(m-1)/2 do k=0,(m-1)/2
do i=0,l-m do i=0,l-m
cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom(m,2*k+1)*fact(m+i)/fact(i)*coef_pm(l,i+m) cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k+1)*fact(m+i)/fact(i)*coef_pm(l,i+m)
sum=sum+cylm*x**(m-(2*k+1))*y**(2*k+1)*z**i sum=sum+cylm*x**(m-(2*k+1))*y**(2*k+1)*z**i
enddo enddo
enddo enddo
@ -1800,13 +1800,6 @@ end
!! !!
!! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n ) !! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n )
double precision function binom(i,j)
implicit none
integer :: i,j
double precision :: fact
binom = fact(i)/(fact(j)*fact(i-j))
end
double precision function binom_gen(alpha,n) double precision function binom_gen(alpha,n)
implicit none implicit none
integer :: n,k integer :: n,k
@ -2087,4 +2080,90 @@ end
end end
! l,m : Y(l,m) parameters
! c(3) : pseudopotential center
! a(3) : Atomic Orbital center
! n_a(3) : Powers of x,y,z in the Atomic Orbital
! g_a : Atomic Orbital exponent
! r : Distance between the Atomic Orbital center and the considered point
double precision function ylm_orb(l,m,c,a,n_a,g_a,r)
implicit none
integer lmax_max,ntot_max
parameter (lmax_max=2)
parameter (ntot_max=14)
integer l,m
double precision a(3),g_a,c(3)
double precision prod,binom_func,accu,bigI,ylm,bessel_mod
double precision theta_AC0,phi_AC0,ac,factor,fourpi,arg,r,areal
integer ntotA,mu,k1,k2,k3,lambda
integer n_a(3)
double precision &
array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)
double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max), y
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
arg=g_a*(ac**2+r**2)
fourpi=4.d0*dacos(-1.d0)
factor=fourpi*dexp(-arg)
areal=2.d0*g_a*ac
ntotA=n_a(1)+n_a(2)+n_a(3)
if(ntotA.gt.ntot_max)stop 'increase ntot_max'
if(ac.eq.0.d0)then
ylm_orb=dsqrt(fourpi)*r**ntotA*dexp(-g_a*r**2)*bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
return
else
theta_AC0=dacos( (a(3)-c(3))/ac )
phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac)
do k1=0,n_a(1)
do k2=0,n_a(2)
do k3=0,n_a(3)
array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3) &
*(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) &
*r**(k1+k2+k3)
enddo
enddo
enddo
do lambda=0,l+ntotA
do mu=-lambda,lambda
do k1=0,n_a(1)
do k2=0,n_a(2)
do k3=0,n_a(3)
array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
enddo
enddo
enddo
enddo
enddo
accu=0.d0
do lambda=0,l+ntotA
do mu=-lambda,lambda
y = ylm(lambda,mu,theta_AC0,phi_AC0)
if (y == 0.d0) then
cycle
endif
do k1=0,n_a(1)
do k2=0,n_a(2)
do k3=0,n_a(3)
prod=y*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3)
if (prod == 0.d0) then
cycle
endif
if (areal*r < 100.d0) then ! overflow!
accu=accu+prod*bessel_mod(areal*r,lambda)
endif
enddo
enddo
enddo
enddo
enddo
ylm_orb=factor*accu
return
endif
end

49
src/MRCC/README.rst
View File

@ -20,6 +20,51 @@ Documentation
.. Do not edit this section. It was auto-generated from the .. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES_CHILDREN file by the `update_README.py` script. .. NEEDED_MODULES_CHILDREN file by the `update_README.py` script.
`davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/davidson.irp.f#L51>`_
Davidson diagonalization with specific diagonal elements of the H matrix
.br
H_jj : specific diagonal H matrix elements to diagonalize de Davidson
.br
dets_in : bitmasks corresponding to determinants
.br
u_in : guess coefficients on the various states. Overwritten
on exit
.br
dim_in : leftmost dimension of u_in
.br
sze : Number of determinants
.br
N_st : Number of eigenstates
.br
iunit : Unit for the I/O
.br
Initial guess vectors are not necessarily orthonormal
`davidson_diag_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/davidson.irp.f#L1>`_
Davidson diagonalization.
.br
dets_in : bitmasks corresponding to determinants
.br
u_in : guess coefficients on the various states. Overwritten
on exit
.br
dim_in : leftmost dimension of u_in
.br
sze : Number of determinants
.br
N_st : Number of eigenstates
.br
iunit : Unit number for the I/O
.br
Initial guess vectors are not necessarily orthonormal
`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/davidson.irp.f#L355>`_
Computes v_0 = H|u_0>
.br
n : number of determinants
.br
H_jj : array of <j|H|j>
`mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc.irp.f#L1>`_ `mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc.irp.f#L1>`_
Undocumented Undocumented
@ -53,7 +98,7 @@ Documentation
`ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L78>`_ `ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L78>`_
Eigenvectors/values of the CI matrix Eigenvectors/values of the CI matrix
`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L144>`_ `ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L145>`_
N_states lowest eigenvalues of the dressed CI matrix N_states lowest eigenvalues of the dressed CI matrix
`delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L43>`_ `delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L43>`_
@ -62,7 +107,7 @@ Documentation
`delta_ij_non_cas <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L34>`_ `delta_ij_non_cas <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L34>`_
Dressing matrix in SD basis Dressing matrix in SD basis
`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L159>`_ `diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L160>`_
Replace the coefficients of the CI states by the coefficients of the Replace the coefficients of the CI states by the coefficients of the
eigenstates of the CI matrix eigenstates of the CI matrix

414
src/MRCC/davidson.irp.f Normal file
View File

@ -0,0 +1,414 @@
subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,istate)
use bitmasks
implicit none
BEGIN_DOC
! Davidson diagonalization.
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! iunit : Unit number for the I/O
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint, iunit, istate
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: energies(N_st)
double precision, allocatable :: H_jj(:)
double precision :: diag_h_mat_elem
integer :: i
ASSERT (N_st > 0)
ASSERT (sze > 0)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
PROVIDE mo_bielec_integrals_in_map
allocate(H_jj(sze))
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(sze,H_jj,dets_in,Nint,istate,delta_ij) &
!$OMP PRIVATE(i)
!$OMP DO SCHEDULE(guided)
do i=1,sze
H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) + delta_ij(i,i,istate)
enddo
!$OMP END DO
!$OMP END PARALLEL
call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit,istate)
deallocate (H_jj)
end
subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit,istate)
use bitmasks
implicit none
BEGIN_DOC
! Davidson diagonalization with specific diagonal elements of the H matrix
!
! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! iunit : Unit for the I/O
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint, istate
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(in) :: H_jj(sze)
integer, intent(in) :: iunit
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: energies(N_st)
integer :: iter
integer :: i,j,k,l,m
logical :: converged
double precision :: overlap(N_st,N_st)
double precision :: u_dot_v, u_dot_u
integer, allocatable :: kl_pairs(:,:)
integer :: k_pairs, kl
integer :: iter2
double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:)
double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:)
double precision :: diag_h_mat_elem
double precision :: residual_norm(N_st)
character*(16384) :: write_buffer
double precision :: to_print(2,N_st)
double precision :: cpu, wall
PROVIDE det_connections
call write_time(iunit)
call wall_time(wall)
call cpu_time(cpu)
write(iunit,'(A)') ''
write(iunit,'(A)') 'Davidson Diagonalization'
write(iunit,'(A)') '------------------------'
write(iunit,'(A)') ''
call write_int(iunit,N_st,'Number of states')
call write_int(iunit,sze,'Number of determinants')
write(iunit,'(A)') ''
write_buffer = '===== '
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ ================'
enddo
write(iunit,'(A)') trim(write_buffer)
write_buffer = ' Iter'
do i=1,N_st
write_buffer = trim(write_buffer)//' Energy Residual'
enddo
write(iunit,'(A)') trim(write_buffer)
write_buffer = '===== '
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ ================'
enddo
write(iunit,'(A)') trim(write_buffer)
allocate( &
kl_pairs(2,N_st*(N_st+1)/2), &
W(sze,N_st,davidson_sze_max), &
U(sze,N_st,davidson_sze_max), &
R(sze,N_st), &
h(N_st,davidson_sze_max,N_st,davidson_sze_max), &
y(N_st,davidson_sze_max,N_st,davidson_sze_max), &
lambda(N_st*davidson_sze_max))
ASSERT (N_st > 0)
ASSERT (sze > 0)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
! Initialization
! ==============
k_pairs=0
do l=1,N_st
do k=1,l
k_pairs+=1
kl_pairs(1,k_pairs) = k
kl_pairs(2,k_pairs) = l
enddo
enddo
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, &
!$OMP Nint,dets_in,u_in) &
!$OMP PRIVATE(k,l,kl,i)
! Orthonormalize initial guess
! ============================
!$OMP DO
do kl=1,k_pairs
k = kl_pairs(1,kl)
l = kl_pairs(2,kl)
if (k/=l) then
overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze)
overlap(l,k) = overlap(k,l)
else
overlap(k,k) = u_dot_u(U_in(1,k),sze)
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze)
! Davidson iterations
! ===================
converged = .False.
do while (.not.converged)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st)
do k=1,N_st
!$OMP DO
do i=1,sze
U(i,k,1) = u_in(i,k)
enddo
!$OMP END DO
enddo
!$OMP END PARALLEL
do iter=1,davidson_sze_max-1
! Compute W_k = H |u_k>
! ----------------------
do k=1,N_st
call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate)
enddo
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
do l=1,N_st
do k=1,N_st
do iter2=1,iter-1
h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze)
h(k,iter,l,iter2) = h(k,iter2,l,iter)
enddo
enddo
do k=1,l
h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze)
h(l,iter,k,iter) = h(k,iter,l,iter)
enddo
enddo
!DEBUG H MATRIX
!do i=1,iter
! print '(10(x,F16.10))', h(1,i,1,1:i)
!enddo
!print *, ''
!END
! Diagonalize h
! -------------
call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter)
! Express eigenvectors of h in the determinant basis
! --------------------------------------------------
do k=1,N_st
do i=1,sze
U(i,k,iter+1) = 0.d0
W(i,k,iter+1) = 0.d0
do l=1,N_st
do iter2=1,iter
U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
enddo
enddo
enddo
enddo
! Compute residual vector
! -----------------------
do k=1,N_st
do i=1,sze
R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1)
enddo
residual_norm(k) = u_dot_u(R(1,k),sze)
to_print(1,k) = lambda(k) + nuclear_repulsion
to_print(2,k) = residual_norm(k)
enddo
write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st)
call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
if (converged) then
exit
endif
! Davidson step
! -------------
do k=1,N_st
do i=1,sze
U(i,k,iter+1) = -1.d0/max(H_jj(i) - lambda(k),1.d-2) * R(i,k)
enddo
enddo
! Gram-Schmidt
! ------------
double precision :: c
do k=1,N_st
do iter2=1,iter
do l=1,N_st
c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze)
do i=1,sze
U(i,k,iter+1) -= c * U(i,l,iter2)
enddo
enddo
enddo
do l=1,k-1
c = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze)
do i=1,sze
U(i,k,iter+1) -= c * U(i,l,iter+1)
enddo
enddo
call normalize( U(1,k,iter+1), sze )
enddo
!DEBUG : CHECK OVERLAP
!print *, '==='
!do k=1,iter+1
! do l=1,k
! c = u_dot_v(U(1,1,k),U(1,1,l),sze)
! print *, k,l, c
! enddo
!enddo
!print *, '==='
!pause
!END DEBUG
enddo
if (.not.converged) then
iter = davidson_sze_max-1
endif
! Re-contract to u_in
! -----------
do k=1,N_st
energies(k) = lambda(k)
do i=1,sze
u_in(i,k) = 0.d0
do iter2=1,iter
do l=1,N_st
u_in(i,k) += U(i,l,iter2)*y(l,iter2,k,1)
enddo
enddo
enddo
enddo
enddo
write_buffer = '===== '
do i=1,N_st
write_buffer = trim(write_buffer)//' ================ ================'
enddo
write(iunit,'(A)') trim(write_buffer)
write(iunit,'(A)') ''
call write_time(iunit)
deallocate ( &
kl_pairs, &
W, &
U, &
R, &
h, &
y, &
lambda &
)
abort_here = abort_all
end
subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
END_DOC
integer, intent(in) :: n,Nint,istate
double precision, intent(out) :: v_0(n)
double precision, intent(in) :: u_0(n)
double precision, intent(in) :: H_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
integer, allocatable :: idx(:)
double precision :: hij
double precision, allocatable :: vt(:)
integer :: i,j,k,l, jj
integer :: i0, j0
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy delta_ij
integer, parameter :: block_size = 157
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,j,k,idx,jj,vt) &
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij)
!$OMP DO SCHEDULE(static)
do i=1,n
v_0(i) = H_jj(i) * u_0(i)
enddo
!$OMP END DO
allocate(idx(0:n), vt(n))
Vt = 0.d0
!$OMP DO SCHEDULE(guided)
do i=1,n
idx(0) = i
call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
do jj=1,idx(0)
j = idx(jj)
if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
hij = hij + delta_ij(j,i,istate)
vt (i) = vt (i) + hij*u_0(j)
vt (j) = vt (j) + hij*u_0(i)
endif
enddo
enddo
!$OMP END DO
!$OMP CRITICAL
do i=1,n
v_0(i) = v_0(i) + vt(i)
enddo
!$OMP END CRITICAL
deallocate(idx,vt)
!$OMP END PARALLEL
end

7
src/MRCC/mrcc_utils.irp.f
View File

@ -92,9 +92,10 @@ END_PROVIDER
if (diag_algorithm == "Davidson") then if (diag_algorithm == "Davidson") then
stop 'use Lapack' integer :: istate
! call davidson_diag(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, & istate = 1
! size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants) call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, &
size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants,istate)
else if (diag_algorithm == "Lapack") then else if (diag_algorithm == "Lapack") then

53
src/Makefile.common Normal file
View File

@ -0,0 +1,53 @@
.PHONY: default silent
export
ifneq ($(IN_MAKE),1)
default:
@$(MAKE) -C $(QPACKAGE_ROOT)/src $$(basename $(PWD))
veryclean:
clean_modules.sh
clean:
IN_MAKE=1 make clean
else # Called by scripts/build_module.sh
default: all .gitignore
# Include the user's config
include $(QPACKAGE_ROOT)/src/Makefile.config
# Create the NEEDED_CHILDREN_MODULES variable, needed for IRPF90
NEEDED_CHILDREN_MODULES=$(shell module_handler.py print_genealogy)
# Check and update dependencies
include Makefile.depend
# Define the Makefile common variables
EZFIO_DIR=$(QPACKAGE_ROOT)/EZFIO
EZFIO=$(EZFIO_DIR)/lib/libezfio_irp.a
INCLUDE_DIRS=$(NEEDED_CHILDREN_MODULES) include
clean_links:
rm -f $(INCLUDE_DIRS) $$(basename $$PWD)
LIB+=$(EZFIO) $(MKL)
IRPF90+=$(patsubst %, -I %, $(INCLUDE_DIRS)) $(IRPF90_FLAGS)
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO) NEEDED_CHILDREN_MODULES $(wildcard *.py)
- $(IRPF90)
- update_README.py
include irpf90.make
endif
.gitignore:
$(QPACKAGE_ROOT)/scripts/create/create_gitignore.sh
# Frequent typos
clena: clean
veryclena: veryclean
vercylean: veryclean

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

@ -0,0 +1,34 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90_FLAGS+= --align=32
FC = ifort -g
#FC = cache_compile.py ifort -g # Accelerates compilation
FCFLAGS=
FCFLAGS+= -xHost
#FCFLAGS+= -xAVX
FCFLAGS+= -O2
FCFLAGS+= -ip
FCFLAGS+= -opt-prefetch
FCFLAGS+= -ftz
MKL=-mkl=parallel
NINJA=ninja
ifeq ($(PROFILE),1)
FC += -p -g
CXX += -pg
endif
ifeq ($(OPENMP),1)
FC += -openmp
IRPF90_FLAGS += --openmp
CXX += -fopenmp
endif
ifeq ($(DEBUG),1)
FC += -C -traceback -fpe0
FCFLAGS+= -axSSE2
IRPF90_FLAGS += -a
#FCFLAGS =-O0
endif

BIN
src/Molden/tree_dependency.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 32 KiB

After

Width:  |  Height:  |  Size: 30 KiB

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DensityFit DDCI_selected Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils QmcChem

2
src/Nuclei/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Ezfio_files Ezfio_files Utils

1
src/Nuclei/README.rst
View File

@ -15,6 +15,7 @@ Needed Modules
.. image:: tree_dependancy.png .. image:: tree_dependancy.png
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_ * `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation Documentation
============= =============

31
src/Pseudo/EZFIO.cfg
View File

@ -55,3 +55,34 @@ doc: Using pseudo potential integral of not
interface: input interface: input
default: False default: False
[pseudo_grid_size]
type: integer
doc: Nb of points of the QMC grid
interface: input
default: 1000
[pseudo_grid_rmax]
type: double precision
doc: R_maxof the QMC grid
interface: input
default: 10.0
[ao_pseudo_grid]
type: double precision
doc: QMC grid
interface: output
size: (ao_basis.ao_num,-pseudo.pseudo_lmax:pseudo.pseudo_lmax,0:pseudo.pseudo_lmax,nuclei.nucl_num,pseudo.pseudo_grid_size)
[mo_pseudo_grid]
type: double precision
doc: QMC grid
interface: output
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: output
size: (aux_basis.aux_basis_num_sqrt,aux_basis.aux_basis_num_sqrt)

4
src/Pseudo/README.rst
View File

@ -18,5 +18,9 @@ Needed Modules
.. image:: tree_dependancy.png .. image:: tree_dependancy.png
<<<<<<< HEAD
=======
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
>>>>>>> LCPQ-master
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_ * `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_

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

6
src/QmcChem/Makefile Normal file
View File

@ -0,0 +1,6 @@
# 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/QmcChem/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Determinants DensityFit

56
src/QmcChem/README.rst Normal file
View File

@ -0,0 +1,56 @@
==============
QmcChem Module
==============
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`ao_pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pot_ao_pseudo_ints.irp.f#L225>`_
Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
.br
<img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
`aux_pseudo_integral <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pot_ao_pseudo_ints.irp.f#L1>`_
Pseudo-potential
`aux_pseudo_integral_local <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pot_ao_pseudo_ints.irp.f#L15>`_
Local pseudo-potential
`aux_pseudo_integral_non_local <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pot_ao_pseudo_ints.irp.f#L121>`_
Local pseudo-potential
`mo_pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pot_ao_pseudo_ints.irp.f#L276>`_
Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
.br
<img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
`test_pseudo_grid_ao <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pot_ao_pseudo_ints.irp.f#L321>`_
Undocumented
`pseudo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pseudo.irp.f#L12>`_
Pseudo-potential expressed in the basis of ao products
`write_pseudopotential <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/pseudo.irp.f#L1>`_
Write the pseudo_potential into the EZFIO file
`save_for_qmc <http://github.com/LCPQ/quantum_package/tree/master/src/QmcChem/save_for_qmcchem.irp.f#L1>`_
Undocumented
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
.. image:: tree_dependancy.png
* `Determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants>`_
* `DensityFit <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit>`_

341
src/QmcChem/pot_ao_pseudo_ints.irp.f Normal file
View File

@ -0,0 +1,341 @@
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
BEGIN_DOC
! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
!
! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
END_DOC
! l,m : Y(l,m) parameters
! c(3) : pseudopotential center
! a(3) : Atomic Orbital center
! n_a(3) : Powers of x,y,z in the Atomic Orbital
! g_a : Atomic Orbital exponent
! r : Distance between the Atomic Orbital center and the considered point
double precision, external :: ylm_orb
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 :: y
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
do j=2,pseudo_grid_size
r(j) = r(j-1) + dr
enddo
ao_pseudo_grid = 0.d0
do j=1,pseudo_grid_size
do k=1,nucl_num
c(1:3) = nucl_coord(k,1:3)
do l=0,pseudo_lmax
do i=1,ao_num
a(1:3) = nucl_coord(ao_nucl(i),1:3)
n_a(1:3) = ao_power(i,1:3)
do p=1,ao_prim_num(i)
g_a = ao_expo_ordered_transp(p,i)
do m=-l,l
y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
ao_pseudo_grid(i,m,l,k,j) = ao_pseudo_grid(i,m,l,k,j) + &
ao_coef_normalized_ordered_transp(p,i)*y
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
!
! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
END_DOC
! l,m : Y(l,m) parameters
! c(3) : pseudopotential center
! a(3) : Atomic Orbital center
! n_a(3) : Powers of x,y,z in the Atomic Orbital
! g_a : Atomic Orbital exponent
! r : Distance between the Atomic Orbital center and the considered point
double precision, external :: ylm_orb
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 :: y
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
do j=2,pseudo_grid_size
r(j) = r(j-1) + dr
enddo
mo_pseudo_grid = 0.d0
do n=1,pseudo_grid_size
do k=1,nucl_num
do l=0,pseudo_lmax
do m=-l,l
do j=1,mo_tot_num
do i=1,ao_num
mo_pseudo_grid(j,m,l,k,n) = mo_pseudo_grid(j,m,l,k,n) + &
ao_pseudo_grid(i,m,l,k,n) * mo_coef(i,j)
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
double precision function test_pseudo_grid_ao(i,j)
implicit none
integer, intent(in) :: i,j
integer :: k,l,m,n
double precision :: r, dr,u
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
test_pseudo_grid_ao = 0.d0
r = 0.d0
do k=1,pseudo_grid_size
do n=1,nucl_num
do l = 0,pseudo_lmax
u = pseudo_v_kl(n,l,1) * exp(-pseudo_dz_kl(n,l,1)*r*r)* r*r*dr
do m=-l,l
test_pseudo_grid_ao += ao_pseudo_grid(i,m,l,n,k) * ao_pseudo_grid(j,m,l,n,k) * u
enddo
enddo
enddo
r = r+dr
enddo
end

171
src/QmcChem/pseudo.irp.f Normal file
View File

@ -0,0 +1,171 @@
subroutine write_pseudopotential
implicit none
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

8
src/QmcChem/save_for_qmcchem.irp.f Normal file
View File

@ -0,0 +1,8 @@
program save_for_qmc
read_wf = .True.
TOUCH read_wf
call write_spindeterminants
if (do_pseudo) then
call write_pseudopotential
endif
end

6
src/Utils/integration.irp.f
View File

@ -6,7 +6,7 @@ subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alph
! fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2) ! fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2)
END_DOC END_DOC
implicit none implicit none
include 'constants.F' include 'include/constants.F'
integer, intent(in) :: dim integer, intent(in) :: dim
integer, intent(in) :: a,b ! powers : (x-xa)**a_x = (x-A(1))**a(1) 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) :: alpha, beta ! exponents
@ -53,7 +53,7 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
! * [ 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 ) ! * [ 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 )
END_DOC END_DOC
implicit none implicit none
include 'constants.F' include 'include/constants.F'
integer, intent(in) :: dim integer, intent(in) :: dim
integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) 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) :: alpha, beta ! exponents
@ -131,7 +131,7 @@ subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder
! * [ 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 ) ! * [ 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 )
END_DOC END_DOC
implicit none implicit none
include 'constants.F' include 'include/constants.F'
integer, intent(in) :: dim integer, intent(in) :: dim
integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
double precision, intent(in) :: alpha, beta, gama ! exponents double precision, intent(in) :: alpha, beta, gama ! exponents

2
src/Utils/need.irp.f
View File

@ -46,7 +46,7 @@
double precision function rinteg(n,u) double precision function rinteg(n,u)
implicit double precision(a-h,o-z) implicit double precision(a-h,o-z)
include 'constants.F' include '../include/constants.F'
! pi=dacos(-1.d0) ! pi=dacos(-1.d0)
ichange=1 ichange=1
factor=1.d0 factor=1.d0

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

@ -34,7 +34,7 @@ end
subroutine overlap_A_B_C(dim,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap) subroutine overlap_A_B_C(dim,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap)
implicit none implicit none
include 'constants.F' include 'include/constants.F'
integer, intent(in) :: dim integer, intent(in) :: dim
integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
double precision, intent(in) :: alpha, beta, gama ! exponents double precision, intent(in) :: alpha, beta, gama ! exponents
@ -131,13 +131,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
integer :: iorder_p(3) 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) 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.1d-10)then ! if(fact_p.lt.1d-20)then
overlap_x = 0.d0 ! overlap_x = 0.d0
overlap_y = 0.d0 ! overlap_y = 0.d0
overlap_z = 0.d0 ! overlap_z = 0.d0
overlap = 0.d0 ! overlap = 0.d0
return ! return
endif ! endif
integer :: nmax integer :: nmax
double precision :: F_integral double precision :: F_integral
nmax = maxval(iorder_p) nmax = maxval(iorder_p)

4
tests/unit_test/unit_test.py
View File

@ -20,7 +20,7 @@ Energy = namedtuple('Energy', ['without_pseudo', 'with_pseudo'])
# O p t # # O p t #
# ~#~#~ # # ~#~#~ #
precision = 1.e-7 precision = 5.e-7
# A test get a geo file and a basis file. # A test get a geo file and a basis file.
# A global dict containt the result for this test # A global dict containt the result for this test
@ -169,7 +169,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None) ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174) ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
ref_energy["vdz"]["HBO"] = Energy(None, -19.11982540413317) ref_energy["vdz"]["HBO"] = Energy(None, -19.119821904)
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
# G l o b a l _ v a r i a b l e # # G l o b a l _ v a r i a b l e #