diff --git a/.travis.yml b/.travis.yml index 9aadf7ef..f9ad2d5f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -7,6 +7,7 @@ python: before_script: - sudo apt-get update - sudo apt-get install gfortran liblapack-dev + - sudo apt-get install graphviz script: - ./setup_environment.sh --robot diff --git a/INSTALL.md b/INSTALL.md index 47319c9c..22d52e39 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -10,6 +10,10 @@ * Bash * Patch (for opam) +## Optional + +* graphviz + ## Standard installation diff --git a/data/ezfio_defaults/WILL_BE_DELETED.ezfio_default b/data/ezfio_defaults/WILL_BE_DELETED.ezfio_default deleted file mode 100644 index 0d9489d5..00000000 --- a/data/ezfio_defaults/WILL_BE_DELETED.ezfio_default +++ /dev/null @@ -1,61 +0,0 @@ -bielec_integrals - read_ao_integrals false - read_mo_integrals false - write_ao_integrals false - write_mo_integrals false - threshold_ao 1.e-15 - threshold_mo 1.e-15 - direct false - -cis_dressed - n_state_cis 10 - n_core_cis 0 - n_act_cis mo_basis_mo_tot_num - mp2_dressing false - standard_doubles true - en_2_2 false - -determinants - n_states 1 - n_states_diag determinants_n_states - n_det_max_jacobi 1000 - threshold_generators 0.99 - threshold_selectors 0.999 - read_wf false - s2_eig false - only_single_double_dm false - -full_ci - n_det_max_fci 10000 - n_det_max_fci_property 50000 - pt2_max 1.e-4 - do_pt2_end true - var_pt2_ratio 0.75 - -cas_sd - n_det_max_cas_sd 100000 - pt2_max 1.e-4 - do_pt2_end true - var_pt2_ratio 0.75 - -all_singles - n_det_max_fci 50000 - pt2_max 1.e-8 - do_pt2_end false - -hartree_fock - n_it_scf_max 200 - thresh_scf 1.e-10 - guess "Huckel" - -cisd_selected - n_det_max_cisd 10000 - pt2_max 1.e-4 - -cisd_sc2_selected - n_det_max_cisd_sc2 10000 - pt2_max 1.e-4 - do_pt2_end true - -properties - z_one_point 3.9 diff --git a/data/ezfio_defaults/cisd.ezfio_default b/data/ezfio_defaults/cisd.ezfio_default deleted file mode 100644 index 95991f8f..00000000 --- a/data/ezfio_defaults/cisd.ezfio_default +++ /dev/null @@ -1,4 +0,0 @@ -cisd_selected - n_det_max_cisd 10000 - pt2_max 1.e-4 - 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 4667478f..4cba032f 100755 --- a/scripts/ezfio_interface/ei_handler.py +++ b/scripts/ezfio_interface/ei_handler.py @@ -197,7 +197,7 @@ def get_dict_config_file(config_file_path, module_lower): * equal to MODULE_lower name by default. - interface : The provider is a imput or a output - default : The default value /!\ stored in a Type named type! - if interface == output + if interface == input - size : Is the string read in ezfio.cgf who containt the size information (like 1 or =sum(ao_num)) """ @@ -230,7 +230,8 @@ def get_dict_config_file(config_file_path, module_lower): # Create the dictionary who containt the value per default d_default = {"ezfio_name": pvd, - "ezfio_dir": module_lower} + "ezfio_dir": module_lower, + "size": "1"} # Check if type if avalaible type_ = config_file.get(section, "type") @@ -274,6 +275,8 @@ def get_dict_config_file(config_file_path, module_lower): def create_ezfio_provider(dict_ezfio_cfg): + import re + """ From dict d[provider_name] = {type, doc, @@ -286,12 +289,13 @@ def create_ezfio_provider(dict_ezfio_cfg): output = output_dict_info['ezfio_dir' return [code, ...] """ + from ezfio_generate_provider import EZFIO_Provider dict_code_provider = dict() ez_p = EZFIO_Provider() for provider_name, dict_info in dict_ezfio_cfg.iteritems(): - if "default" in dict_info: + if "input" in dict_info["interface"]: ez_p.set_type(dict_info['type'].fortran) ez_p.set_name(provider_name) ez_p.set_doc(dict_info['doc']) @@ -299,6 +303,9 @@ def create_ezfio_provider(dict_ezfio_cfg): ez_p.set_ezfio_name(dict_info['ezfio_name']) ez_p.set_output("output_%s" % dict_info['ezfio_dir']) + # (nuclei.nucl_num,pseudo.klocmax) => (nucl_num,klocmax) + ez_p.set_size(re.sub(r'\w+\.', "", dict_info['size'])) + dict_code_provider[provider_name] = str(ez_p) + "\n" return dict_code_provider @@ -342,15 +349,32 @@ def create_ezfio_stuff(dict_ezfio_cfg, config_or_default="config"): def size_format_to_ezfio(size_raw): """ If size_raw == "=" is a formula -> do nothing; return + Else convert the born of a multidimential array + (12,begin:end) into (12,begin+end+1) for example If the value are between parenthses -> do nothing; return - Else put it in parenthsesis """ size_raw = str(size_raw) - if any([size_raw.startswith('='), - size_raw.startswith("(") and size_raw.endswith(")")]): + if size_raw.startswith('='): size_convert = size_raw else: + size_raw = provider_info["size"].translate(None, "()") + size_raw = size_raw.replace('.', '_') + + a_size_raw = [] + for dim in size_raw.split(","): + try: + (begin, end) = map(str.strip, dim.split(":")) + except ValueError: + a_size_raw.append(dim) + else: + 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_convert = "({0})".format(size_raw) return size_convert @@ -719,7 +743,6 @@ def save_ocaml_qp_edit(str_ocaml_qp_edit): if __name__ == "__main__": arguments = docopt(__doc__) - # ___ # | ._ o _|_ # _|_ | | | |_ diff --git a/scripts/ezfio_interface/ezfio_generate_provider.py b/scripts/ezfio_interface/ezfio_generate_provider.py index 7f3c8441..6cd919dc 100755 --- a/scripts/ezfio_interface/ezfio_generate_provider.py +++ b/scripts/ezfio_interface/ezfio_generate_provider.py @@ -10,10 +10,11 @@ fetched from the EZFIO file. import sys + class EZFIO_Provider(object): data = """ -BEGIN_PROVIDER [ %(type)s, %(name)s ] +BEGIN_PROVIDER [ %(type)s, %(name)s %(size)s ] implicit none BEGIN_DOC ! %(doc)s @@ -48,6 +49,9 @@ END_PROVIDER msg = "Error : %s is not set in EZFIO.cfg" % (v) print >>sys.stderr, msg sys.exit(1) + if "size" not in self.__dict__: + self.__dict__["size"] = "" + return self.data % self.__dict__ def set_write(self): @@ -83,6 +87,12 @@ END_PROVIDER def set_output(self, t): self.output = t + def set_size(self, t): + + if t != "1": + self.size = ", " + t + else: + self.size = "" def test_module(): T = EZFIO_Provider() diff --git a/scripts/ezfio_interface/qp_convert_ezfio_v1_to_v2.sh b/scripts/ezfio_interface/qp_convert_ezfio_v1_to_v2.sh index 3ef9a5f3..9c2829bc 100755 --- a/scripts/ezfio_interface/qp_convert_ezfio_v1_to_v2.sh +++ b/scripts/ezfio_interface/qp_convert_ezfio_v1_to_v2.sh @@ -15,10 +15,10 @@ mv $1/hartree_Fock/thresh_SCF $1/hartree_fock/thresh_scf 2> /dev/null # Set disk_acess echo "Change {read,write}_ao_integrals > disk_access_ao_integrals" -biint=$1/bielec_integrals +biint=$1/Integrals_bielec if [[ -f $biint/read_ao_integrals ]]; then - if [[ `cat $1/bielec_integrals/read_ao_integrals` -eq "T" ]] + if [[ `cat $1/Integrals_bielec/read_ao_integrals` -eq "T" ]] then echo "Read" > $biint/disk_access_ao_integrals @@ -33,4 +33,6 @@ if [[ -f $biint/read_ao_integrals ]]; then rm $biint/read_ao_integrals $biint/write_ao_integrals $biint/write_ao_intergals 2> /dev/null fi +mv $1/MonoInts $1/Integrals_Monoelec + echo "Done" \ No newline at end of file diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py index 400def37..0857f598 100755 --- a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py +++ b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py @@ -273,7 +273,7 @@ def write_ezfio(res, filename): # \_| |___/\___|\__,_|\__,_|\___/ # - ezfio.set_pseudo_integrals_do_pseudo(False) + ezfio.set_pseudo_do_pseudo(False) def get_full_path(file_path): diff --git a/scripts/install/install_m4.sh b/scripts/install/install_m4.sh index 6f012d04..8bf514fd 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/module/module_handler.py b/scripts/module/module_handler.py index a35dbe2c..99021f48 100755 --- a/scripts/module/module_handler.py +++ b/scripts/module/module_handler.py @@ -121,7 +121,12 @@ def create_png_from_path(path): "path = /home/razoa/quantum_package/src/Molden/NEEDED_CHILDREN_MODULES" l_module = os.path.split(path)[0].split("/")[-1] - create_png([l_module]) + + import pydot + try: + create_png([l_module]) + except pydot.InvocationException: + pass def create_png(l_module): @@ -154,8 +159,10 @@ def create_png(l_module): draw_module_edge(module, d_ref[module]) # Save - path = '{0}.png'.format("_".join(l_module)) - print "png saved in {0}".format(path) + path = '{0}.png'.format("tree_dependancy") + # path = '{0}.png'.format("_".join(l_module)) + # print "png saved in {0}".format(path) + graph.write_png(path) if __name__ == '__main__': diff --git a/scripts/module/pydot.py b/scripts/module/pydot.py index c788c137..c4a84bb6 100644 --- a/scripts/module/pydot.py +++ b/scripts/module/pydot.py @@ -31,7 +31,8 @@ import sys try: import dot_parser except Exception as e: - print >> sys.stderr, "Couldn't import dot_parser, loading of dot files will not be possible." + pass + # print >> sys.stderr, "Couldn't import dot_parser, loading of dot files will not be possible." GRAPH_ATTRIBUTES = set(['Damping', 'K', 'URL', 'aspect', 'bb', 'bgcolor', diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 1f8594c5..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]]] @@ -311,11 +311,11 @@ if __name__ == "__main__": # ~#~#~#~#~ # klocmax = max([len(i) for i in v_k]) - ezfio.pseudo_integrals_klocmax = klocmax + ezfio.pseudo_pseudo_klocmax = klocmax - ezfio.pseudo_integrals_v_k = zip(*v_k) - ezfio.pseudo_integrals_n_k = zip(*n_k) - ezfio.pseudo_integrals_dz_k = zip(*dz_k) + ezfio.pseudo_pseudo_v_k = zip(*v_k) + ezfio.pseudo_pseudo_n_k = zip(*n_k) + ezfio.pseudo_pseudo_dz_k = zip(*dz_k) # ~#~#~#~#~#~#~#~#~ # # N o n _ L o c a l # @@ -324,15 +324,15 @@ if __name__ == "__main__": lmax = max([len(i) for i in v_kl]) kmax = max([len(sublist) for list_ in v_kl for sublist in list_]) - ezfio.pseudo_integrals_lmaxpo = lmax - ezfio.pseudo_integrals_kmax = kmax + ezfio.pseudo_pseudo_lmax = lmax - 1 + ezfio.pseudo_pseudo_kmax = kmax v_kl = make_it_square(v_kl, [lmax, kmax]) n_kl = make_it_square(n_kl, [lmax, kmax], int) dz_kl = make_it_square(dz_kl, [lmax, kmax]) - ezfio.pseudo_integrals_v_kl = zip(*v_kl) - ezfio.pseudo_integrals_n_kl = zip(*n_kl) - ezfio.pseudo_integrals_dz_kl = zip(*dz_kl) + ezfio.pseudo_pseudo_v_kl = zip(*v_kl) + ezfio.pseudo_pseudo_n_kl = zip(*n_kl) + ezfio.pseudo_pseudo_dz_kl = zip(*dz_kl) - ezfio.pseudo_integrals_do_pseudo = True + ezfio.pseudo_do_pseudo = True diff --git a/scripts/run_Makefile_common.sh b/scripts/run_Makefile_common.sh index 859cc8d0..8673611e 100755 --- a/scripts/run_Makefile_common.sh +++ b/scripts/run_Makefile_common.sh @@ -66,3 +66,6 @@ ${QPACKAGE_ROOT}/scripts/module/create_Makefile_depend.sh # Update EZFIO interface ${QPACKAGE_ROOT}/scripts/ezfio_interface/ei_handler.py + +# Create png +${QPACKAGE_ROOT}/scripts/module/module_handler.py create_png \ No newline at end of file diff --git a/scripts/update_README.py b/scripts/update_README.py index 177d57a4..0bec525e 100755 --- a/scripts/update_README.py +++ b/scripts/update_README.py @@ -2,7 +2,7 @@ """Updates the README.rst file as the include directive is disabled on GitHub.""" __date__ = "Thu Apr 3 23:06:18 CEST 2014" -__author__ = "Anthony Scemama " +__author__ = "Anthony Scemama & TApplencourt " README = "README.rst" @@ -21,12 +21,6 @@ header = """ """ -try: - # subprocess.check_output("git status".split()) - has_git = True -except OSError: - has_git = False - def fetch_splitted_data(): """Read the README.rst file and split it in strings: @@ -87,10 +81,12 @@ def update_needed(data): modules = file.read() file.close() + header_image = ".. image:: tree_dependancy.png\n\n" + if modules.strip() != "": modules = ['* `%s <%s%s>`_' % (x, URL, x) for x in modules.split()] modules = "\n".join(modules) - modules = Needed_key + header + modules + '\n\n' + modules = Needed_key + header + header_image + modules + '\n\n' has_modules = False for i in range(len(data)): diff --git a/setup_environment.sh b/setup_environment.sh index 878c23e1..d62bd198 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -28,7 +28,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/AOs/NEEDED_CHILDREN_MODULES b/src/AOs/NEEDED_CHILDREN_MODULES index 2c80725f..ac848c17 100644 --- a/src/AOs/NEEDED_CHILDREN_MODULES +++ b/src/AOs/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Nuclei +Nuclei Utils diff --git a/src/AOs/README.rst b/src/AOs/README.rst index 3b1675ba..b411cd20 100644 --- a/src/AOs/README.rst +++ b/src/AOs/README.rst @@ -39,7 +39,10 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Nuclei `_ +* `Utils `_ Documentation ============= @@ -144,5 +147,11 @@ Documentation Per convention, for P,D,F and G AOs, we take the index of the AO with the the corresponding power in the "X" axis +`n_pt_max_i_x `_ + Undocumented + +`n_pt_max_integrals `_ + Undocumented + diff --git a/src/MonoInts/dimensions.irp.f b/src/AOs/dimensions_integrals.irp.f similarity index 100% rename from src/MonoInts/dimensions.irp.f rename to src/AOs/dimensions_integrals.irp.f diff --git a/src/AOs/tree_dependancy.png b/src/AOs/tree_dependancy.png new file mode 100644 index 00000000..3bf621ee Binary files /dev/null and b/src/AOs/tree_dependancy.png differ diff --git a/src/Bielec_integrals/NEEDED_CHILDREN_MODULES b/src/Bielec_integrals/NEEDED_CHILDREN_MODULES deleted file mode 100644 index 24979463..00000000 --- a/src/Bielec_integrals/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -MonoInts Bitmask diff --git a/src/Bitmask/README.rst b/src/Bitmask/README.rst index 7eeac8c9..4c73d802 100644 --- a/src/Bitmask/README.rst +++ b/src/Bitmask/README.rst @@ -40,6 +40,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `MOs `_ Documentation 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/Bitmask/tree_dependancy.png b/src/Bitmask/tree_dependancy.png new file mode 100644 index 00000000..c931c9fa Binary files /dev/null and b/src/Bitmask/tree_dependancy.png differ diff --git a/src/CAS_SD/EZFIO.cfg b/src/CAS_SD/EZFIO.cfg index 7f27f95a..6fb8c8a0 100644 --- a/src/CAS_SD/EZFIO.cfg +++ b/src/CAS_SD/EZFIO.cfg @@ -1,29 +1,3 @@ -[N_det_max_cas_sd] -type: Det_number_max -doc: Max number of determinants in the wave function -interface: input -default: 10000 - -[do_pt2_end] -type: logical -doc: If true, compute the PT2 at the end of the selection -interface: input -default: True - -[PT2_max] -type: PT2_energy -doc: The selection process stops when the largest PT2 (for all the state is lower - than pt2_max in absolute value -interface: input -default: 0.0001 - -[var_pt2_ratio] -type: Normalized_float -doc: The selection process stops when the energy ratio variational/(variational+PT2) - is equal to var_pt2_ratio -interface: input -default: 0.75 - [energy] type: double precision doc: "Calculated CAS-SD energy" diff --git a/src/CAS_SD/README.rst b/src/CAS_SD/README.rst index 2e56e56e..2a533523 100644 --- a/src/CAS_SD/README.rst +++ b/src/CAS_SD/README.rst @@ -24,6 +24,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Perturbation `_ * `Selectors_full `_ * `Generators_CAS `_ diff --git a/src/CAS_SD/cas_sd.irp.f b/src/CAS_SD/cas_sd.irp.f index e28334b2..c823489f 100644 --- a/src/CAS_SD/cas_sd.irp.f +++ b/src/CAS_SD/cas_sd.irp.f @@ -13,12 +13,12 @@ program full_ci N_det_old = 0 pt2 = 1.d0 diag_algorithm = "Lapack" - if (N_det > n_det_max_cas_sd) then + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_cas_sd + N_det = N_det_max soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -30,7 +30,7 @@ program full_ci print *, '-----' endif - do while (N_det < n_det_max_cas_sd.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) N_det_old = N_det call H_apply_CAS_SD(pt2, norm_pert, H_pert_diag, N_st) @@ -38,10 +38,10 @@ program full_ci PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > n_det_max_cas_sd) then + if (N_det > N_det_max) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_cas_sd + N_det = N_det_max soft_touch N_det psi_det psi_coef endif call diagonalize_CI diff --git a/src/CAS_SD/cas_sd_selected.irp.f b/src/CAS_SD/cas_sd_selected.irp.f index c3abadf4..08cfcf41 100644 --- a/src/CAS_SD/cas_sd_selected.irp.f +++ b/src/CAS_SD/cas_sd_selected.irp.f @@ -12,12 +12,12 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" - if (N_det > n_det_max_cas_sd) then + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_cas_sd + N_det = N_det_max soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -29,17 +29,17 @@ program full_ci print *, '-----' endif - do while (N_det < n_det_max_cas_sd.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > n_det_max_cas_sd) then + if (N_det > N_det_max) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_cas_sd + N_det = N_det_max soft_touch N_det psi_det psi_coef endif call diagonalize_CI diff --git a/src/CAS_SD/tree_dependancy.png b/src/CAS_SD/tree_dependancy.png new file mode 100644 index 00000000..0942fadc Binary files /dev/null and b/src/CAS_SD/tree_dependancy.png differ diff --git a/src/CID/README.rst b/src/CID/README.rst index 5d2fa851..b74099b5 100644 --- a/src/CID/README.rst +++ b/src/CID/README.rst @@ -15,6 +15,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Selectors_full `_ * `SingleRefMethod `_ diff --git a/src/CID/tree_dependancy.png b/src/CID/tree_dependancy.png new file mode 100644 index 00000000..5e421085 Binary files /dev/null and b/src/CID/tree_dependancy.png differ diff --git a/src/CID_SC2_selected/README.rst b/src/CID_SC2_selected/README.rst index ec9e7a3f..8227a14f 100644 --- a/src/CID_SC2_selected/README.rst +++ b/src/CID_SC2_selected/README.rst @@ -19,5 +19,7 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `CID_selected `_ diff --git a/src/CID_SC2_selected/tree_dependancy.png b/src/CID_SC2_selected/tree_dependancy.png new file mode 100644 index 00000000..86e753a4 Binary files /dev/null and b/src/CID_SC2_selected/tree_dependancy.png differ diff --git a/src/CID_selected/README.rst b/src/CID_selected/README.rst index 50cc701f..5ef3b185 100644 --- a/src/CID_selected/README.rst +++ b/src/CID_selected/README.rst @@ -22,6 +22,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Perturbation `_ * `CID `_ diff --git a/src/CID_selected/tree_dependancy.png b/src/CID_selected/tree_dependancy.png new file mode 100644 index 00000000..f8ad68e3 Binary files /dev/null and b/src/CID_selected/tree_dependancy.png differ diff --git a/src/CIS/README.rst b/src/CIS/README.rst index e3b2478a..f48ea5c1 100644 --- a/src/CIS/README.rst +++ b/src/CIS/README.rst @@ -31,6 +31,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Selectors_full `_ * `SingleRefMethod `_ diff --git a/src/CIS/tree_dependancy.png b/src/CIS/tree_dependancy.png new file mode 100644 index 00000000..bf88fde6 Binary files /dev/null and b/src/CIS/tree_dependancy.png differ diff --git a/src/CISD/README.rst b/src/CISD/README.rst index 68ab4cfb..b052b71f 100644 --- a/src/CISD/README.rst +++ b/src/CISD/README.rst @@ -15,6 +15,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Selectors_full `_ * `SingleRefMethod `_ diff --git a/src/CISD/tree_dependancy.png b/src/CISD/tree_dependancy.png new file mode 100644 index 00000000..2d040122 Binary files /dev/null and b/src/CISD/tree_dependancy.png differ diff --git a/src/CISD_SC2_selected/EZFIO.cfg b/src/CISD_SC2_selected/EZFIO.cfg index b76eac5d..3b49c52b 100644 --- a/src/CISD_SC2_selected/EZFIO.cfg +++ b/src/CISD_SC2_selected/EZFIO.cfg @@ -1,22 +1,3 @@ -[N_det_max_cisd_sc2] -type: Det_number_max -doc: Max number of determinants -interface: input -default: 10000 - -[do_pt2_end] -type: logical -doc: If true, compute the PT2 at the end of the selection -interface: input -default: True - -[PT2_max] -type: PT2_energy -doc: The selection process stops when the largest PT2 (for all the states) is lower - than abs(pt2_max) -interface: input -default: 0.0001 - [energy] type: double precision doc: Calculated CISD_SC2 energy of ground_state diff --git a/src/CISD_SC2_selected/README.rst b/src/CISD_SC2_selected/README.rst index b4f4fac1..40bf3bfc 100644 --- a/src/CISD_SC2_selected/README.rst +++ b/src/CISD_SC2_selected/README.rst @@ -19,5 +19,7 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `CISD_selected `_ diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f index 55234fcb..f4c96f54 100644 --- a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -15,12 +15,12 @@ program cisd_sc2_selected E_old(1) = HF_energy davidson_threshold = 1.d-10 - if (N_det > n_det_max_cisd_sc2) then + if (N_det > N_det_max) then call diagonalize_CI_SC2 call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_cisd_sc2 + N_det = N_det_max soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -34,7 +34,7 @@ program cisd_sc2_selected integer :: i_count i_count = 0 - do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) print*,'----' print*,'' call H_apply_SC2_selected(pt2, norm_pert, H_pert_diag, N_st) @@ -70,7 +70,7 @@ program cisd_sc2_selected call ezfio_set_full_ci_energy(CI_SC2_energy(1)) enddo - N_det = min(n_det_max_cisd_sc2,N_det) + N_det = min(N_det_max,N_det) davidson_threshold = 1.d-10 touch N_det psi_det psi_coef davidson_threshold davidson_criterion call diagonalize_CI_SC2 diff --git a/src/CISD_SC2_selected/tree_dependancy.png b/src/CISD_SC2_selected/tree_dependancy.png new file mode 100644 index 00000000..ac57a98d Binary files /dev/null and b/src/CISD_SC2_selected/tree_dependancy.png differ diff --git a/src/CISD_selected/README.rst b/src/CISD_selected/README.rst index a72c5a21..d2730f2b 100644 --- a/src/CISD_selected/README.rst +++ b/src/CISD_selected/README.rst @@ -14,12 +14,6 @@ Documentation `cisd `_ Undocumented -`n_det_max_cisd `_ - Get n_det_max_cisd from EZFIO file - -`pt2_max `_ - Get pt2_max from EZFIO file - Needed Modules @@ -28,6 +22,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Perturbation `_ * `CISD `_ diff --git a/src/CISD_selected/cisd.ezfio_config b/src/CISD_selected/cisd.ezfio_config deleted file mode 100644 index be3854b1..00000000 --- a/src/CISD_selected/cisd.ezfio_config +++ /dev/null @@ -1,3 +0,0 @@ -cisd_selected - n_det_max_cisd integer - pt2_max double precision diff --git a/src/CISD_selected/cisd_selection.irp.f b/src/CISD_selected/cisd_selection.irp.f index c74d0f3a..b05e9ea4 100644 --- a/src/CISD_selected/cisd_selection.irp.f +++ b/src/CISD_selected/cisd_selection.irp.f @@ -18,7 +18,7 @@ program cisd print *, 'E = ', CI_energy(i) enddo E_old = CI_energy - do while (maxval(abs(pt2(1:N_st))) > pt2_max.and.n_det < n_det_max_cisd) + do while (maxval(abs(pt2(1:N_st))) > pt2_max.and.n_det < N_det_max) print*,'----' print*,'' call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) @@ -38,7 +38,7 @@ program cisd exit endif enddo - N_det = min(N_det,n_det_max_cisd) + N_det = min(N_det,N_det_max) touch N_det psi_det psi_coef call diagonalize_CI deallocate(pt2,norm_pert,H_pert_diag) diff --git a/src/CISD_selected/options.irp.f b/src/CISD_selected/options.irp.f deleted file mode 100644 index 53aa0ef6..00000000 --- a/src/CISD_selected/options.irp.f +++ /dev/null @@ -1,34 +0,0 @@ - BEGIN_PROVIDER [ integer, n_det_max_cisd ] - implicit none - BEGIN_DOC - ! Get n_det_max_cisd from EZFIO file - END_DOC - logical :: has_n_det_max_cisd - PROVIDE ezfio_filename - call ezfio_has_cisd_selected_n_det_max_cisd(has_n_det_max_cisd) - if (has_n_det_max_cisd) then - call ezfio_get_cisd_selected_n_det_max_cisd(n_det_max_cisd) - else - n_det_max_cisd = 30000 - call ezfio_set_cisd_selected_n_det_max_cisd(n_det_max_cisd) - endif - print*,'n_det_max_cisd = ',n_det_max_cisd - END_PROVIDER - - BEGIN_PROVIDER [ double precision , pt2_max ] - implicit none - BEGIN_DOC - ! Get pt2_max from EZFIO file - END_DOC - logical :: has_pt2_max - PROVIDE ezfio_filename - call ezfio_has_cisd_selected_pt2_max(has_pt2_max) - if (has_pt2_max) then - call ezfio_get_cisd_selected_pt2_max(pt2_max) - else - pt2_max = 1.d-9 - call ezfio_set_cisd_selected_pt2_max(pt2_max) - endif - print*,'pt2_max = ',pt2_max - END_PROVIDER - diff --git a/src/CISD_selected/tree_dependancy.png b/src/CISD_selected/tree_dependancy.png new file mode 100644 index 00000000..c33e2e66 Binary files /dev/null and b/src/CISD_selected/tree_dependancy.png differ diff --git a/src/DDCI_selected/README.rst b/src/DDCI_selected/README.rst index 9cfdbefa..d631d53a 100644 --- a/src/DDCI_selected/README.rst +++ b/src/DDCI_selected/README.rst @@ -19,6 +19,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Perturbation `_ * `Selectors_full `_ * `Generators_CAS `_ diff --git a/src/DDCI_selected/ddci.irp.f b/src/DDCI_selected/ddci.irp.f index 3f3b6159..1a44e3b1 100644 --- a/src/DDCI_selected/ddci.irp.f +++ b/src/DDCI_selected/ddci.irp.f @@ -11,12 +11,12 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" - if (N_det > n_det_max_fci) then + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = N_det_max soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -28,17 +28,17 @@ program full_ci print *, '-----' endif - do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) call H_apply_DDCI_selection(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > n_det_max_fci) then + if (N_det > N_det_max) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = N_det_max soft_touch N_det psi_det psi_coef endif call diagonalize_CI diff --git a/src/DDCI_selected/options.irp.f b/src/DDCI_selected/options.irp.f deleted file mode 100644 index 553fbe9f..00000000 --- a/src/DDCI_selected/options.irp.f +++ /dev/null @@ -1,32 +0,0 @@ -BEGIN_SHELL [ /usr/bin/python ] -from ezfio_with_default import EZFIO_Provider -T = EZFIO_Provider() -T.set_type ( "integer" ) -T.set_name ( "N_det_max_fci" ) -T.set_doc ( "Max number of determinants in the wave function" ) -T.set_ezfio_dir ( "full_ci" ) -T.set_ezfio_name( "N_det_max_fci" ) -T.set_output ( "output_full_ci" ) -print T - -T.set_type ( "logical" ) -T.set_name ( "do_pt2_end" ) -T.set_doc ( "If true, compute the PT2 at the end of the selection" ) -T.set_ezfio_name( "do_pt2_end" ) -print T - -T.set_type ( "double precision" ) -T.set_name ( "pt2_max" ) -T.set_doc ( """The selection process stops when the largest PT2 (for all the states) -is lower than pt2_max in absolute value""" ) -T.set_ezfio_name( "pt2_max" ) -print T - -T.set_type ( "double precision" ) -T.set_name ( "var_pt2_ratio" ) -T.set_doc ( """The selection process stops when the energy ratio variational/(variational+PT2) -is equal to var_pt2_ratio""" ) -T.set_ezfio_name( "var_pt2_ratio" ) -print T -END_SHELL - diff --git a/src/DDCI_selected/tree_dependancy.png b/src/DDCI_selected/tree_dependancy.png new file mode 100644 index 00000000..9975ff6c Binary files /dev/null and b/src/DDCI_selected/tree_dependancy.png differ diff --git a/src/Bielec_integrals/ASSUMPTIONS.rst b/src/DensityFit/ASSUMPTIONS.rst similarity index 100% rename from src/Bielec_integrals/ASSUMPTIONS.rst rename to src/DensityFit/ASSUMPTIONS.rst diff --git a/src/Bielec_integrals/Makefile b/src/DensityFit/Makefile similarity index 100% rename from src/Bielec_integrals/Makefile rename to src/DensityFit/Makefile 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/DensityFit/tree_dependancy.png b/src/DensityFit/tree_dependancy.png new file mode 100644 index 00000000..62bad0ff Binary files /dev/null and b/src/DensityFit/tree_dependancy.png differ diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index eab4ea0e..4a2270c7 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -1,15 +1,27 @@ +[N_det_max] +type: Det_number_max +doc: Max number of determinants in the wave function +interface: input +default: 10000 + +[N_det_max_property] +type: Det_number_max +doc: Max number of determinants in the wave function when you select for a given property +interface: input +default: 10000 + +[N_det_max_jacobi] +type: Det_number_max +doc: Maximum number of determinants diagonalized by Jacobi +interface: input +default: 1000 + [N_states] type: States_number doc: Number of states to consider interface: input default: 1 -[N_det_max_jacobi] -type: Strictly_positive_int -doc: Maximum number of determinants diagonalized by Jacobi -interface: input -default: 1000 - [read_wf] type: logical doc: If true, read the wave function from the EZFIO file @@ -40,9 +52,6 @@ doc: Thresholds on selectors (fraction of the norm) interface: input default: 0.999 - -# Only create the ezfio_config, (no Input_* and no PROVIDER) - [n_states_diag] type: integer doc: n_states_diag @@ -72,13 +81,13 @@ type: integer interface: OCaml doc: psi_coef type: double precision -size: (determinants_n_det,determinants_n_states) +size: (determinants.n_det,determinants.n_states) [psi_det] interface: OCaml doc: psi_det type: integer*8 -size: (determinants_n_int*determinants_bit_kind/8,2,determinants_n_det) +size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det) [det_num] interface: OCaml @@ -89,13 +98,13 @@ type: integer interface: OCaml doc: det_occ type: integer -size: (electrons_elec_alpha_num,determinants_det_num,2) +size: (electrons.elec_alpha_num,determinants.det_num,2) [det_coef] interface: OCaml doc: det_coef type: double precision -size: (determinants_det_num) +size: (determinants.det_num) [expected_s2] interface: OCaml diff --git a/src/Determinants/NEEDED_CHILDREN_MODULES b/src/Determinants/NEEDED_CHILDREN_MODULES index afc397de..5505ce78 100644 --- a/src/Determinants/NEEDED_CHILDREN_MODULES +++ b/src/Determinants/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Bielec_integrals +Integrals_Monoelec Integrals_Bielec \ No newline at end of file diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index e9307357..a9455783 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -32,7 +32,10 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* `Bielec_integrals `_ +.. image:: tree_dependancy.png + +* `Integrals_Monoelec `_ +* `Integrals_Bielec `_ Documentation ============= @@ -529,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/Determinants/tree_dependancy.png b/src/Determinants/tree_dependancy.png new file mode 100644 index 00000000..1227a472 Binary files /dev/null and b/src/Determinants/tree_dependancy.png differ diff --git a/src/Electrons/NEEDED_CHILDREN_MODULES b/src/Electrons/NEEDED_CHILDREN_MODULES index 83260f86..d9eaf57c 100644 --- a/src/Electrons/NEEDED_CHILDREN_MODULES +++ b/src/Electrons/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Output +Ezfio_files diff --git a/src/Electrons/README.rst b/src/Electrons/README.rst index e71a552d..da93bf74 100644 --- a/src/Electrons/README.rst +++ b/src/Electrons/README.rst @@ -24,7 +24,9 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* `Output `_ +.. image:: tree_dependancy.png + +* `Ezfio_files `_ Documentation ============= diff --git a/src/Electrons/tree_dependancy.png b/src/Electrons/tree_dependancy.png new file mode 100644 index 00000000..b90a1f83 Binary files /dev/null and b/src/Electrons/tree_dependancy.png differ diff --git a/src/Ezfio_files/README.rst b/src/Ezfio_files/README.rst index e0ef23da..4a4ce480 100644 --- a/src/Ezfio_files/README.rst +++ b/src/Ezfio_files/README.rst @@ -28,6 +28,24 @@ Documentation 'x' : READ/WRITE, FORMATTED .br +`output_cpu_time_0 `_ + Initial CPU and wall times when printing in the output files + +`output_wall_time_0 `_ + Initial CPU and wall times when printing in the output files + +`write_bool `_ + Write an logical value in output + +`write_double `_ + Write a double precision value in output + +`write_int `_ + Write an integer value in output + +`write_time `_ + Write a time stamp in the output for chronological reconstruction + diff --git a/src/Output/output.irp.f b/src/Ezfio_files/output.irp.f similarity index 96% rename from src/Output/output.irp.f rename to src/Ezfio_files/output.irp.f index 85f5cc0e..1a45e865 100644 --- a/src/Output/output.irp.f +++ b/src/Ezfio_files/output.irp.f @@ -20,8 +20,8 @@ BEGIN_SHELL [ /bin/bash ] ! Output file for $NAME END_DOC PROVIDE output_wall_time_0 output_cpu_time_0 ezfio_filename - integer :: getUnitAndOpen - call ezfio_set_output_empty(.False.) +! integer :: getUnitAndOpen +! call ezfio_set_output_empty(.False.) IRP_IF COARRAY if (this_image() == 1) then output_$NAME = 6 !getUnitAndOpen(trim(ezfio_filename)//'/output/'//'$NAME.rst','a') diff --git a/src/Ezfio_files/tree_dependancy.png b/src/Ezfio_files/tree_dependancy.png new file mode 100644 index 00000000..48f53991 Binary files /dev/null and b/src/Ezfio_files/tree_dependancy.png differ diff --git a/src/FCIdump/README.rst b/src/FCIdump/README.rst index 8a467b16..1f25ca6e 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 file. -`fcidump `_ - Undocumented - Needed Modules @@ -21,5 +18,7 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Determinants `_ diff --git a/src/FCIdump/tree_dependancy.png b/src/FCIdump/tree_dependancy.png new file mode 100644 index 00000000..5563cfce Binary files /dev/null and b/src/FCIdump/tree_dependancy.png differ diff --git a/src/Full_CI/EZFIO.cfg b/src/Full_CI/EZFIO.cfg index 89679d1a..37f25eda 100644 --- a/src/Full_CI/EZFIO.cfg +++ b/src/Full_CI/EZFIO.cfg @@ -1,35 +1,3 @@ -[N_det_max_fci] -type: Det_number_max -doc: Max number of determinants in the wave function -interface: input -default: 10000 - -[N_det_max_fci_property] -type: Det_number_max -doc: Max number of determinants in the wave function when you select for a given property -interface: input -default: 10000 - -[do_pt2_end] -type: logical -doc: If true, compute the PT2 at the end of the selection -interface: input -default: True - -[PT2_max] -type: PT2_energy -doc: The selection process stops when the largest PT2 (for all the state) is lower - than pt2_max in absolute value -interface: input -default: 0.0001 - -[var_pt2_ratio] -type: Normalized_float -doc: The selection process stops when the energy ratio variational/(variational+PT2) - is equal to var_pt2_ratio -interface: input -default: 0.75 - [energy] type: double precision doc: Calculated Selected FCI energy diff --git a/src/Full_CI/README.rst b/src/Full_CI/README.rst index 51bb05d2..1186feb1 100644 --- a/src/Full_CI/README.rst +++ b/src/Full_CI/README.rst @@ -24,6 +24,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Perturbation `_ * `Selectors_full `_ * `Generators_full `_ diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f index 594cc474..e2a9700e 100644 --- a/src/Full_CI/full_ci.irp.f +++ b/src/Full_CI/full_ci.irp.f @@ -11,12 +11,12 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" - if (N_det > n_det_max_fci) then + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = N_det_max soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -38,7 +38,7 @@ program full_ci integer :: n_det_before print*,'Beginning the selection ...' - do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) n_det_before = N_det call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) @@ -46,10 +46,10 @@ program full_ci PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > n_det_max_fci) then + if (N_det > N_det_max) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = N_det_max soft_touch N_det psi_det psi_coef endif call diagonalize_CI @@ -68,7 +68,7 @@ program full_ci exit endif enddo - N_det = min(n_det_max_fci,N_det) + N_det = min(N_det_max,N_det) touch N_det psi_det psi_coef call diagonalize_CI if(do_pt2_end)then diff --git a/src/Full_CI/full_ci_no_skip.irp.f b/src/Full_CI/full_ci_no_skip.irp.f index aa84fb9d..73958bf9 100644 --- a/src/Full_CI/full_ci_no_skip.irp.f +++ b/src/Full_CI/full_ci_no_skip.irp.f @@ -11,12 +11,12 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" - if (N_det > n_det_max_fci) then + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = N_det_max soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -38,7 +38,7 @@ program full_ci integer :: n_det_before print*,'Beginning the selection ...' - do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) n_det_before = N_det call H_apply_FCI_no_skip(pt2, norm_pert, H_pert_diag, N_st) @@ -46,10 +46,10 @@ program full_ci PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > n_det_max_fci) then + if (N_det > N_det_max) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = N_det_max soft_touch N_det psi_det psi_coef endif call diagonalize_CI @@ -68,7 +68,7 @@ program full_ci exit endif enddo - N_det = min(n_det_max_fci,N_det) + N_det = min(N_det_max,N_det) touch N_det psi_det psi_coef call diagonalize_CI if(do_pt2_end)then diff --git a/src/Full_CI/tree_dependancy.png b/src/Full_CI/tree_dependancy.png new file mode 100644 index 00000000..766eec3a Binary files /dev/null and b/src/Full_CI/tree_dependancy.png differ diff --git a/src/Generators_CAS/README.rst b/src/Generators_CAS/README.rst index b729212c..e489a995 100644 --- a/src/Generators_CAS/README.rst +++ b/src/Generators_CAS/README.rst @@ -43,5 +43,7 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Determinants `_ diff --git a/src/Generators_CAS/tree_dependancy.png b/src/Generators_CAS/tree_dependancy.png new file mode 100644 index 00000000..49fcefcd Binary files /dev/null and b/src/Generators_CAS/tree_dependancy.png differ diff --git a/src/Generators_full/README.rst b/src/Generators_full/README.rst index e558f48b..b9260a55 100644 --- a/src/Generators_full/README.rst +++ b/src/Generators_full/README.rst @@ -40,6 +40,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Determinants `_ * `Hartree_Fock `_ diff --git a/src/Generators_full/tree_dependancy.png b/src/Generators_full/tree_dependancy.png new file mode 100644 index 00000000..a489df79 Binary files /dev/null and b/src/Generators_full/tree_dependancy.png differ diff --git a/src/Hartree_Fock/NEEDED_CHILDREN_MODULES b/src/Hartree_Fock/NEEDED_CHILDREN_MODULES index b779faec..784cb0fb 100644 --- a/src/Hartree_Fock/NEEDED_CHILDREN_MODULES +++ b/src/Hartree_Fock/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Bielec_integrals MOGuess +Integrals_Bielec MOGuess diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 2b470e9f..70544122 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -10,7 +10,9 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* `Bielec_integrals `_ +.. image:: tree_dependancy.png + +* `Integrals_Bielec `_ * `MOGuess `_ Documentation diff --git a/src/Hartree_Fock/tree_dependancy.png b/src/Hartree_Fock/tree_dependancy.png new file mode 100644 index 00000000..38a98b49 Binary files /dev/null and b/src/Hartree_Fock/tree_dependancy.png differ diff --git a/src/MonoInts/ASSUMPTIONS.rst b/src/Integrals_Bielec/ASSUMPTIONS.rst similarity index 100% rename from src/MonoInts/ASSUMPTIONS.rst rename to src/Integrals_Bielec/ASSUMPTIONS.rst diff --git a/src/Bielec_integrals/EZFIO.cfg b/src/Integrals_Bielec/EZFIO.cfg similarity index 100% rename from src/Bielec_integrals/EZFIO.cfg rename to src/Integrals_Bielec/EZFIO.cfg diff --git a/src/Output/Makefile b/src/Integrals_Bielec/Makefile similarity index 100% rename from src/Output/Makefile rename to src/Integrals_Bielec/Makefile diff --git a/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES b/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..62fd5829 --- /dev/null +++ b/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Pseudo MOs Bitmask diff --git a/src/Bielec_integrals/README.rst b/src/Integrals_Bielec/README.rst similarity index 71% rename from src/Bielec_integrals/README.rst rename to src/Integrals_Bielec/README.rst index d4bbb16b..7ec077b5 100644 --- a/src/Bielec_integrals/README.rst +++ b/src/Integrals_Bielec/README.rst @@ -16,7 +16,10 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* `MonoInts `_ +.. image:: tree_dependancy.png + +* `Pseudo `_ +* `MOs `_ * `Bitmask `_ Documentation @@ -25,192 +28,192 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`ao_bielec_integral `_ +`ao_bielec_integral `_ integral of the AO basis or (ij|kl) i(r1) j(r1) 1/r12 k(r2) l(r2) -`ao_bielec_integral_schwartz `_ +`ao_bielec_integral_schwartz `_ Needed to compute Schwartz inequalities -`ao_bielec_integral_schwartz_accel `_ +`ao_bielec_integral_schwartz_accel `_ integral of the AO basis or (ij|kl) i(r1) j(r1) 1/r12 k(r2) l(r2) -`ao_bielec_integrals_in_map `_ +`ao_bielec_integrals_in_map `_ Map of Atomic integrals i(r1) j(r2) 1/r12 k(r1) l(r2) -`ao_l4 `_ +`ao_l4 `_ Computes the product of l values of i,j,k,and l -`compute_ao_bielec_integrals `_ +`compute_ao_bielec_integrals `_ Compute AO 1/r12 integrals for all i and fixed j,k,l -`eri `_ +`eri `_ ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) -`general_primitive_integral `_ +`general_primitive_integral `_ Computes the integral where p,q,r,s are Gaussian primitives -`give_polynom_mult_center_x `_ +`give_polynom_mult_center_x `_ subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw : I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) -`i_x1_new `_ +`i_x1_new `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult `_ +`i_x1_pol_mult `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a1 `_ +`i_x1_pol_mult_a1 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a2 `_ +`i_x1_pol_mult_a2 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_recurs `_ +`i_x1_pol_mult_recurs `_ recursive function involved in the bielectronic integral -`i_x2_new `_ +`i_x2_new `_ recursive function involved in the bielectronic integral -`i_x2_pol_mult `_ +`i_x2_pol_mult `_ recursive function involved in the bielectronic integral -`integrale_new `_ +`integrale_new `_ calculate the integral of the polynom :: I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q) between ( 0 ; 1) -`n_pt_sup `_ +`n_pt_sup `_ Returns the upper boundary of the degree of the polynomial involved in the bielctronic integral : Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z) -`gauleg `_ +`gauleg `_ Gauss-Legendre -`gauleg_t2 `_ +`gauleg_t2 `_ t_w(i,1,k) = w(i) t_w(i,2,k) = t(i) -`gauleg_w `_ +`gauleg_w `_ t_w(i,1,k) = w(i) t_w(i,2,k) = t(i) -`n_pt_max_integrals_16 `_ +`n_pt_max_integrals_16 `_ Aligned n_pt_max_integrals -`ao_integrals_map `_ +`ao_integrals_map `_ AO integrals -`bielec_integrals_index `_ +`bielec_integrals_index `_ Undocumented -`bielec_integrals_index_reverse `_ +`bielec_integrals_index_reverse `_ Undocumented -`clear_ao_map `_ +`clear_ao_map `_ Frees the memory of the AO map -`clear_mo_map `_ +`clear_mo_map `_ Frees the memory of the MO map -`get_ao_bielec_integral `_ +`get_ao_bielec_integral `_ Gets one AO bi-electronic integral from the AO map -`get_ao_bielec_integrals `_ +`get_ao_bielec_integrals `_ Gets multiple AO bi-electronic integral from the AO map . All i are retrieved for j,k,l fixed. -`get_ao_bielec_integrals_non_zero `_ +`get_ao_bielec_integrals_non_zero `_ Gets multiple AO bi-electronic integral from the AO map . All non-zero i are retrieved for j,k,l fixed. -`get_ao_map_size `_ +`get_ao_map_size `_ Returns the number of elements in the AO map -`get_mo_bielec_integral `_ +`get_mo_bielec_integral `_ Returns one integral in the MO basis -`get_mo_bielec_integrals `_ +`get_mo_bielec_integrals `_ Returns multiple integrals in the MO basis, all i for j,k,l fixed. -`get_mo_bielec_integrals_existing_ik `_ +`get_mo_bielec_integrals_existing_ik `_ Returns multiple integrals in the MO basis, all i(1)j(1) 1/r12 k(2)l(2) i for j,k,l fixed. -`get_mo_map_size `_ +`get_mo_map_size `_ Return the number of elements in the MO map -`insert_into_ao_integrals_map `_ +`insert_into_ao_integrals_map `_ Create new entry into AO map -`insert_into_mo_integrals_map `_ +`insert_into_mo_integrals_map `_ Create new entry into MO map, or accumulate in an existing entry -`mo_bielec_integral `_ +`mo_bielec_integral `_ Returns one integral in the MO basis -`mo_integrals_map `_ +`mo_integrals_map `_ MO integrals -`add_integrals_to_map `_ +`add_integrals_to_map `_ Adds integrals to tha MO map according to some bitmask -`mo_bielec_integral_jj `_ +`mo_bielec_integral_jj `_ mo_bielec_integral_jj(i,j) = J_ij mo_bielec_integral_jj_exchange(i,j) = K_ij mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_anti `_ +`mo_bielec_integral_jj_anti `_ mo_bielec_integral_jj(i,j) = J_ij mo_bielec_integral_jj_exchange(i,j) = K_ij mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_anti_from_ao `_ +`mo_bielec_integral_jj_anti_from_ao `_ mo_bielec_integral_jj_from_ao(i,j) = J_ij mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_exchange `_ +`mo_bielec_integral_jj_exchange `_ mo_bielec_integral_jj(i,j) = J_ij mo_bielec_integral_jj_exchange(i,j) = K_ij mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_exchange_from_ao `_ +`mo_bielec_integral_jj_exchange_from_ao `_ mo_bielec_integral_jj_from_ao(i,j) = J_ij mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_from_ao `_ +`mo_bielec_integral_jj_from_ao `_ mo_bielec_integral_jj_from_ao(i,j) = J_ij mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij -`mo_bielec_integrals_in_map `_ +`mo_bielec_integrals_in_map `_ If True, the map of MO bielectronic integrals is provided -`mo_bielec_integrals_index `_ +`mo_bielec_integrals_index `_ Computes an unique index for i,j,k,l integrals -`read_ao_integrals `_ +`read_ao_integrals `_ One level of abstraction for disk_access_ao_integrals and disk_access_mo_integrals -`read_mo_integrals `_ +`read_mo_integrals `_ One level of abstraction for disk_access_ao_integrals and disk_access_mo_integrals -`write_ao_integrals `_ +`write_ao_integrals `_ One level of abstraction for disk_access_ao_integrals and disk_access_mo_integrals -`write_mo_integrals `_ +`write_mo_integrals `_ One level of abstraction for disk_access_ao_integrals and disk_access_mo_integrals diff --git a/src/Bielec_integrals/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f similarity index 97% rename from src/Bielec_integrals/ao_bi_integrals.irp.f rename to src/Integrals_Bielec/ao_bi_integrals.irp.f index 77f9a098..a9fc3eb1 100644 --- a/src/Bielec_integrals/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -349,7 +349,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integer :: load_ao_integrals print*,'Reading the AO integrals' if (load_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') == 0) then - write(output_bielec_integrals,*) 'AO integrals provided' + print*, 'AO integrals provided' ao_bielec_integrals_in_map = .True. return endif @@ -367,7 +367,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] PROVIDE progress_bar call omp_init_lock(lock) lmax = ao_num*(ao_num+1)/2 - write(output_bielec_integrals,*) 'Providing the AO integrals' + print*, 'Providing the AO integrals' call wall_time(wall_0) call wall_time(wall_1) call cpu_time(cpu_1) @@ -378,7 +378,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] !$OMP DEFAULT(NONE) & !$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, & !$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, & - !$OMP ao_overlap_abs,ao_overlap,output_bielec_integrals,abort_here, & + !$OMP ao_overlap_abs,ao_overlap,abort_here, & !$OMP wall_0,progress_bar,progress_value) allocate(buffer_i(size_buffer)) @@ -445,7 +445,7 @@ IRP_ENDIF if (thread_num == 0) then if (wall_2 - wall_0 > 1.d0) then wall_0 = wall_2 - write(output_bielec_integrals,*) 100.*float(kk)/float(lmax), '% in ', & + print*, 100.*float(kk)/float(lmax), '% in ', & wall_2-wall_1, 's', map_mb(ao_integrals_map) ,'MB' progress_value = dble(map_mb(ao_integrals_map)) endif @@ -462,32 +462,30 @@ IRP_ENDIF stop 'Aborting in AO integrals calculation' endif IRP_IF COARRAY - write(output_bielec_integrals,*) 'Communicating the map' + print*, 'Communicating the map' call communicate_ao_integrals() IRP_ENDIF COARRAY - write(output_bielec_integrals,*) 'Sorting the map' + print*, 'Sorting the map' call map_sort(ao_integrals_map) call cpu_time(cpu_2) call wall_time(wall_2) integer(map_size_kind) :: get_ao_map_size, ao_map_size ao_map_size = get_ao_map_size() - write(output_bielec_integrals,*) 'AO integrals provided:' - write(output_bielec_integrals,*) ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' - write(output_bielec_integrals,*) ' Number of AO integrals :', ao_map_size - write(output_bielec_integrals,*) ' cpu time :',cpu_2 - cpu_1, 's' - write(output_bielec_integrals,*) ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' + print*, 'AO integrals provided:' + print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' + print*, ' Number of AO integrals :', ao_map_size + print*, ' cpu time :',cpu_2 - cpu_1, 's' + print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' ao_bielec_integrals_in_map = .True. if (write_ao_integrals) then call dump_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') - call ezfio_set_bielec_integrals_disk_access_ao_integrals(.True.) + call ezfio_set_integrals_bielec_disk_access_ao_integrals(.True.) endif END_PROVIDER - - - + BEGIN_PROVIDER [ double precision, ao_bielec_integral_schwartz,(ao_num,ao_num) ] implicit none BEGIN_DOC diff --git a/src/Bielec_integrals/gauss_legendre.irp.f b/src/Integrals_Bielec/gauss_legendre.irp.f similarity index 100% rename from src/Bielec_integrals/gauss_legendre.irp.f rename to src/Integrals_Bielec/gauss_legendre.irp.f diff --git a/src/Bielec_integrals/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f similarity index 99% rename from src/Bielec_integrals/map_integrals.irp.f rename to src/Integrals_Bielec/map_integrals.irp.f index 8fdff799..83950c5c 100644 --- a/src/Bielec_integrals/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -13,7 +13,7 @@ BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max) sze = key_max call map_init(ao_integrals_map,sze) - write(output_bielec_integrals,*) 'AO map initialized' + print*, 'AO map initialized' END_PROVIDER subroutine bielec_integrals_index(i,j,k,l,i1) @@ -244,7 +244,7 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_map ] call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max) sze = key_max call map_init(mo_integrals_map,sze) - write(output_bielec_integrals,*) 'MO map initialized' + print*, 'MO map initialized' END_PROVIDER subroutine insert_into_ao_integrals_map(n_integrals, & diff --git a/src/Bielec_integrals/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f similarity index 94% rename from src/Bielec_integrals/mo_bi_integrals.irp.f rename to src/Integrals_Bielec/mo_bi_integrals.irp.f index 4e4c70aa..bc724e9b 100644 --- a/src/Bielec_integrals/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -31,7 +31,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] integer :: load_mo_integrals print*,'Reading the MO integrals' if (load_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') == 0) then - write(output_bielec_integrals,*) 'MO integrals provided' + print*, 'MO integrals provided' return endif endif @@ -84,8 +84,8 @@ subroutine add_integrals_to_map(mask_ijkl) call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) size_buffer = min(ao_num*ao_num*ao_num,16000000) - write(output_bielec_integrals,*) 'Providing the molecular integrals ' - write(output_bielec_integrals,*) 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + print*, 'Providing the molecular integrals ' + print*, 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' call wall_time(wall_1) @@ -99,7 +99,7 @@ subroutine add_integrals_to_map(mask_ijkl) !$OMP wall_0,thread_num) & !$OMP DEFAULT(NONE) & !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,& - !$OMP mo_coef_transp,output_bielec_integrals, & + !$OMP mo_coef_transp, & !$OMP mo_coef_transp_is_built, list_ijkl, & !$OMP mo_coef_is_built, wall_1, abort_here, & !$OMP mo_coef,mo_integrals_threshold,ao_integrals_map,mo_integrals_map,progress_bar,progress_value) @@ -272,7 +272,7 @@ IRP_ENDIF if (thread_num == 0) then if (wall_2 - wall_0 > 1.d0) then wall_0 = wall_2 - write(output_bielec_integrals,*) 100.*float(l1)/float(ao_num), '% in ', & + print*, 100.*float(l1)/float(ao_num), '% in ', & wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' progress_value = dble(map_mb(mo_integrals_map)) @@ -291,7 +291,7 @@ IRP_ENDIF stop 'Aborting in MO integrals calculation' endif IRP_IF COARRAY - write(output_bielec_integrals,*) 'Communicating the map' + print*, 'Communicating the map' call communicate_mo_integrals() IRP_ENDIF call map_unique(mo_integrals_map) @@ -304,15 +304,15 @@ IRP_ENDIF deallocate(list_ijkl) - write(output_bielec_integrals,*)'Molecular integrals provided:' - write(output_bielec_integrals,*)' Size of MO map ', map_mb(mo_integrals_map) ,'MB' - write(output_bielec_integrals,*)' Number of MO integrals: ', mo_map_size - write(output_bielec_integrals,*)' cpu time :',cpu_2 - cpu_1, 's' - write(output_bielec_integrals,*)' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + print*,' Number of MO integrals: ', mo_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' if (write_mo_integrals) then call dump_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') - call ezfio_set_bielec_integrals_disk_access_mo_integrals(.True.) + call ezfio_set_integrals_bielec_disk_access_mo_integrals(.True.) endif diff --git a/src/Bielec_integrals/read_write.irp.f b/src/Integrals_Bielec/read_write.irp.f similarity index 100% rename from src/Bielec_integrals/read_write.irp.f rename to src/Integrals_Bielec/read_write.irp.f diff --git a/src/Integrals_Bielec/tree_dependancy.png b/src/Integrals_Bielec/tree_dependancy.png new file mode 100644 index 00000000..d7298f31 Binary files /dev/null and b/src/Integrals_Bielec/tree_dependancy.png differ diff --git a/src/Output/ASSUMPTIONS.rst b/src/Integrals_Monoelec/ASSUMPTIONS.rst similarity index 100% rename from src/Output/ASSUMPTIONS.rst rename to src/Integrals_Monoelec/ASSUMPTIONS.rst diff --git a/src/Integrals_Monoelec/Makefile b/src/Integrals_Monoelec/Makefile new file mode 100644 index 00000000..abe31fed --- /dev/null +++ b/src/Integrals_Monoelec/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=pseudopot.f90 +OBJ=IRPF90_temp/pseudopot.o + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/Integrals_Monoelec/NEEDED_CHILDREN_MODULES b/src/Integrals_Monoelec/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..31cf3c24 --- /dev/null +++ b/src/Integrals_Monoelec/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +MOs Pseudo diff --git a/src/MonoInts/README.rst b/src/Integrals_Monoelec/README.rst similarity index 65% rename from src/MonoInts/README.rst rename to src/Integrals_Monoelec/README.rst index ac4983a9..b7054f47 100644 --- a/src/MonoInts/README.rst +++ b/src/Integrals_Monoelec/README.rst @@ -4,8 +4,10 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `MOs `_ -* `Pseudo_integrals `_ +* `Pseudo `_ Documentation ============= @@ -13,205 +15,211 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`ao_mono_elec_integral `_ +`ao_mono_elec_integral `_ array of the mono electronic hamiltonian on the AOs basis : sum of the kinetic and nuclear electronic potential -`check_ortho `_ +`check_ortho `_ Undocumented -`do_print `_ +`do_print `_ Undocumented -`n_pt_max_i_x `_ - Undocumented - -`n_pt_max_integrals `_ - Undocumented - -`ao_deriv2_x `_ +`ao_deriv2_x `_ second derivatives matrix elements in the ao basis .. math:: .br {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle -`ao_deriv2_y `_ +`ao_deriv2_y `_ second derivatives matrix elements in the ao basis .. math:: .br {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle -`ao_deriv2_z `_ +`ao_deriv2_z `_ second derivatives matrix elements in the ao basis .. math:: .br {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle -`ao_kinetic_integral `_ +`ao_kinetic_integral `_ array of the priminitve basis kinetic integrals \langle \chi_i |\hat{T}| \chi_j \rangle -`mo_kinetic_integral `_ +`mo_kinetic_integral `_ Undocumented -`mo_mono_elec_integral `_ +`mo_mono_elec_integral `_ array of the mono electronic hamiltonian on the MOs basis : sum of the kinetic and nuclear electronic potential -`orthonormalize_mos `_ +`orthonormalize_mos `_ Undocumented -`ao_nucl_elec_integral `_ +`ao_nucl_elec_integral `_ interaction nuclear electron -`ao_nucl_elec_integral_per_atom `_ +`ao_nucl_elec_integral_per_atom `_ ao_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom -`give_polynom_mult_center_mono_elec `_ +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented -`nai_pol_mult `_ +`nai_pol_mult `_ Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented -`mo_nucl_elec_integral `_ +`ao_pseudo_integral `_ + Pseudo-potential + +`ao_pseudo_integral_local `_ + Local pseudo-potential + +`ao_pseudo_integral_non_local `_ + Local pseudo-potential + +`mo_nucl_elec_integral `_ interaction nuclear electron on the MO basis -`mo_nucl_elec_integral_per_atom `_ +`mo_nucl_elec_integral_per_atom `_ mo_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom -`save_ortho_mos `_ +`mo_pseudo_integral `_ + interaction nuclear electron on the MO basis + +`save_ortho_mos `_ Undocumented -`ao_deriv_1_x `_ +`ao_deriv_1_x `_ array of the integrals of AO_i * d/dx AO_j array of the integrals of AO_i * d/dy AO_j array of the integrals of AO_i * d/dz AO_j -`ao_deriv_1_y `_ +`ao_deriv_1_y `_ array of the integrals of AO_i * d/dx AO_j array of the integrals of AO_i * d/dy AO_j array of the integrals of AO_i * d/dz AO_j -`ao_deriv_1_z `_ +`ao_deriv_1_z `_ array of the integrals of AO_i * d/dx AO_j array of the integrals of AO_i * d/dy AO_j array of the integrals of AO_i * d/dz AO_j -`ao_dipole_x `_ +`ao_dipole_x `_ array of the integrals of AO_i * x AO_j array of the integrals of AO_i * y AO_j array of the integrals of AO_i * z AO_j -`ao_dipole_y `_ +`ao_dipole_y `_ array of the integrals of AO_i * x AO_j array of the integrals of AO_i * y AO_j array of the integrals of AO_i * z AO_j -`ao_dipole_z `_ +`ao_dipole_z `_ array of the integrals of AO_i * x AO_j array of the integrals of AO_i * y AO_j array of the integrals of AO_i * z AO_j -`ao_spread_x `_ +`ao_spread_x `_ array of the integrals of AO_i * x^2 AO_j array of the integrals of AO_i * y^2 AO_j array of the integrals of AO_i * z^2 AO_j -`ao_spread_y `_ +`ao_spread_y `_ array of the integrals of AO_i * x^2 AO_j array of the integrals of AO_i * y^2 AO_j array of the integrals of AO_i * z^2 AO_j -`ao_spread_z `_ +`ao_spread_z `_ array of the integrals of AO_i * x^2 AO_j array of the integrals of AO_i * y^2 AO_j array of the integrals of AO_i * z^2 AO_j -`overlap_bourrin_deriv_x `_ +`overlap_bourrin_deriv_x `_ Undocumented -`overlap_bourrin_dipole `_ +`overlap_bourrin_dipole `_ Undocumented -`overlap_bourrin_spread `_ +`overlap_bourrin_spread `_ Undocumented -`overlap_bourrin_x `_ +`overlap_bourrin_x `_ Undocumented -`overlap_bourrin_x_abs `_ +`overlap_bourrin_x_abs `_ Undocumented -`power `_ +`power `_ Undocumented -`mo_deriv_1_x `_ +`mo_deriv_1_x `_ array of the integrals of MO_i * d/dx MO_j array of the integrals of MO_i * d/dy MO_j array of the integrals of MO_i * d/dz MO_j -`mo_deriv_1_y `_ +`mo_deriv_1_y `_ array of the integrals of MO_i * d/dx MO_j array of the integrals of MO_i * d/dy MO_j array of the integrals of MO_i * d/dz MO_j -`mo_deriv_1_z `_ +`mo_deriv_1_z `_ array of the integrals of MO_i * d/dx MO_j array of the integrals of MO_i * d/dy MO_j array of the integrals of MO_i * d/dz MO_j -`mo_dipole_x `_ +`mo_dipole_x `_ array of the integrals of MO_i * x MO_j array of the integrals of MO_i * y MO_j array of the integrals of MO_i * z MO_j -`mo_dipole_y `_ +`mo_dipole_y `_ array of the integrals of MO_i * x MO_j array of the integrals of MO_i * y MO_j array of the integrals of MO_i * z MO_j -`mo_dipole_z `_ +`mo_dipole_z `_ array of the integrals of MO_i * x MO_j array of the integrals of MO_i * y MO_j array of the integrals of MO_i * z MO_j -`mo_spread_x `_ +`mo_spread_x `_ array of the integrals of MO_i * x^2 MO_j array of the integrals of MO_i * y^2 MO_j array of the integrals of MO_i * z^2 MO_j -`mo_spread_y `_ +`mo_spread_y `_ array of the integrals of MO_i * x^2 MO_j array of the integrals of MO_i * y^2 MO_j array of the integrals of MO_i * z^2 MO_j -`mo_spread_z `_ +`mo_spread_z `_ array of the integrals of MO_i * x^2 MO_j array of the integrals of MO_i * y^2 MO_j array of the integrals of MO_i * z^2 MO_j diff --git a/src/MonoInts/ao_mono_ints.irp.f b/src/Integrals_Monoelec/ao_mono_ints.irp.f similarity index 88% rename from src/MonoInts/ao_mono_ints.irp.f rename to src/Integrals_Monoelec/ao_mono_ints.irp.f index 84e46f61..a59ed3df 100644 --- a/src/MonoInts/ao_mono_ints.irp.f +++ b/src/Integrals_Monoelec/ao_mono_ints.irp.f @@ -7,7 +7,7 @@ BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)] END_DOC do j = 1, ao_num do i = 1, ao_num - ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j) + ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j) + ao_pseudo_integral(i,j) enddo enddo END_PROVIDER diff --git a/src/MonoInts/check_orthonormality.irp.f b/src/Integrals_Monoelec/check_orthonormality.irp.f similarity index 100% rename from src/MonoInts/check_orthonormality.irp.f rename to src/Integrals_Monoelec/check_orthonormality.irp.f diff --git a/src/MonoInts/kin_ao_ints.irp.f b/src/Integrals_Monoelec/kin_ao_ints.irp.f similarity index 100% rename from src/MonoInts/kin_ao_ints.irp.f rename to src/Integrals_Monoelec/kin_ao_ints.irp.f diff --git a/src/MonoInts/kin_mo_ints.irp.f b/src/Integrals_Monoelec/kin_mo_ints.irp.f similarity index 100% rename from src/MonoInts/kin_mo_ints.irp.f rename to src/Integrals_Monoelec/kin_mo_ints.irp.f diff --git a/src/MonoInts/mo_mono_ints.irp.f b/src/Integrals_Monoelec/mo_mono_ints.irp.f similarity index 89% rename from src/MonoInts/mo_mono_ints.irp.f rename to src/Integrals_Monoelec/mo_mono_ints.irp.f index 8d8fed6a..714222ec 100644 --- a/src/MonoInts/mo_mono_ints.irp.f +++ b/src/Integrals_Monoelec/mo_mono_ints.irp.f @@ -7,7 +7,7 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to END_DOC do j = 1, mo_tot_num do i = 1, mo_tot_num - mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) + mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j) enddo enddo END_PROVIDER diff --git a/src/MonoInts/orthonormalize.irp.f b/src/Integrals_Monoelec/orthonormalize.irp.f similarity index 100% rename from src/MonoInts/orthonormalize.irp.f rename to src/Integrals_Monoelec/orthonormalize.irp.f diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/Integrals_Monoelec/pot_ao_ints.irp.f similarity index 99% rename from src/MonoInts/pot_ao_ints.irp.f rename to src/Integrals_Monoelec/pot_ao_ints.irp.f index a9e72e57..1e84d3d4 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_ints.irp.f @@ -10,11 +10,7 @@ integer :: i,j,k,l,n_pt_in,m double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - if (do_pseudo) then - ao_nucl_elec_integral = ao_nucl_elec_integral_pseudo - else - ao_nucl_elec_integral = 0.d0 - endif + ao_nucl_elec_integral = 0.d0 ! _ ! /| / |_) diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f new file mode 100644 index 00000000..95023177 --- /dev/null +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -0,0 +1,226 @@ +BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)] + implicit none + BEGIN_DOC +! Pseudo-potential + 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 +END_PROVIDER + + BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)] + 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 + + ao_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 (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & + !$OMP ao_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, ao_num + + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_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 + ao_pseudo_integral_local(i,j) = ao_pseudo_integral_local(i,j) + & + ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_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(ao_num), '% 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, ao_pseudo_integral_non_local, (ao_num_align,ao_num)] + 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 + + ao_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 (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & + !$OMP ao_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, ao_num + + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_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 + ao_pseudo_integral_non_local(i,j) = ao_pseudo_integral_non_local(i,j) + & + ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_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(ao_num), '% 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 + + + diff --git a/src/MonoInts/pot_mo_ints.irp.f b/src/Integrals_Monoelec/pot_mo_ints.irp.f similarity index 100% rename from src/MonoInts/pot_mo_ints.irp.f rename to src/Integrals_Monoelec/pot_mo_ints.irp.f diff --git a/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f new file mode 100644 index 00000000..6c412e4b --- /dev/null +++ b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f @@ -0,0 +1,33 @@ +BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_num)] + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + BEGIN_DOC + ! interaction nuclear electron on the MO basis + END_DOC + + mo_pseudo_integral = 0.d0 + + if (.not.do_pseudo) then + return + endif + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_pseudo_integral, ao_pseudo_integral) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_pseudo_integral(j,i) = mo_pseudo_integral(j,i) + & + c_j1 * ao_pseudo_integral(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + + diff --git a/src/Pseudo_integrals/int.f90 b/src/Integrals_Monoelec/pseudopot.f90 similarity index 91% rename from src/Pseudo_integrals/int.f90 rename to src/Integrals_Monoelec/pseudopot.f90 index 4eb9d2ae..b605a45e 100644 --- a/src/Pseudo_integrals/int.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 @@ -1003,11 +1003,11 @@ end double precision function crochet(n,g) implicit none integer n -double precision g,pi,dble_fact,expo -pi=dacos(-1.d0) +double precision g,dble_fact,expo +double precision, parameter :: sq_pi_ov_2=dsqrt(dacos(-1.d0)*0.5d0) expo=0.5d0*dfloat(n+1) crochet=dble_fact(n-1)/(2.d0*g)**expo -if(mod(n,2).eq.0)crochet=crochet*dsqrt(pi/2.d0) +if(mod(n,2).eq.0)crochet=crochet*sq_pi_ov_2 end !! @@ -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 @@ -2043,10 +2036,10 @@ double precision function int_prod_bessel_loc(l,gam,n,a) int=0 ! Int f_0 - coef_nk=1.d0/dble_fact(2*n+1) + coef_nk=1.d0/dble_fact( n+n+1 ) expo=0.5d0*dfloat(n+l+1) crochet=dble_fact(n+l-1)/(2.d0*gam)**expo - if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(pi/2.d0) + if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(0.5d0*pi) f_0 = coef_nk*a**n*crochet @@ -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/MonoInts/save_ortho_mos.irp.f b/src/Integrals_Monoelec/save_ortho_mos.irp.f similarity index 100% rename from src/MonoInts/save_ortho_mos.irp.f rename to src/Integrals_Monoelec/save_ortho_mos.irp.f diff --git a/src/MonoInts/spread_dipole_ao.irp.f b/src/Integrals_Monoelec/spread_dipole_ao.irp.f similarity index 100% rename from src/MonoInts/spread_dipole_ao.irp.f rename to src/Integrals_Monoelec/spread_dipole_ao.irp.f diff --git a/src/MonoInts/spread_dipole_mo.irp.f b/src/Integrals_Monoelec/spread_dipole_mo.irp.f similarity index 100% rename from src/MonoInts/spread_dipole_mo.irp.f rename to src/Integrals_Monoelec/spread_dipole_mo.irp.f diff --git a/src/Integrals_Monoelec/tree_dependancy.png b/src/Integrals_Monoelec/tree_dependancy.png new file mode 100644 index 00000000..b12bf22c Binary files /dev/null and b/src/Integrals_Monoelec/tree_dependancy.png differ diff --git a/src/MOGuess/NEEDED_CHILDREN_MODULES b/src/MOGuess/NEEDED_CHILDREN_MODULES index 88d7352e..0a8143ff 100644 --- a/src/MOGuess/NEEDED_CHILDREN_MODULES +++ b/src/MOGuess/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -MonoInts +Integrals_Monoelec diff --git a/src/MOGuess/README.rst b/src/MOGuess/README.rst index 771825ee..2427d23b 100644 --- a/src/MOGuess/README.rst +++ b/src/MOGuess/README.rst @@ -8,7 +8,9 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* `MonoInts `_ +.. image:: tree_dependancy.png + +* `Integrals_Monoelec `_ Documentation ============= diff --git a/src/MOGuess/tree_dependancy.png b/src/MOGuess/tree_dependancy.png new file mode 100644 index 00000000..fa7ab7f1 Binary files /dev/null and b/src/MOGuess/tree_dependancy.png differ diff --git a/src/MOs/README.rst b/src/MOs/README.rst index 5ddd7928..577bf23c 100644 --- a/src/MOs/README.rst +++ b/src/MOs/README.rst @@ -36,6 +36,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `AOs `_ * `Electrons `_ diff --git a/src/MOs/tree_dependancy.png b/src/MOs/tree_dependancy.png new file mode 100644 index 00000000..07d88313 Binary files /dev/null and b/src/MOs/tree_dependancy.png differ diff --git a/src/MP2/README.rst b/src/MP2/README.rst index f68d5936..73494ab2 100644 --- a/src/MP2/README.rst +++ b/src/MP2/README.rst @@ -19,6 +19,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Perturbation `_ * `Selectors_full `_ * `SingleRefMethod `_ diff --git a/src/MP2/tree_dependancy.png b/src/MP2/tree_dependancy.png new file mode 100644 index 00000000..f4bd74fa Binary files /dev/null and b/src/MP2/tree_dependancy.png differ diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst index 5bb2009c..ece272a1 100644 --- a/src/MRCC/README.rst +++ b/src/MRCC/README.rst @@ -8,6 +8,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Perturbation `_ * `Selectors_full `_ * `Generators_full `_ @@ -51,7 +53,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 `_ @@ -60,7 +62,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/MRCC/tree_dependancy.png b/src/MRCC/tree_dependancy.png new file mode 100644 index 00000000..ccd5da28 Binary files /dev/null and b/src/MRCC/tree_dependancy.png differ diff --git a/src/Makefile.config.example b/src/Makefile.config.example index fe6cee8b..410bd0e2 100644 --- a/src/Makefile.config.example +++ b/src/Makefile.config.example @@ -4,6 +4,7 @@ DEBUG = 0 IRPF90_FLAGS+= --align=32 FC = ifort -g +#FC = cache_compile.py ifort -g # Accelerates compilation FCFLAGS= FCFLAGS+= -xHost #FCFLAGS+= -xAVX diff --git a/src/Molden/NEEDED_CHILDREN_MODULES b/src/Molden/NEEDED_CHILDREN_MODULES index b936db90..64795a92 100644 --- a/src/Molden/NEEDED_CHILDREN_MODULES +++ b/src/Molden/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -MOs +MOs Utils diff --git a/src/Molden/README.rst b/src/Molden/README.rst index 128a020a..ec738c9a 100644 --- a/src/Molden/README.rst +++ b/src/Molden/README.rst @@ -31,5 +31,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* `MOs `_ +.. image:: tree_dependancy.png + +* `MOs `_ +* `Utils `_ diff --git a/src/Molden/tree_dependancy.png b/src/Molden/tree_dependancy.png new file mode 100644 index 00000000..23007e0f Binary files /dev/null and b/src/Molden/tree_dependancy.png differ diff --git a/src/MonoInts/EZFIO.cfg b/src/MonoInts/EZFIO.cfg deleted file mode 100644 index 33a8d619..00000000 --- a/src/MonoInts/EZFIO.cfg +++ /dev/null @@ -1,6 +0,0 @@ -[do_pseudo] -type: logical -doc: Using pseudo potential integral of not -interface: input -default: False - diff --git a/src/MonoInts/NEEDED_CHILDREN_MODULES b/src/MonoInts/NEEDED_CHILDREN_MODULES deleted file mode 100644 index be46a359..00000000 --- a/src/MonoInts/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -MOs Pseudo_integrals diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 4533ccfe..6e14bf1d 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Determinants Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD DDCI_selected MRCC Pseudo_integrals +AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Determinants Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils DensityFit QmcChem diff --git a/src/Nuclei/NEEDED_CHILDREN_MODULES b/src/Nuclei/NEEDED_CHILDREN_MODULES index 83260f86..dcdb5f86 100644 --- a/src/Nuclei/NEEDED_CHILDREN_MODULES +++ b/src/Nuclei/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Output +Ezfio_files Utils diff --git a/src/Nuclei/README.rst b/src/Nuclei/README.rst index aaad706d..2e99e8fb 100644 --- a/src/Nuclei/README.rst +++ b/src/Nuclei/README.rst @@ -12,7 +12,10 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -* `Output `_ +.. image:: tree_dependancy.png + +* `Ezfio_files `_ +* `Utils `_ Documentation ============= diff --git a/src/Nuclei/tree_dependancy.png b/src/Nuclei/tree_dependancy.png new file mode 100644 index 00000000..72cfaeee Binary files /dev/null and b/src/Nuclei/tree_dependancy.png differ diff --git a/src/Output/NEEDED_CHILDREN_MODULES b/src/Output/NEEDED_CHILDREN_MODULES deleted file mode 100644 index dcdb5f86..00000000 --- a/src/Output/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -Ezfio_files Utils diff --git a/src/Output/README.rst b/src/Output/README.rst deleted file mode 100644 index 5fe93f50..00000000 --- a/src/Output/README.rst +++ /dev/null @@ -1,62 +0,0 @@ -============= -Output Module -============= - -This module deals with the program I/O in log files. -All output should be printed using routines present in this module. - - - - - - - - - - - - - - - - - - - - - -Needed Modules -============== - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - -* `Ezfio_files `_ -* `Utils `_ - -Documentation -============= - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - -`output_cpu_time_0 `_ - Initial CPU and wall times when printing in the output files - -`output_wall_time_0 `_ - Initial CPU and wall times when printing in the output files - -`write_bool `_ - Write an logical value in output - -`write_double `_ - Write a double precision value in output - -`write_int `_ - Write an integer value in output - -`write_time `_ - Write a time stamp in the output for chronological reconstruction - - - diff --git a/src/Output/output.ezfio_config b/src/Output/output.ezfio_config deleted file mode 100644 index 9e90e610..00000000 --- a/src/Output/output.ezfio_config +++ /dev/null @@ -1,3 +0,0 @@ -output - empty logical - diff --git a/src/Perturbation/EZFIO.cfg b/src/Perturbation/EZFIO.cfg new file mode 100644 index 00000000..28629104 --- /dev/null +++ b/src/Perturbation/EZFIO.cfg @@ -0,0 +1,19 @@ +[do_pt2_end] +type: logical +doc: If true, compute the PT2 at the end of the selection +interface: input +default: True + +[PT2_max] +type: PT2_energy +doc: The selection process stops when the largest PT2 (for all the state) is lower + than pt2_max in absolute value +interface: input +default: 0.0001 + +[var_pt2_ratio] +type: Normalized_float +doc: The selection process stops when the energy ratio variational/(variational+PT2) + is equal to var_pt2_ratio +interface: input +default: 0.75 \ No newline at end of file diff --git a/src/Pseudo_integrals/ASSUMPTIONS.rst b/src/Pseudo/ASSUMPTIONS.rst similarity index 100% rename from src/Pseudo_integrals/ASSUMPTIONS.rst rename to src/Pseudo/ASSUMPTIONS.rst diff --git a/src/Pseudo/EZFIO.cfg b/src/Pseudo/EZFIO.cfg new file mode 100644 index 00000000..bf3909d1 --- /dev/null +++ b/src/Pseudo/EZFIO.cfg @@ -0,0 +1,82 @@ +[pseudo_klocmax] +doc: test +type:integer +interface: input_without_default + +[pseudo_n_k] +doc: test +type: integer +interface: input_without_default +size: (nuclei.nucl_num,pseudo.pseudo_klocmax) + +[pseudo_v_k] +doc: test +type: double precision +interface: input_without_default +size: (nuclei.nucl_num,pseudo.pseudo_klocmax) + +[pseudo_dz_k] +doc: test +type: double precision +interface: input_without_default +size: (nuclei.nucl_num,pseudo.pseudo_klocmax) + +[pseudo_lmax] +doc: test +type:integer +interface: input_without_default + +[pseudo_kmax] +doc: test +type:integer +interface: input_without_default + +[pseudo_n_kl] +doc: test +type: integer +interface: input_without_default +size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax) + +[pseudo_v_kl] +doc: test +type: double precision +interface: input_without_default +size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax) + +[pseudo_dz_kl] +doc: test +type: double precision +interface: input_without_default +size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax) + +[do_pseudo] +type: logical +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: 100 + +[pseudo_grid_rmax] +type: double precision +doc: R_maxof the QMC grid +interface: input +default: 4.0 + +[pseudo_grid] +type: double precision +doc: QMC grid +interface: output +size: (pseudo.pseudo_grid_size,ao_basis.ao_num,-pseudo.pseudo_lmax:pseudo.pseudo_lmax,0:pseudo.pseudo_lmax,nuclei.nucl_num) + +[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_integrals/Makefile b/src/Pseudo/Makefile similarity index 82% rename from src/Pseudo_integrals/Makefile rename to src/Pseudo/Makefile index 5cf11b78..06dc50ff 100644 --- a/src/Pseudo_integrals/Makefile +++ b/src/Pseudo/Makefile @@ -1,6 +1,6 @@ # Define here all new external source files and objects.Don't forget to prefix the # object files with IRPF90_temp/ -SRC=int.f90 -OBJ=IRPF90_temp/int.o +SRC= +OBJ= include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/Pseudo/NEEDED_CHILDREN_MODULES b/src/Pseudo/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..d47a550b --- /dev/null +++ b/src/Pseudo/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Ezfio_files Nuclei diff --git a/src/Pseudo/README.rst b/src/Pseudo/README.rst new file mode 100644 index 00000000..f6e7fc6e --- /dev/null +++ b/src/Pseudo/README.rst @@ -0,0 +1,23 @@ +============= +Pseudo Module +============= + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +.. image:: tree_dependancy.png + +* `Ezfio_files `_ +* `Nuclei `_ + diff --git a/src/Pseudo/tree_dependancy.png b/src/Pseudo/tree_dependancy.png new file mode 100644 index 00000000..98e8be37 Binary files /dev/null and b/src/Pseudo/tree_dependancy.png differ diff --git a/src/Pseudo_integrals/EZFIO.cfg b/src/Pseudo_integrals/EZFIO.cfg deleted file mode 100644 index f56bc325..00000000 --- a/src/Pseudo_integrals/EZFIO.cfg +++ /dev/null @@ -1,5 +0,0 @@ -[do_pseudo] -type: logical -doc: Using pseudo potential integral of not -interface: input -default: False \ No newline at end of file diff --git a/src/Pseudo_integrals/NEEDED_CHILDREN_MODULES b/src/Pseudo_integrals/NEEDED_CHILDREN_MODULES deleted file mode 100644 index 88c6f86b..00000000 --- a/src/Pseudo_integrals/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -AOs Electrons diff --git a/src/Pseudo_integrals/README.rst b/src/Pseudo_integrals/README.rst deleted file mode 100644 index 08076556..00000000 --- a/src/Pseudo_integrals/README.rst +++ /dev/null @@ -1,24 +0,0 @@ -======================= -Pseudo_integrals Module -======================= - -Documentation -============= - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - -`ao_nucl_elec_integral_pseudo `_ - interaction nuclear electron - - - -Needed Modules -============== - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - -* `AOs `_ -* `Electrons `_ - diff --git a/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f b/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f deleted file mode 100644 index 37dc7591..00000000 --- a/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f +++ /dev/null @@ -1,167 +0,0 @@ - BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo, (ao_num_align,ao_num)] - BEGIN_DOC -! interaction nuclear electron - END_DOC - implicit none - 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 - - ao_nucl_elec_integral_pseudo = 0.d0 - - ! - ! | _ _ _. | - ! |_ (_) (_ (_| | - ! - !! Parameters of the local part of pseudo: - - integer klocmax - integer, allocatable :: n_k(:,:) - double precision, allocatable :: v_k(:,:), dz_k(:,:) - - call ezfio_get_pseudo_integrals_klocmax(klocmax) - - allocate(n_k(nucl_num,klocmax),v_k(nucl_num,klocmax), dz_k(nucl_num,klocmax)) - - call ezfio_get_pseudo_integrals_v_k(v_k) - call ezfio_get_pseudo_integrals_n_k(n_k) - call ezfio_get_pseudo_integrals_dz_k(dz_k) - - !! Dump array - integer, allocatable :: n_k_dump(:) - double precision, allocatable :: v_k_dump(:), dz_k_dump(:) - - allocate(n_k_dump(1:klocmax), v_k_dump(1:klocmax), dz_k_dump(1:klocmax)) - - - ! - ! |\ | _ ._ | _ _ _. | - ! | \| (_) | | | (_) (_ (_| | - ! - !! Parameters of non local part of pseudo: - - integer :: kmax,lmax - integer, allocatable :: n_kl(:,:,:) - double precision, allocatable :: v_kl(:,:,:), dz_kl(:,:,:) - - call ezfio_get_pseudo_integrals_lmaxpo(lmax) - call ezfio_get_pseudo_integrals_kmax(kmax) - !lmax plus one -> lmax - lmax = lmax - 1 - - allocate(n_kl(nucl_num,kmax,0:lmax), v_kl(nucl_num,kmax,0:lmax), dz_kl(nucl_num,kmax,0:lmax)) - - call ezfio_get_pseudo_integrals_n_kl(n_kl) - call ezfio_get_pseudo_integrals_v_kl(v_kl) - call ezfio_get_pseudo_integrals_dz_kl(dz_kl) - - - !! Dump array - integer, allocatable :: n_kl_dump(:,:) - double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:) - - allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax)) - - ! _ - ! / _. | _ | - ! \_ (_| | (_ |_| | - ! - - write(output_monoints,*) '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, n_kl_dump, v_kl_dump, dz_kl_dump, & - !$OMP wall_0,wall_2,thread_num, output_monoints) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & - !$OMP ao_nucl_elec_integral_pseudo,nucl_num,nucl_charge, & - !$OMP klocmax,lmax,kmax,v_k,n_k, dz_k, n_kl, v_kl, dz_kl, & - !$OMP wall_1) - - !$OMP DO SCHEDULE (guided) - - do j = 1, ao_num - - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - num_B = ao_nucl(i) - power_B(1:3)= ao_power(i,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) - - do l=1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(l,j) - - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_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 = v_k(k,1:klocmax) - n_k_dump = n_k(k,1:klocmax) - dz_k_dump = dz_k(k,1:klocmax) - - c = c + Vloc(klocmax, v_k_dump,n_k_dump, dz_k_dump, & - A_center,power_A,alpha,B_center,power_B,beta,C_center) - - - n_kl_dump = n_kl(k,1:kmax,0:lmax) - v_kl_dump = v_kl(k,1:kmax,0:lmax) - dz_kl_dump = dz_kl(k,1:kmax,0:lmax) - - c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) - - enddo - ao_nucl_elec_integral_pseudo(i,j) = ao_nucl_elec_integral_pseudo(i,j) + & - ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_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 - write(output_monoints,*) 100.*float(j)/float(ao_num), '% in ', & - wall_2-wall_1, 's' - endif - endif - enddo - - !$OMP END DO - !$OMP END PARALLEL - - -! _ -! | \ _ _. | | _ _ _. _|_ _ -! |_/ (/_ (_| | | (_) (_ (_| |_ (/_ -! - - deallocate(n_k,v_k, dz_k) - deallocate(n_k_dump,v_k_dump, dz_k_dump) - - deallocate(n_kl,v_kl, dz_kl) - deallocate(n_kl_dump,v_kl_dump, dz_kl_dump) - - -END_PROVIDER diff --git a/src/Pseudo_integrals/pseudo.ezfio_config b/src/Pseudo_integrals/pseudo.ezfio_config deleted file mode 100644 index a01bc2c0..00000000 --- a/src/Pseudo_integrals/pseudo.ezfio_config +++ /dev/null @@ -1,10 +0,0 @@ -pseudo_integrals - klocmax integer - v_k double precision (nuclei_nucl_num,pseudo_integrals_klocmax) - n_k integer (nuclei_nucl_num,pseudo_integrals_klocmax) - dz_k double precision (nuclei_nucl_num,pseudo_integrals_klocmax) - lmaxpo integer - kmax integer - v_kl double precision (nuclei_nucl_num,pseudo_integrals_kmax,pseudo_integrals_lmaxpo) - n_kl integer (nuclei_nucl_num,pseudo_integrals_kmax,pseudo_integrals_lmaxpo) - dz_kl double precision (nuclei_nucl_num,pseudo_integrals_kmax,pseudo_integrals_lmaxpo) diff --git a/src/QmcChem/ASSUMPTIONS.rst b/src/QmcChem/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/MonoInts/Makefile b/src/QmcChem/Makefile similarity index 73% rename from src/MonoInts/Makefile rename to src/QmcChem/Makefile index b1f3b02c..06dc50ff 100644 --- a/src/MonoInts/Makefile +++ b/src/QmcChem/Makefile @@ -3,4 +3,4 @@ SRC= OBJ= -include $(QPACKAGE_ROOT)/src/Makefile.common \ No newline at end of file +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..c4b467d7 --- /dev/null +++ b/src/QmcChem/README.rst @@ -0,0 +1,23 @@ +============== +QmcChem Module +============== + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + + +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..24d75504 --- /dev/null +++ b/src/QmcChem/pot_ao_pseudo_ints.irp.f @@ -0,0 +1,275 @@ +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, pseudo_grid, (pseudo_grid_size,ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num) ] + 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 + + pseudo_grid = 0.d0 + 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 j=1,pseudo_grid_size + 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)) + pseudo_grid(j,i,m,l,k) = pseudo_grid(j,i,m,l,k) + & + ao_coef_normalized_ordered_transp(p,i)*y + enddo + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + + diff --git a/src/QmcChem/pseudo.irp.f b/src/QmcChem/pseudo.irp.f new file mode 100644 index 00000000..b0cd57a5 --- /dev/null +++ b/src/QmcChem/pseudo.irp.f @@ -0,0 +1,170 @@ +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_pseudo_grid(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..cfc18338 --- /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 save_dets_qmcchem + call write_spindeterminants +! call write_pseudopotential +! call save_aux_basis +end diff --git a/src/Selectors_full/README.rst b/src/Selectors_full/README.rst index 11286fce..0ffabb6c 100644 --- a/src/Selectors_full/README.rst +++ b/src/Selectors_full/README.rst @@ -165,6 +165,8 @@ Needed Modules .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +.. image:: tree_dependancy.png + * `Determinants `_ * `Hartree_Fock `_ diff --git a/src/Selectors_full/tree_dependancy.png b/src/Selectors_full/tree_dependancy.png new file mode 100644 index 00000000..08ba5eda Binary files /dev/null and b/src/Selectors_full/tree_dependancy.png differ diff --git a/src/Utils/README.rst b/src/Utils/README.rst index ff16f3ed..e36ab45e 100644 --- a/src/Utils/README.rst +++ b/src/Utils/README.rst @@ -242,7 +242,7 @@ Documentation `align_double `_ Compute 1st dimension such that it is aligned for vectorization. -`approx_dble `_ +`approx_dble `_ Undocumented `binom `_ @@ -257,10 +257,16 @@ Documentation `binom_transp `_ Binomial coefficients -`dble_fact `_ +`dble_fact `_ + Undocumented + +`dble_fact_even `_ n!! -`dble_logfact `_ +`dble_fact_odd `_ + n!! + +`dble_logfact `_ n!! `fact `_ @@ -269,29 +275,29 @@ Documentation `fact_inv `_ 1/n! -`inv_int `_ +`inv_int `_ 1/i `logfact `_ n! -`normalize `_ +`normalize `_ Normalizes vector u u is expected to be aligned in memory. -`nproc `_ +`nproc `_ Number of current OpenMP threads -`u_dot_u `_ +`u_dot_u `_ Compute -`u_dot_v `_ +`u_dot_v `_ Compute -`wall_time `_ +`wall_time `_ The equivalent of cpu_time, but for the wall time. -`write_git_log `_ +`write_git_log `_ Write the last git commit in file iunit. diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index ec1e4200..d356454b 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -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/src/Utils/tree_dependancy.png b/src/Utils/tree_dependancy.png new file mode 100644 index 00000000..38b04785 Binary files /dev/null and b/src/Utils/tree_dependancy.png differ diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index 7a78ad59..705f4f7b 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -138,21 +138,21 @@ END_PROVIDER double precision function dble_fact(n) implicit none integer :: n - double precision :: dble_fact_peer, dble_fact_odd + double precision :: dble_fact_even, dble_fact_odd dble_fact = 1.d0 if(n.lt.0) return if(iand(n,1).eq.0)then - dble_fact = dble_fact_peer(n) + dble_fact = dble_fact_even(n) else dble_fact= dble_fact_odd(n) endif end function -double precision function dble_fact_peer(n) result(fact2) +double precision function dble_fact_even(n) result(fact2) implicit none BEGIN_DOC ! n!! diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py index 23e831ec..598300c5 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 = 2.e-7 # A test get a geo file and a basis file. # A global dict containt the result for this test @@ -88,12 +88,12 @@ def check_disk_acess(geo, basis, mult=1): # ~#~#~#~#~#~#~#~#~#~#~#~#~ # # Test 1 - ezfio.bielec_integrals_disk_access_ao_integrals = "Write" + ezfio.integrals_bielec_disk_access_ao_integrals = "Write" cmd = "qp_edit -c {0}.ezfio".format(filename) subprocess.check_call([cmd], shell=True) # Test 2 - ezfio.bielec_integrals_disk_access_ao_integrals = "IculeAcess" + ezfio.integrals_bielec_disk_access_ao_integrals = "IculeAcess" cmd = "qp_edit -c {0}.ezfio".format(filename) try: @@ -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.11982530413317) # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # # G l o b a l _ v a r i a b l e # @@ -199,7 +199,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True): ezfio.hartree_fock_thresh_scf = 1.e-10 ezfio.hartree_fock_n_it_scf_max = 100 - ezfio.pseudo_integrals_do_pseudo = pseudo + ezfio.pseudo_do_pseudo = pseudo # ~#~#~ # # R u n # @@ -244,8 +244,8 @@ def run_full_ci_10k_pt2_end(geo, basis, pseudo): ref_energy_var = defaultdict(dict) ref_energy_pt2 = defaultdict(dict) - ref_energy_var["sto-3g"]["methane"] = Energy(-0.398058753535695E+02, None) - ref_energy_pt2["sto-3g"]["methane"] = Energy(-0.398059182483741E+02, None) + ref_energy_var["sto-3g"]["methane"] = Energy(-39.8058687211, None) + ref_energy_pt2["sto-3g"]["methane"] = Energy(-39.8059180427, None) # ~#~#~#~ # # I n i t #