diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 6ec6e0d0..538c5f63 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -14,7 +14,7 @@ let spec = +> flag "m" (optional_with_default 1 int) ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1." +> flag "p" no_arg - ~doc:"Using pseudo." + ~doc:"Using pseudopotentials" +> anon ("xyz_file" %: string) ;; diff --git a/scripts/cache_compile.py b/scripts/cache_compile.py new file mode 100755 index 00000000..9898eae5 --- /dev/null +++ b/scripts/cache_compile.py @@ -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\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() diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py index fe7a619e..3bc73be6 100755 --- a/scripts/ezfio_interface/ei_handler.py +++ b/scripts/ezfio_interface/ei_handler.py @@ -365,7 +365,10 @@ def create_ezfio_stuff(dict_ezfio_cfg, config_or_default="config"): except ValueError: a_size_raw.append(dim) 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) diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index a2f07e4e..5db71852 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -51,5 +51,7 @@ then ${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" else ${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 diff --git a/scripts/install/install_docopt.sh b/scripts/install/install_docopt.sh index 65a01f0f..862e61cd 100755 --- a/scripts/install/install_docopt.sh +++ b/scripts/install/install_docopt.sh @@ -18,4 +18,4 @@ cd ${QPACKAGE_ROOT} rm -f -- scripts/${DOCOPT}{,c} ${QPACKAGE_ROOT}/scripts/install/fetch_from_web.py ${DOCOPT_URL} ${DOCOPT} -mv ${DOCOPT} scripts/utility/${DOCOPT} \ No newline at end of file +mv ${DOCOPT} scripts/utility/${DOCOPT} diff --git a/scripts/install/install_m4.sh b/scripts/install/install_m4.sh index 2ccf00a1..6495b82e 100755 --- a/scripts/install/install_m4.sh +++ b/scripts/install/install_m4.sh @@ -16,6 +16,7 @@ fi cd ${QPACKAGE_ROOT} +rm -f l${QPACKAGE_ROOT}/bin/m4 if [[ -z ${M4} ]] then rm -f -- bin/m4 diff --git a/scripts/install/install_ninja.sh b/scripts/install/install_ninja.sh index ec73708c..e54c5366 100755 --- a/scripts/install/install_ninja.sh +++ b/scripts/install/install_ninja.sh @@ -1,7 +1,7 @@ #!/bin/bash # -# Installs m4 for ocaml -# Thu Oct 23 22:02:08 CEST 2014 +# Installs the ninja build system +# Thu May 28 13:21:16 CEST 2015 BASE="ninja" URL="https://github.com/martine/ninja/archive/master.tar.gz" @@ -21,5 +21,3 @@ rm -rf ${BASE} mv ${BASE}-master ${BASE} cd ${BASE} ./configure.py --bootstrap - - diff --git a/scripts/module/create_Makefile_depend.sh b/scripts/module/create_Makefile_depend.sh new file mode 100755 index 00000000..2485bb97 --- /dev/null +++ b/scripts/module/create_Makefile_depend.sh @@ -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 + + diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 98641a18..d5dfb584 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -186,7 +186,7 @@ def add_zero(array, size, type): def make_it_square(matrix, dim, type=float): """ - matix the matrix to squate + matix the matrix to square dim array [lmax, kmax] type the null value you want [[[28.59107316], [19.37583724]], [[50.25646328]]] diff --git a/setup_environment.sh b/setup_environment.sh index 0d7001c4..0eb5e20d 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -29,7 +29,7 @@ EOF source quantum_package.rc - +mkdir -p install_logs echo "${BLUE}===== Installing IRPF90 ===== ${BLACK}" ${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 ]] diff --git a/src/Bitmask/bitmasks_module.f90 b/src/Bitmask/bitmasks_module.f90 index c542fac6..c5a9093f 100644 --- a/src/Bitmask/bitmasks_module.f90 +++ b/src/Bitmask/bitmasks_module.f90 @@ -8,4 +8,4 @@ module bitmasks integer, parameter :: d_part2 = 4 integer, parameter :: s_hole = 5 integer, parameter :: s_part = 6 -end module +end module bitmasks diff --git a/src/DensityFit/ASSUMPTIONS.rst b/src/DensityFit/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/DensityFit/Makefile b/src/DensityFit/Makefile new file mode 100644 index 00000000..06dc50ff --- /dev/null +++ b/src/DensityFit/Makefile @@ -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 diff --git a/src/DensityFit/NEEDED_CHILDREN_MODULES b/src/DensityFit/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..d6315de9 --- /dev/null +++ b/src/DensityFit/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +AOs Pseudo diff --git a/src/DensityFit/README.rst b/src/DensityFit/README.rst new file mode 100644 index 00000000..9d3c4d9b --- /dev/null +++ b/src/DensityFit/README.rst @@ -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 `_ + Exponents and coefficients of the auxiliary basis + +`aux_basis_coef_transp `_ + Exponents of the auxiliary basis + +`aux_basis_expo `_ + Exponents and coefficients of the auxiliary basis + +`aux_basis_expo_transp `_ + Exponents of the auxiliary basis + +`aux_basis_idx `_ + aux_basis_idx(k) -> i,j + +`aux_basis_nucl `_ + Exponents of the auxiliary basis + +`aux_basis_num `_ + Number of auxiliary basis functions + +`aux_basis_num_8 `_ + Number of auxiliary basis functions + +`aux_basis_num_sqrt `_ + Number of auxiliary basis functions + +`aux_basis_overlap_matrix `_ + Auxiliary basis set + +`aux_basis_power `_ + Exponents of the auxiliary basis + +`aux_basis_prim_num `_ + Exponents of the auxiliary basis + +`aux_basis_prim_num_max `_ + = ao_prim_num_max + +`save_aux_basis `_ + Undocumented + +`aux_basis_four_overlap `_ + \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 `_ +* `Pseudo `_ + diff --git a/src/DensityFit/aux_basis.ezfio_config b/src/DensityFit/aux_basis.ezfio_config new file mode 100644 index 00000000..cd56d1c7 --- /dev/null +++ b/src/DensityFit/aux_basis.ezfio_config @@ -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) + + diff --git a/src/DensityFit/aux_basis.irp.f b/src/DensityFit/aux_basis.irp.f new file mode 100644 index 00000000..009c59cf --- /dev/null +++ b/src/DensityFit/aux_basis.irp.f @@ -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 diff --git a/src/DensityFit/overlap.irp.f b/src/DensityFit/overlap.irp.f new file mode 100644 index 00000000..836df432 --- /dev/null +++ b/src/DensityFit/overlap.irp.f @@ -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 + + + diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index cb141fe3..9bf75ee3 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -532,12 +532,6 @@ Documentation `save_casino `_ Undocumented -`save_dets_qmcchem `_ - Undocumented - -`save_for_qmc `_ - Undocumented - `save_natorb `_ Undocumented diff --git a/src/Determinants/save_for_qmcchem.irp.f b/src/Determinants/save_for_qmcchem.irp.f deleted file mode 100644 index 93a9778d..00000000 --- a/src/Determinants/save_for_qmcchem.irp.f +++ /dev/null @@ -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 diff --git a/src/FCIdump/README.rst b/src/FCIdump/README.rst index 68d3a7d8..0efb57a5 100644 --- a/src/FCIdump/README.rst +++ b/src/FCIdump/README.rst @@ -10,9 +10,6 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES_CHILDREN file by the `update_README.py` script. -`fcidump `_ - Undocumented - Needed Modules diff --git a/src/Integrals_Monoelec/README.rst b/src/Integrals_Monoelec/README.rst index 4bc612a2..720929ef 100644 --- a/src/Integrals_Monoelec/README.rst +++ b/src/Integrals_Monoelec/README.rst @@ -97,10 +97,10 @@ Documentation `ao_pseudo_integral `_ Pseudo-potential -`ao_pseudo_integral_local `_ +`ao_pseudo_integral_local `_ Local pseudo-potential -`ao_pseudo_integral_non_local `_ +`ao_pseudo_integral_non_local `_ Local pseudo-potential `mo_nucl_elec_integral `_ diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f index 093b9fc8..95023177 100644 --- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -5,6 +5,8 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)] END_DOC 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 = ao_pseudo_integral_non_local else ao_pseudo_integral = 0.d0 endif @@ -221,3 +223,4 @@ END_PROVIDER END_PROVIDER + diff --git a/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f index b6851f21..6c412e4b 100644 --- a/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f @@ -30,3 +30,4 @@ BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_n !$OMP END PARALLEL DO END_PROVIDER + diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90 index 0cea5f66..b605a45e 100644 --- a/src/Integrals_Monoelec/pseudopot.f90 +++ b/src/Integrals_Monoelec/pseudopot.f90 @@ -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 :: areal,freal,breal,t1,t2,int_prod_bessel 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 k2=0,n_a(2) 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) 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 k2p=0,n_b(2) 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) 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 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) 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 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) 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,allocatable :: array_R_loc(:,:,:) 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) 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 k2p=0,n_b(2) 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) & - *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) enddo enddo @@ -823,7 +823,7 @@ double precision function bigI(lambda,mu,l,m,k1,k2,k3) implicit none integer lambda,mu,l,m,k1,k2,k3 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) if(mu.gt.0.and.m.gt.0)then @@ -834,8 +834,8 @@ do k=0,mu/2 do i=0,lambda-mu do kp=0,m/2 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) - cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + 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_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) enddo enddo @@ -868,7 +868,7 @@ do i=0,lambda do kp=0,m/2 do ip=0,l-m 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) enddo enddo @@ -884,7 +884,7 @@ factor2=dsqrt((2*l+1)/(4.d0*pi)) do k=0,mu/2 do i=0,lambda-mu 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) sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3) enddo @@ -904,8 +904,8 @@ do k=0,(mu-1)/2 do i=0,lambda-mu do kp=0,(m-1)/2 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) - cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + 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_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) enddo enddo @@ -926,7 +926,7 @@ do i=0,lambda do kp=0,(m-1)/2 do ip=0,l-m 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) enddo enddo @@ -944,7 +944,7 @@ factor2=dsqrt((2*l+1)/(4.d0*pi)) do k=0,(mu-1)/2 do i=0,lambda-mu 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) sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3) enddo @@ -964,8 +964,8 @@ do k=0,mu/2 do i=0,lambda-mu do kp=0,(m-1)/2 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) - cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + 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_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) enddo enddo @@ -985,8 +985,8 @@ do k=0,(mu-1)/2 do i=0,lambda-mu do kp=0,m/2 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) - cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + 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_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) enddo enddo @@ -1699,14 +1699,14 @@ end double precision function coef_pm(n,k) implicit none 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.eq.0)then coef_pm=1.d0 return endif 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 !! 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) implicit none 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) x=dsin(theta)*dcos(phi) y=dsin(theta)*dsin(phi) @@ -1745,7 +1745,7 @@ if(m.gt.0)then sum=0.d0 do k=0,m/2 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 enddo enddo @@ -1765,7 +1765,7 @@ m=-m sum=0.d0 do k=0,(m-1)/2 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 enddo enddo @@ -1800,13 +1800,6 @@ end !! !! 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) implicit none integer :: n,k @@ -2087,4 +2080,90 @@ 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 diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst index 25560b3a..2b649954 100644 --- a/src/MRCC/README.rst +++ b/src/MRCC/README.rst @@ -20,6 +20,51 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES_CHILDREN file by the `update_README.py` script. +`davidson_diag_hjj_mrcc `_ + 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 `_ + 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 `_ + Computes v_0 = H|u_0> + .br + n : number of determinants + .br + H_jj : array of + `mrcc `_ Undocumented @@ -53,7 +98,7 @@ Documentation `ci_electronic_energy_dressed `_ Eigenvectors/values of the CI matrix -`ci_energy_dressed `_ +`ci_energy_dressed `_ N_states lowest eigenvalues of the dressed CI matrix `delta_ij `_ @@ -62,7 +107,7 @@ Documentation `delta_ij_non_cas `_ Dressing matrix in SD basis -`diagonalize_ci_dressed `_ +`diagonalize_ci_dressed `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix diff --git a/src/MRCC/davidson.irp.f b/src/MRCC/davidson.irp.f new file mode 100644 index 00000000..52f03ff3 --- /dev/null +++ b/src/MRCC/davidson.irp.f @@ -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 = = + ! ------------------------------------------- + + 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 + 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 + + diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 611a70be..872588cf 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -92,9 +92,10 @@ END_PROVIDER if (diag_algorithm == "Davidson") then - stop 'use Lapack' -! call davidson_diag(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, & -! size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants) + integer :: istate + istate = 1 + 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 diff --git a/src/Makefile.common b/src/Makefile.common new file mode 100644 index 00000000..d6c9dd90 --- /dev/null +++ b/src/Makefile.common @@ -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 diff --git a/src/Makefile.config.example b/src/Makefile.config.example new file mode 100644 index 00000000..14272070 --- /dev/null +++ b/src/Makefile.config.example @@ -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 diff --git a/src/Molden/tree_dependency.png b/src/Molden/tree_dependency.png index fd419b86..23007e0f 100644 Binary files a/src/Molden/tree_dependency.png and b/src/Molden/tree_dependency.png differ diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index a525c891..77be07e1 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -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 \ No newline at end of file +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 \ No newline at end of file diff --git a/src/Nuclei/NEEDED_CHILDREN_MODULES b/src/Nuclei/NEEDED_CHILDREN_MODULES index d9eaf57c..dcdb5f86 100644 --- a/src/Nuclei/NEEDED_CHILDREN_MODULES +++ b/src/Nuclei/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Ezfio_files +Ezfio_files Utils diff --git a/src/Nuclei/README.rst b/src/Nuclei/README.rst index 577a8b34..0b2ae381 100644 --- a/src/Nuclei/README.rst +++ b/src/Nuclei/README.rst @@ -15,6 +15,7 @@ Needed Modules .. image:: tree_dependancy.png * `Ezfio_files `_ +* `Utils `_ Documentation ============= diff --git a/src/Pseudo/EZFIO.cfg b/src/Pseudo/EZFIO.cfg index 8a6afa2c..9804c807 100644 --- a/src/Pseudo/EZFIO.cfg +++ b/src/Pseudo/EZFIO.cfg @@ -55,3 +55,34 @@ doc: Using pseudo potential integral of not interface: input 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) + + diff --git a/src/Pseudo/README.rst b/src/Pseudo/README.rst index f27d7c20..45a589d4 100644 --- a/src/Pseudo/README.rst +++ b/src/Pseudo/README.rst @@ -18,5 +18,9 @@ Needed Modules .. image:: tree_dependancy.png +<<<<<<< HEAD +======= +* `Ezfio_files `_ +>>>>>>> LCPQ-master * `Nuclei `_ diff --git a/src/QmcChem/ASSUMPTIONS.rst b/src/QmcChem/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/QmcChem/Makefile b/src/QmcChem/Makefile new file mode 100644 index 00000000..06dc50ff --- /dev/null +++ b/src/QmcChem/Makefile @@ -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 diff --git a/src/QmcChem/NEEDED_CHILDREN_MODULES b/src/QmcChem/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..f7ed9913 --- /dev/null +++ b/src/QmcChem/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants DensityFit diff --git a/src/QmcChem/README.rst b/src/QmcChem/README.rst new file mode 100644 index 00000000..b1f1b738 --- /dev/null +++ b/src/QmcChem/README.rst @@ -0,0 +1,56 @@ +============== +QmcChem Module +============== + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`ao_pseudo_grid `_ + 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 + + +`aux_pseudo_integral `_ + Pseudo-potential + +`aux_pseudo_integral_local `_ + Local pseudo-potential + +`aux_pseudo_integral_non_local `_ + Local pseudo-potential + +`mo_pseudo_grid `_ + 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 + + +`test_pseudo_grid_ao `_ + Undocumented + +`pseudo_matrix `_ + Pseudo-potential expressed in the basis of ao products + +`write_pseudopotential `_ + Write the pseudo_potential into the EZFIO file + +`save_for_qmc `_ + Undocumented + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +.. image:: tree_dependancy.png + +* `Determinants `_ +* `DensityFit `_ + diff --git a/src/QmcChem/pot_ao_pseudo_ints.irp.f b/src/QmcChem/pot_ao_pseudo_ints.irp.f new file mode 100644 index 00000000..0aa2a98c --- /dev/null +++ b/src/QmcChem/pot_ao_pseudo_ints.irp.f @@ -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 +! +! + 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 +! +! + 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 diff --git a/src/QmcChem/pseudo.irp.f b/src/QmcChem/pseudo.irp.f new file mode 100644 index 00000000..abf4a5b8 --- /dev/null +++ b/src/QmcChem/pseudo.irp.f @@ -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 + diff --git a/src/QmcChem/save_for_qmcchem.irp.f b/src/QmcChem/save_for_qmcchem.irp.f new file mode 100644 index 00000000..4b028a7c --- /dev/null +++ b/src/QmcChem/save_for_qmcchem.irp.f @@ -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 diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f index a98efb3d..704ad631 100644 --- a/src/Utils/integration.irp.f +++ b/src/Utils/integration.irp.f @@ -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) END_DOC implicit none - include 'constants.F' + include 'include/constants.F' integer, intent(in) :: dim integer, intent(in) :: a,b ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta ! exponents @@ -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 ) END_DOC implicit none - include 'constants.F' + include 'include/constants.F' integer, intent(in) :: dim integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta ! exponents @@ -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 ) END_DOC implicit none - include 'constants.F' + include 'include/constants.F' integer, intent(in) :: dim integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta, gama ! exponents diff --git a/src/Utils/need.irp.f b/src/Utils/need.irp.f index 22cb6a48..f34dbad6 100644 --- a/src/Utils/need.irp.f +++ b/src/Utils/need.irp.f @@ -46,7 +46,7 @@ double precision function rinteg(n,u) implicit double precision(a-h,o-z) - include 'constants.F' + include '../include/constants.F' ! pi=dacos(-1.d0) ichange=1 factor=1.d0 diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index 3926bbc5..2e21b1c4 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -34,7 +34,7 @@ end subroutine overlap_A_B_C(dim,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap) implicit none - include 'constants.F' + include 'include/constants.F' integer, intent(in) :: dim integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) double precision, intent(in) :: alpha, beta, gama ! exponents @@ -131,13 +131,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& integer :: iorder_p(3) call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) - if(fact_p.lt.1d-10)then - overlap_x = 0.d0 - overlap_y = 0.d0 - overlap_z = 0.d0 - overlap = 0.d0 - return - endif +! if(fact_p.lt.1d-20)then +! overlap_x = 0.d0 +! overlap_y = 0.d0 +! overlap_z = 0.d0 +! overlap = 0.d0 +! return +! endif integer :: nmax double precision :: F_integral nmax = maxval(iorder_p) diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py index 53086de9..119d7676 100755 --- a/tests/unit_test/unit_test.py +++ b/tests/unit_test/unit_test.py @@ -20,7 +20,7 @@ Energy = namedtuple('Energy', ['without_pseudo', 'with_pseudo']) # O p t # # ~#~#~ # -precision = 1.e-7 +precision = 5.e-7 # A test get a geo file and a basis file. # 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["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 #