diff --git a/data/ezfio_defaults/WILL_BE_DELETED.ezfio_default b/data/ezfio_defaults/WILL_BE_DELETED.ezfio_default new file mode 100644 index 00000000..79f0f114 --- /dev/null +++ b/data/ezfio_defaults/WILL_BE_DELETED.ezfio_default @@ -0,0 +1,66 @@ +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 + +cassd + n_det_max_cassd 10000 + pt2_max 1.e-4 + do_pt2_end True + +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 \ No newline at end of file diff --git a/data/ezfio_defaults/q_package.ezfio_default b/data/ezfio_defaults/q_package.ezfio_default index 5da4d813..347e3751 100644 --- a/data/ezfio_defaults/q_package.ezfio_default +++ b/data/ezfio_defaults/q_package.ezfio_default @@ -32,6 +32,12 @@ full_ci 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 diff --git a/ocaml/test_input.ml b/ocaml/test_input.ml deleted file mode 100644 index 956fd745..00000000 --- a/ocaml/test_input.ml +++ /dev/null @@ -1,177 +0,0 @@ -open Qptypes;; - -let test_ao () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Ao_basis.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Ao_basis.to_string b); - print_endline (Input.Ao_basis.to_rst b |> Rst_string.to_string); -;; - -let test_bielec_intergals () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Bielec_integrals.read () with - | Some x -> x - | None -> assert false - in - let output = Input.Bielec_integrals.to_string b - in - print_endline output; - let rst = Input.Bielec_integrals.to_rst b in - let b2 = match Input.Bielec_integrals.of_rst rst with - | Some x -> x - | None -> assert false - in - if (b = b2) then - print_endline "OK" - else - print_endline "rst failed"; -;; - -let test_bitmasks () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Bitmasks.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Bitmasks.to_string b); -;; - - -let test_dets () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Determinants.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Determinants.to_rst b |> Rst_string.to_string ) ; - print_endline (Input.Determinants.sexp_of_t b |> Sexplib.Sexp.to_string ) ; - let rst = Input.Determinants.to_rst b in - let b2 = match Input.Determinants.of_rst rst with - | Some x -> x - | None -> assert false - in - if (b2 = b) then - print_endline "OK" - else - print_endline "Failed" -;; - -let test_cisd_sc2 () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Cisd_sc2_selected.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Cisd_sc2_selected.to_string b); - let rst = Input.Cisd_sc2_selected.to_rst b in - let b2 = match Input.Cisd_sc2_selected.of_rst rst with - | Some x -> x - | None -> assert false - in - if (b = b2) then - print_endline "OK" - else - print_endline "rst failed"; - -;; - -let test_electrons () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Electrons.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Electrons.to_string b); - let rst = Input.Electrons.to_rst b in - let b2 = match Input.Electrons.of_rst rst with - | Some x -> x - | None -> assert false - in - if (b = b2) then - print_endline "OK" - else - print_endline "Failed in rst" -;; - -let test_fci () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Full_ci.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Full_ci.to_string b); - let rst = Input.Full_ci.to_rst b in - let b2 = match Input.Full_ci.of_rst rst with - | Some x -> x - | None -> assert false - in - print_endline (Input.Full_ci.to_string b); - if (b = b2) then - print_endline "OK" - else - print_endline "Failed in rst" -;; - -let test_hf () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Hartree_fock.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Hartree_fock.to_string b); - let rst = Input.Hartree_fock.to_rst b in - let b2 = match Input.Hartree_fock.of_rst rst with - | Some x -> x - | None -> assert false - in - print_endline (Input.Hartree_fock.to_string b); - if (b = b2) then - print_endline "OK" - else - print_endline "Failed in rst" -;; - -let test_mo () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Mo_basis.read () with - | Some x -> x - | None -> assert false - in - print_endline (Input.Mo_basis.to_string b); -;; - -let test_nucl () = - Ezfio.set_file "F2.ezfio" ; - let b = match Input.Nuclei.read () with - | Some x -> x - | None -> assert false - in - let rst = Input.Nuclei.to_rst b in - let b2 = match Input.Nuclei.of_rst rst with - | Some x -> x - | None -> assert false - in - print_endline (Input.Nuclei.to_string b); - if (b = b2) then - print_endline "OK" - else - print_endline "Failed in rst" -;; - -(* -test_ao ();; -test_bitmasks (); -test_cis (); -test_cisd_sc2 (); -test_dets (); -test_hf ();; -test_mo ();; -test_nucl (); -test_bielec_intergals ();; -test_electrons(); -*) -test_dets (); - diff --git a/scripts/ezfio_interface.py b/scripts/ezfio_interface.py new file mode 100755 index 00000000..7c196cf8 --- /dev/null +++ b/scripts/ezfio_interface.py @@ -0,0 +1,333 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Take a path in argv +Check if EZFIO.cfg exists. +EZFIO.cfg are in MODULE directories. +create : ezfio_interface.irp.f + folder_ezfio_inteface_config + +Example EZFIO.cfg: +``` +[thresh_SCF] +doc: Threshold on the convergence of the Hartree Fock energy +type: Threshold +default: 1.e-10 + +[do_pt2_end] +type: logical +doc: If true, compute the PT2 at the end of the selection +default: true +``` + +""" + +import sys +import os +import os.path + +import ConfigParser + +from collections import defaultdict +from collections import namedtuple + +Type = namedtuple('Type', 'ocaml fortran') +def bool_convertor(b): + return ( b.lower() in [ "true", ".true." ] ) + + +def get_type_dict(): + """ + This function makes the correspondance between the type of value read in + ezfio.cfg into the f90 and Ocam Type + return fancy_type[fancy_type] = namedtuple('Type', 'ocaml fortran') + For example fancy_type['Ndet'].fortran = interger + .ocaml = int + """ + + # ~#~#~#~ # + # I n i t # + # ~#~#~#~ # + + # Dict to change ocaml LowLevel type into FortranLowLevel type + ocaml_to_fortran = {"int": "integer", + "float": "double precision", + "logical": "logical", + "string": "character*60"} + + fancy_type = defaultdict(dict) + + # ~#~#~#~#~#~#~#~ # + # R a w _ t y p e # + # ~#~#~#~#~#~#~#~ # + + fancy_type['integer'] = Type("int", "integer") + fancy_type['int'] = Type("int", "integer") + + fancy_type['float'] = Type("float", "double precision") + fancy_type['double precision'] = Type("float", "double precision") + + fancy_type['logical'] = Type("bool", "logical") + fancy_type['bool'] = Type("bool", "logical") + + # ~#~#~#~#~#~#~#~ # + # q p _ t y p e s # + # ~#~#~#~#~#~#~#~ # + + src = os.environ['QPACKAGE_ROOT'] + "/ocaml/qptypes_generator.ml" + + with open(src, "r") as f: + l = [i for i in f.read().splitlines() if i.strip().startswith("*")] + + for i in l: + ocaml_fancy_type = i.split()[1].strip() + ocaml_type = i.split()[3] + fortran_type = ocaml_to_fortran[ocaml_type] + + fancy_type[ocaml_fancy_type] = Type(ocaml_type, fortran_type) + + return dict(fancy_type) + + +type_dict = get_type_dict() + + +def get_dict_config_file(config_file_path, module_lower): + """ + Input: + config_file_path is the config file path + (for example FULL_PATH/EZFIO.cfg) + module_lower is the MODULE name lowered + (Ex fullci) + + Return a dict d[provider_name] = {type, + doc, + ezfio_name, + ezfio_dir, + interface, + default} + + Type : Is a fancy_type named typle who containt fortran and ocaml type + doc : Is the doc + ezfio_name : Will be the name of the file + ezfio_dir : Will be the folder who containt the ezfio_name + * /ezfio_dir/ezfio_name + * equal to MODULE_lower name for the moment. + interface : The provider is a imput or a output + if is a output: + default : The default value + + """ + # ~#~#~#~ # + # I n i t # + # ~#~#~#~ # + + d = defaultdict(dict) + l_info_required = ["doc", "interface"] + l_info_optional = ["ezfio_name"] + + # ~#~#~#~#~#~#~#~#~#~#~ # + # L o a d _ C o n f i g # + # ~#~#~#~#~#~#~#~#~#~#~ # + + config_file = ConfigParser.ConfigParser() + config_file.readfp(open(config_file_path)) + + # ~#~#~#~#~#~#~#~#~ # + # F i l l _ d i c t # + # ~#~#~#~#~#~#~#~#~ # + + def error(o, p, c): + "o option ; p provider_name ;c config_file_path" + print "You need a {0} for {1} in {2}".format(o, p, c) + + for section in config_file.sections(): + # pvd = provider + pvd = section.lower() + + # Create the dictionary who containt the value per default + d_default = {"ezfio_name": pvd} + + # Set the ezfio_dir + d[pvd]["ezfio_dir"] = module_lower + + # Check if type if avalaible + type_ = config_file.get(section, "type") + if type_ not in type_dict: + print "{0} not avalaible. Choose in:".format(type_) + print ", ".join([i for i in type_dict]) + sys.exit(1) + else: + d[pvd]["type"] = type_dict[type_] + + # Fill the dict with REQUIRED information + for option in l_info_required: + try: + d[pvd][option] = config_file.get(section, option) + except ConfigParser.NoOptionError: + error(option, pvd, config_file_path) + sys.exit(1) + + # Fill the dict with OPTIONAL information + for option in l_info_optional: + try: + d[pvd][option] = config_file.get(section, option).lower() + except ConfigParser.NoOptionError: + d[pvd][option] = d_default[option] + + # If interface is output we need a default value information + if d[pvd]["interface"] == "input": + try: + d[pvd]["default"] = config_file.get(section, "default") + except ConfigParser.NoOptionError: + error("default", pvd, config_file_path) + sys.exit(1) + + return dict(d) + + +def create_ezfio_provider(dict_ezfio_cfg): + """ + From dict d[provider_name] = {type, + doc, + ezfio_name, + ezfio_dir, + interface, + default} + create the a list who containt all the code for the provider + return [code, ...] + """ + from ezfio_with_default 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: + ez_p.set_type(dict_info['type'].fortran) + ez_p.set_name(provider_name) + ez_p.set_doc(dict_info['doc']) + ez_p.set_ezfio_dir(dict_info['ezfio_dir']) + ez_p.set_ezfio_name(dict_info['ezfio_name']) + ez_p.set_default(dict_info['default']) + ez_p.set_output("output_%s" % dict_info['ezfio_dir']) + + dict_code_provider[provider_name] = str(ez_p) + + return dict_code_provider + + +def save_ezfio_provider(path_head, dict_code_provider): + """ + Write in path_head/"ezfio_interface.irp.f" the value of dict_code_provider + """ + + path = "{0}/ezfio_interface.irp.f".format(path_head) + +# print "Path = {}".format(path) + + try: + f = open(path, "r") + except IOError: + old_output = "" + else: + old_output = f.read() + f.close() + + output = "! DO NOT MODIFY BY HAND\n" + \ + "! Created by $QPACKAGE_ROOT/scripts/ezfio_interface.py\n" + \ + "! from file {0}/EZFIO.cfg\n".format(path_head) + \ + "\n" + for provider_name, code in dict_code_provider.iteritems(): + output += str(code) + "\n" + + if output != old_output: + with open(path, "w") as f: + f.write(output) + + +def create_ezfio_config(dict_ezfio_cfg, opt, module_lower): + """ + From dict_ezfio_cfg[provider_name] = {type, default, ezfio_name,ezfio_dir,doc} + Return the string ezfio_interface_config + """ + + result = [module_lower] + lenmax = max([len(i) for i in dict_ezfio_cfg]) + 2 + l = sorted(dict_ezfio_cfg.keys()) + for provider_name in l: + provider_info = dict_ezfio_cfg[provider_name] + s = " {0} {1}".format( + provider_name.lower().ljust(lenmax), + provider_info["type"].fortran) + result.append(s) + return "\n".join(result) + + +def save_ezfio_config(module_lower, str_ezfio_config): + """ + Write the str_ezfio_config in + $QPACKAGE_ROOT/EZFIO/{0}.ezfio_interface_config".format(module_lower) + """ + + ezfio_dir = "{0}/EZFIO".format(os.environ['QPACKAGE_ROOT']) + path = "{0}/config/{1}.ezfio_interface_config".format(ezfio_dir, + module_lower) + +# print "Path = {}".format(path) + + try: + f = open(path, "r") + except IOError: + old_output = "" + else: + old_output = f.read() + f.close() + + if str_ezfio_config != old_output: + with open(path, "w") as f: + f.write(str_ezfio_config) + + +def main(): + """ + Two condition: + -Take the EZFIO.cfg path in arg + or + -Look if EZFIO.cfg is present in the pwd + """ + + try: + config_file_path = sys.argv[1] + except: + config_file_path = "EZFIO.cfg" + if "EZFIO.cfg" not in os.listdir(os.getcwd()): + sys.exit(0) + + config_file_path = os.path.expanduser(config_file_path) + config_file_path = os.path.expandvars(config_file_path) + config_file_path = os.path.abspath(config_file_path) +# print config_file_path + + path_dirname = os.path.dirname(config_file_path) + module = [i for i in path_dirname.split("/") if i][-1] + module_lower = module.lower() + +# print "Read {0}".format(config_file_path) + dict_info_provider = get_dict_config_file(config_file_path, module_lower) + +# print "Generating the ezfio_interface.irp.f: \n" + d_config = create_ezfio_provider(dict_info_provider) + +# print "Saving the ezfio_interface.irp.f" + save_ezfio_provider(path_dirname, d_config) + +# print "Generating the ezfio_config" + config_ezfio = create_ezfio_config(dict_info_provider, "config", module_lower) + +# print "Saving ezfio_config" + save_ezfio_config(module_lower, config_ezfio) + + +if __name__ == "__main__": + main() diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py index 8fad8511..cc9a04ea 100755 --- a/scripts/ezfio_interface/ei_handler.py +++ b/scripts/ezfio_interface/ei_handler.py @@ -507,6 +507,9 @@ def create_ocaml_input(dict_ezfio_cfg, module_lower): l_type.append(v["type"]) l_doc.append(v["doc"]) + if not l_ezfio_name: + raise ValueError + e_glob = EZFIO_ocaml(l_ezfio_name=l_ezfio_name, l_type=l_type, l_doc=l_doc) @@ -602,6 +605,49 @@ def save_ocaml_input(module_lower, str_ocaml_input): f.write(str_ocaml_input) +def get_l_module_lower(): + """ + Get all module who have EZFIO.cfg with input data + (NB `search` in all the ligne and `match` only in one) + """ + + # ~#~#~#~ # + # I n i t # + # ~#~#~#~ # + + mypath = "{0}/src".format(os.environ['QPACKAGE_ROOT']) + + # ~#~#~#~#~#~#~#~ # + # L _ f o l d e r # + # ~#~#~#~#~#~#~#~ # + + from os import listdir + from os.path import isdir, join, exists + + l_folder = [f for f in listdir(mypath) if isdir(join(mypath, f))] + + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + # L _ m o d u l e _ l o w e r # + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + + l_module_lower = [] + import re + p = re.compile(ur'interface:\s+input') + + for f in l_folder: + path = "{0}/{1}/EZFIO.cfg".format(mypath, f) + if exists(path): + with open(path, 'r') as file_: + if p.search(file_.read()): + l_module_lower.append(f.lower()) + + # ~#~#~#~#~#~ # + # R e t u r n # + # ~#~#~#~#~#~ # + + return l_module_lower + + def create_ocaml_input_global(): """ Check for all the EZFIO.cfg get the module lower @@ -612,15 +658,7 @@ def create_ocaml_input_global(): # I n i t # # ~#~#~#~# # - from os import listdir - from os.path import isdir, join, exists - - mypath = "{0}/src".format(os.environ['QPACKAGE_ROOT']) - - onlyfodler = [f for f in listdir(mypath) if isdir(join(mypath, f))] - - l_module_lower = [f.lower() for f in onlyfodler - if exists("{0}/{1}/EZFIO.cfg".format(mypath, f))] + l_module_lower = get_l_module_lower() # ~#~#~#~#~#~#~#~# # # C r e a t i o n # @@ -761,8 +799,12 @@ if __name__ == "__main__": # O c a m l # # ~#~#~#~#~#~# if do_all or arguments["--ocaml"]: - str_ocaml_input = create_ocaml_input(dict_ezfio_cfg, module_lower) - save_ocaml_input(module_lower, str_ocaml_input) + try: + str_ocaml_input = create_ocaml_input(dict_ezfio_cfg, module_lower) + except ValueError: + pass + else: + save_ocaml_input(module_lower, str_ocaml_input) # ~#~#~#~#~#~#~#~#~#~#~#~#~ # # e z f i o _ d e f a u l t # diff --git a/scripts/ezfio_with_default.py b/scripts/ezfio_with_default.py index 46981132..907ca45e 100755 --- a/scripts/ezfio_with_default.py +++ b/scripts/ezfio_with_default.py @@ -90,19 +90,14 @@ END_PROVIDER self.default = t def get_default(self): + filename = '/'.join( [os.environ['QPACKAGE_ROOT'], 'data', + 'ezfio_defaults', + 'WILL_BE_DELETED.ezfio_default'] ) - from os import listdir - from os.path import isfile, join - - mypath = '/'.join( [os.environ['QPACKAGE_ROOT'], 'data', 'ezfio_defaults','/']) - onlyfiles = [ f for f in listdir(mypath) if isfile(join(mypath,f)) ] - - lines= [] - for filename in onlyfiles: - file = open(mypath+filename,'r') - lines.extend(file.readlines()[:]) - file.close() - + file = open(filename,'r') + lines = file.readlines() + file.close() + k=-1 # Search directory for k,line in enumerate(lines): if line[0] != ' ': @@ -120,16 +115,19 @@ END_PROVIDER break v = buffer[1] name = self.name + true = True + false= False try: v_eval = eval(v) + except: + v = "call ezfio_get_%(v)s(%(name)s)"%locals() + else: if type(v_eval) == bool: v = '.%s.'%(v) elif type(v_eval) == float: v = v.replace('e','d') v = v.replace('E','D') v = "%(name)s = %(v)s"%locals() - except: - v = "call ezfio_get_%(v)s(%(name)s)"%locals() self.default = v diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index ad6a57ee..280c9f72 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -53,7 +53,7 @@ class H_apply(object): !$OMP N_elec_in_key_hole_2,ia_ja_pairs) & !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, & !$OMP hole_1, particl_1, hole_2, particl_2, & - !$OMP elec_alpha_num,i_generator)""" + !$OMP elec_alpha_num,i_generator) FIRSTPRIVATE(iproc)""" s["omp_end_parallel"] = "!$OMP END PARALLEL" s["omp_master"] = "!$OMP MASTER" s["omp_end_master"] = "!$OMP END MASTER" diff --git a/scripts/qp_convert.py b/scripts/qp_convert.py index 54a6d2da..ab008e9e 100755 --- a/scripts/qp_convert.py +++ b/scripts/qp_convert.py @@ -53,16 +53,68 @@ def write_ezfioFile(res,filename): basis = res.uncontracted_basis geom = res.geometry + res.clean_contractions() + # AO Basis + import string + at = [] + num_prim = [] + magnetic_number = [] + angular_number = [] + power_x = [] + power_y = [] + power_z = [] + coefficient = [] + exponent = [] + res.convert_to_cartesian() + for b in res.basis: + c = b.center + for i,atom in enumerate(res.geometry): + if atom.coord == c: + at.append(i+1) + num_prim.append(len(b.prim)) + s = b.sym + power_x.append( string.count(s,"x") ) + power_y.append( string.count(s,"y") ) + power_z.append( string.count(s,"z") ) + coefficient.append( b.coef ) + exponent.append( [ p.expo for p in b.prim ] ) + ezfio.set_ao_basis_ao_num(len(res.basis)) + ezfio.set_ao_basis_ao_nucl(at) + ezfio.set_ao_basis_ao_prim_num(num_prim) + ezfio.set_ao_basis_ao_power(power_x+power_y+power_z) + prim_num_max = ezfio.get_ao_basis_ao_prim_num_max() + len_res_basis = len(res.basis) + for i in range(len(res.basis)): + coefficient[i] += [ 0. for j in range(len(coefficient[i]),prim_num_max) ] + exponent[i] += [ 0. for j in range(len(exponent[i]),prim_num_max) ] + coefficient = reduce(lambda x, y: x+y, coefficient, []) + exponent = reduce(lambda x, y: x+y, exponent, []) + coef = [] + expo = [] + for i in range(prim_num_max): + for j in range(i,len(coefficient),prim_num_max): + coef.append ( coefficient[j] ) + expo.append ( exponent[j] ) + ezfio.set_ao_basis_ao_coef(coef) + ezfio.set_ao_basis_ao_expo(expo) + ezfio.set_ao_basis_ao_basis("Read by resultsFile") + + # MO MoTag = res.determinants_mo_type ezfio.set_mo_basis_mo_label('Orthonormalized') MO_type = MoTag - allMOs = res.uncontracted_mo_sets[MO_type] + allMOs = res.mo_sets[MO_type] - closed = [ (allMOs[i].eigenvalue,i) for i in res.closed_mos ] - active = [ (allMOs[i].eigenvalue,i) for i in res.active_mos ] - virtual =[ (allMOs[i].eigenvalue,i) for i in res.virtual_mos ] + try: + closed = [ (allMOs[i].eigenvalue,i) for i in res.closed_mos ] + active = [ (allMOs[i].eigenvalue,i) for i in res.active_mos ] + virtual =[ (allMOs[i].eigenvalue,i) for i in res.virtual_mos ] + except: + closed = [] + virtual = [] + active = [ (allMOs[i].eigenvalue,i) for i in range(len(allMOs)) ] # closed.sort() # active.sort() @@ -111,117 +163,6 @@ def write_ezfioFile(res,filename): while len(MoMatrix) < len(MOs[0].vector)**2: MoMatrix.append(0.) - ezfio.set_mo_basis_mo_tot_num(mo_tot_num) - ezfio.set_mo_basis_mo_occ(OccNum) - - - res.clean_contractions() - # AO Basis - import string - at = [] - num_prim = [] - magnetic_number = [] - angular_number = [] - power_x = [] - power_y = [] - power_z = [] - coefficient = [] - exponent = [] - res.convert_to_cartesian() - for b in res.basis: - c = b.center - for i,atom in enumerate(res.geometry): - if atom.coord == c: - at.append(i+1) - num_prim.append(len(b.prim)) - s = b.sym - power_x.append( string.count(s,"x") ) - power_y.append( string.count(s,"y") ) - power_z.append( string.count(s,"z") ) - coefficient.append( b.coef ) - exponent.append( [ p.expo for p in b.prim ] ) - ezfio.set_ao_basis_ao_num(len(res.basis)) - ezfio.set_ao_basis_ao_nucl(at) - ezfio.set_ao_basis_ao_prim_num(num_prim) - ezfio.set_ao_basis_ao_power(power_x+power_y+power_z) - prim_num_max = ezfio.get_ao_basis_ao_prim_num_max() - len_res_basis = len(res.basis) - for i in range(len(res.basis)): - coefficient[i] += [ 0. for j in range(len(coefficient[i]),prim_num_max) ] - exponent[i] += [ 0. for j in range(len(exponent[i]),prim_num_max) ] - coefficient = reduce(lambda x, y: x+y, coefficient, []) - exponent = reduce(lambda x, y: x+y, exponent, []) - coef = [] - expo = [] - for i in range(prim_num_max): - for j in range(i,len(coefficient),prim_num_max): - coef.append ( coefficient[j] ) - expo.append ( exponent[j] ) - ezfio.set_ao_basis_ao_coef(coef) - ezfio.set_ao_basis_ao_expo(expo) - ezfio.set_ao_basis_ao_basis("Read by resultsFile") - -# Apply threshold to determinants - if len(res.determinants) == 1: - sorted_determinants = [ (-1.,1.,res.determinants[0]) ] - else: - sorted_determinants = [] - for i,j in zip(res.det_coefficients[0],res.determinants): - sorted_determinants.append((-abs(i),i,j)) - sorted_determinants.sort() - norm = 0.0 - for length, (a,b,c) in enumerate(sorted_determinants): - if -a < det_threshold: - length -=1 - break - norm += a**2 - norm = sqrt(norm) - length += 1 - for i in xrange(length): - a = sorted_determinants[i] - sorted_determinants[i] = (a[0],a[1]/norm,a[2]) - sorted_determinants = sorted_determinants[:length] - -# MOs - mo_tot_num = len(res.mo_sets[MoTag]) - closed_mos = res.closed_mos - active_mos = res.active_mos - virtual_mos = res.virtual_mos - to_remove = [] - to_add = [] - for i in active_mos: - found = False - for (a,b,c) in sorted_determinants: - if i in c['alpha']+c['beta']: - found = True - break - if not found: - to_remove.append(i) - to_add.append(i) - virtual_mos = to_add + virtual_mos - for i in active_mos: - always = True - for (a,b,c) in sorted_determinants: - if not (i in c['alpha'] and i in c['beta']): - always = False - break - if always: - to_remove.append(i) - closed_mos.append(i) - for i in to_remove: - active_mos.remove(i) - - - MOindices = closed_mos + active_mos + virtual_mos - while len(MOindices) < mo_tot_num: - MOindices.append(len(MOindices)) - MOmap = list(MOindices) - for i in range(len(MOindices)): - MOmap[i] = MOindices.index(i) - - - ezfio.set_mo_basis_mo_tot_num(mo_tot_num) - mo = [] for i in MOindices: mo.append(res.mo_sets[MoTag][i]) @@ -234,35 +175,11 @@ def write_ezfioFile(res,filename): while len(mo) < mo_tot_num: mo.append(newmo) Energies = [ m.eigenvalue for m in mo ] - - if res.occ_num is not None: - OccNum = [] - for i in MOindices: - OccNum.append(res.occ_num[MoTag][i]) - - while len(OccNum) < mo_tot_num: - OccNum.append(0.) - ezfio.set_mo_basis_mo_occ(OccNum) - - cls = [ "v" for i in mo ] - for i in closed_mos: - cls[MOmap[i]] = 'c' - for i in active_mos: - cls[MOmap[i]] = 'a' - - sym0 = [ i.sym for i in res.mo_sets[MoTag] ] - sym = [ i.sym for i in res.mo_sets[MoTag] ] - for i in xrange(len(sym)): - sym[MOmap[i]] = sym0[i] - - MoMatrix = [] - for m in mo: - for coef in m.vector: - MoMatrix.append(coef) - while len(MoMatrix) < len(mo[0].vector)**2: - MoMatrix.append(0.) + + ezfio.set_mo_basis_mo_tot_num(mo_tot_num) + ezfio.set_mo_basis_mo_occ(OccNum) ezfio.set_mo_basis_mo_coef(MoMatrix) - del MoMatrix + diff --git a/src/Bielec_integrals/map_integrals.irp.f b/src/Bielec_integrals/map_integrals.irp.f index 1426288a..8fdff799 100644 --- a/src/Bielec_integrals/map_integrals.irp.f +++ b/src/Bielec_integrals/map_integrals.irp.f @@ -368,11 +368,11 @@ subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map) enddo logical :: integral_is_in_map - if (cache_key_kind == 8) then + if (key_kind == 8) then call i8radix_sort(hash,iorder,kk,-1) - else if (cache_key_kind == 4) then + else if (key_kind == 4) then call iradix_sort(hash,iorder,kk,-1) - else if (cache_key_kind == 2) then + else if (key_kind == 2) then call i2radix_sort(hash,iorder,kk,-1) endif diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index e0e04614..2d044ca5 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -123,7 +123,6 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_ call ezfio_has_bitmasks_generators(exists) if (exists) then - print*,'EXIST !!' call ezfio_get_bitmasks_generators(generators_bitmask) else integer :: k, ispin @@ -181,7 +180,6 @@ BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] PROVIDE ezfio_filename call ezfio_has_bitmasks_cas(exists) - print*,'exists = ',exists if (exists) then call ezfio_get_bitmasks_cas(cas_bitmask) else diff --git a/src/CAS_SD_selected/ASSUMPTIONS.rst b/src/CAS_SD/ASSUMPTIONS.rst similarity index 100% rename from src/CAS_SD_selected/ASSUMPTIONS.rst rename to src/CAS_SD/ASSUMPTIONS.rst diff --git a/src/CAS_SD/EZFIO.cfg b/src/CAS_SD/EZFIO.cfg new file mode 100644 index 00000000..5629e90b --- /dev/null +++ b/src/CAS_SD/EZFIO.cfg @@ -0,0 +1,36 @@ +[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" +interface: output + +[energy_pt2] +type: double precision +doc: "Calculated selected CAS-SD energy with PT2 correction" +interface: output + diff --git a/src/CAS_SD_selected/H_apply.irp.f b/src/CAS_SD/H_apply.irp.f similarity index 68% rename from src/CAS_SD_selected/H_apply.irp.f rename to src/CAS_SD/H_apply.irp.f index 76c046ae..85b2f1c6 100644 --- a/src/CAS_SD_selected/H_apply.irp.f +++ b/src/CAS_SD/H_apply.irp.f @@ -2,11 +2,14 @@ use bitmasks BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * -s = H_apply("FCI") +s = H_apply("CAS_SD") +print s + +s = H_apply("CAS_SD_selected") s.set_selection_pt2("epstein_nesbet_2x2") print s -s = H_apply("FCI_PT2") +s = H_apply("CAS_SD_PT2") s.set_perturbation("epstein_nesbet_2x2") print s diff --git a/src/CAS_SD_selected/Makefile b/src/CAS_SD/Makefile similarity index 100% rename from src/CAS_SD_selected/Makefile rename to src/CAS_SD/Makefile diff --git a/src/CAS_SD_selected/NEEDED_MODULES b/src/CAS_SD/NEEDED_MODULES similarity index 100% rename from src/CAS_SD_selected/NEEDED_MODULES rename to src/CAS_SD/NEEDED_MODULES diff --git a/src/CAS_SD/README.rst b/src/CAS_SD/README.rst new file mode 100644 index 00000000..71f828ab --- /dev/null +++ b/src/CAS_SD/README.rst @@ -0,0 +1,9 @@ +====================== +CAS_SD_selected Module +====================== + +Selected CAS + SD module. + +1) Set the different MO classes using the ``qp_set_mo_class`` command +2) Run the selected CAS+SD program + diff --git a/src/CAS_SD_selected/cas_sd.irp.f b/src/CAS_SD/cas_sd.irp.f similarity index 87% rename from src/CAS_SD_selected/cas_sd.irp.f rename to src/CAS_SD/cas_sd.irp.f index e5e0497a..144eec88 100644 --- a/src/CAS_SD_selected/cas_sd.irp.f +++ b/src/CAS_SD/cas_sd.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_cas_sd) 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_cas_sd 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) - call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) + do while (N_det < n_det_max_cas_sd.and.maxval(abs(pt2(1:N_st))) > pt2_max) + call H_apply_CAS_SD(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_cas_sd) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = n_det_max_cas_sd soft_touch N_det psi_det psi_coef endif call diagonalize_CI @@ -60,7 +60,7 @@ program full_ci integer :: exc_max, degree_min exc_max = 0 print *, 'CAS determinants : ', N_det_generators - do i=1,N_det_generators + do i=1,min(N_det_generators,10) do k=i,N_det_generators call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators(1,1,i),degree,N_int) exc_max = max(exc_max,degree) diff --git a/src/CAS_SD/cas_sd_selected.irp.f b/src/CAS_SD/cas_sd_selected.irp.f new file mode 100644 index 00000000..fc9d4dd2 --- /dev/null +++ b/src/CAS_SD/cas_sd_selected.irp.f @@ -0,0 +1,86 @@ +program full_ci + implicit none + integer :: i,k + + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + character*(64) :: perturbation + + pt2 = 1.d0 + diag_algorithm = "Lapack" + if (N_det > n_det_max_cas_sd) then + call diagonalize_CI + call save_wavefunction + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = n_det_max_cas_sd + soft_touch N_det psi_det psi_coef + call diagonalize_CI + call save_wavefunction + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + endif + + do while (N_det < n_det_max_cas_sd.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 + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = n_det_max_cas_sd + soft_touch N_det psi_det psi_coef + endif + call diagonalize_CI + call save_wavefunction + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + call ezfio_set_full_ci_energy(CI_energy) + if (abort_all) then + exit + endif + enddo + + ! Check that it is a CAS-SD + logical :: in_cas + integer :: exc_max, degree_min + exc_max = 0 + print *, 'CAS determinants : ', N_det_generators + do i=1,min(N_det_generators,10) + do k=i,N_det_generators + call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators(1,1,i),degree,N_int) + exc_max = max(exc_max,degree) + enddo + call debug_det(psi_det_generators(1,1,i),N_int) + print *, '' + enddo + print *, 'Max excitation degree in the CAS :', exc_max + do i=1,N_det + in_cas = .False. + degree_min = 1000 + do k=1,N_det_generators + call get_excitation_degree(psi_det_generators(1,1,k),psi_det(1,1,i),degree,N_int) + degree_min = min(degree_min,degree) + enddo + if (degree_min > 2) then + print *, 'Error : This is not a CAS-SD : ' + print *, 'Excited determinant:', degree_min + call debug_det(psi_det(1,1,k),N_int) + stop + endif + enddo +end diff --git a/src/CAS_SD_selected/README.rst b/src/CAS_SD_selected/README.rst deleted file mode 100644 index f66fe229..00000000 --- a/src/CAS_SD_selected/README.rst +++ /dev/null @@ -1,5 +0,0 @@ -====================== -CAS_SD_selected Module -====================== - -Selected CAS + SD module diff --git a/src/CAS_SD_selected/options.irp.f b/src/CAS_SD_selected/options.irp.f deleted file mode 100644 index 553fbe9f..00000000 --- a/src/CAS_SD_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/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index 0b6bc9fd..6d310d95 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -13,6 +13,7 @@ program cisd print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy enddo + call save_wavefunction ! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) ! do i = 1, N_states ! print*,'eigvalues(i) = ',eigvalues(i) diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index e44562e7..a9a282ae 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -1,4 +1,4 @@ -subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc $parameters ) +subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters ) use omp_lib use bitmasks implicit none @@ -14,7 +14,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2) integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2) - integer, intent(in) :: iproc + integer, intent(in) :: iproc_in integer(bit_kind), allocatable :: hole_save(:,:) integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:) integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:) @@ -30,6 +30,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene integer, allocatable :: ia_ja_pairs(:,:,:) integer, allocatable :: ib_jb_pairs(:,:) double precision :: diag_H_mat_elem + integer :: iproc integer(omp_lock_kind), save :: lck, ifirst=0 if (ifirst == 0) then !$ call omp_init_lock(lck) @@ -38,12 +39,13 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene logical :: check_double_excitation check_double_excitation = .True. - + iproc = iproc_in $initialization $omp_parallel +!$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & @@ -248,7 +250,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene $finalization end -subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $parameters ) +subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $parameters ) use omp_lib use bitmasks implicit none @@ -262,7 +264,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param integer ,intent(in) :: i_generator integer(bit_kind),intent(in) :: key_in(N_int,2) integer(bit_kind),intent(in) :: hole_1(N_int,2), particl_1(N_int,2) - integer, intent(in) :: iproc + integer, intent(in) :: iproc_in integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind),allocatable :: hole_save(:,:) integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:) @@ -281,8 +283,11 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param logical, allocatable :: array_pairs(:,:) double precision :: diag_H_mat_elem integer(omp_lock_kind), save :: lck, ifirst=0 + integer :: iproc logical :: check_double_excitation + iproc = iproc_in + check_double_excitation = .True. $check_double_excitation @@ -295,6 +300,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param $initialization $omp_parallel +!$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & @@ -396,7 +402,8 @@ subroutine $subroutine($params_main) integer :: iproc $initialization - PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators + nmax = mod( N_det_generators,nproc ) @@ -406,6 +413,7 @@ subroutine $subroutine($params_main) call wall_time(wall_0) + iproc = 0 allocate( mask(N_int,2,6) ) do i_generator=1,nmax @@ -443,12 +451,12 @@ subroutine $subroutine($params_main) call $subroutine_diexc(psi_det_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & - i_generator, 0 $params_post) + i_generator, iproc $params_post) endif if($do_mono_excitations)then call $subroutine_monoexc(psi_det_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & - i_generator, 0 $params_post) + i_generator, iproc $params_post) endif call wall_time(wall_1) $printout_always @@ -463,7 +471,6 @@ subroutine $subroutine($params_main) !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc) call wall_time(wall_0) - iproc = 0 !$ iproc = omp_get_thread_num() allocate( mask(N_int,2,6) ) !$OMP DO SCHEDULE(dynamic,1) diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index 07331782..3c7eb581 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -95,9 +95,9 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) enddo i += 1 -! if (i > N_det) then -! return -! endif + if (i > N_det) then + return + endif !DIR$ FORCEINLINE do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref) @@ -116,39 +116,39 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) enddo if (is_in_wavefunction) then get_index_in_psi_det_sorted_bit = i - exit -! return +! exit + return endif endif i += 1 if (i > N_det) then - exit -! return +! exit + return endif enddo ! DEBUG is_in_wf - if (is_in_wavefunction) then - degree = 1 - do i=1,N_det - integer :: degree - call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) - if (degree == 0) then - exit - endif - enddo - if (degree /=0) then - stop 'pouet 1' - endif - else - do i=1,N_det - call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) - if (degree == 0) then - stop 'pouet 2' - endif - enddo - endif +! if (is_in_wavefunction) then +! degree = 1 +! do i=1,N_det +! integer :: degree +! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) +! if (degree == 0) then +! exit +! endif +! enddo +! if (degree /=0) then +! stop 'pouet 1' +! endif +! else +! do i=1,N_det +! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) +! if (degree == 0) then +! stop 'pouet 2' +! endif +! enddo +! endif ! END DEBUG is_in_wf end diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 0a0ec864..00e683fc 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -467,33 +467,52 @@ END_PROVIDER ! to accelerate the search of a random determinant in the wave ! function. END_DOC + + call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, & + psi_det_sorted_bit, psi_coef_sorted_bit) +END_PROVIDER + +subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out) + use bitmasks + implicit none + integer, intent(in) :: Ndet + integer(bit_kind), intent(in) :: det_in (N_int,2,psi_det_size) + double precision , intent(in) :: coef_in(psi_det_size,N_states) + integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size) + double precision , intent(out) :: coef_out(psi_det_size,N_states) + BEGIN_DOC + ! Determinants are sorted are sorted according to their det_search_key. + ! Useful to accelerate the search of a random determinant in the wave + ! function. + END_DOC integer :: i,j,k integer, allocatable :: iorder(:) integer*8, allocatable :: bit_tmp(:) integer*8, external :: det_search_key - allocate ( iorder(N_det), bit_tmp(N_det) ) + allocate ( iorder(Ndet), bit_tmp(Ndet) ) - do i=1,N_det + do i=1,Ndet iorder(i) = i !$DIR FORCEINLINE - bit_tmp(i) = det_search_key(psi_det(1,1,i),N_int) + bit_tmp(i) = det_search_key(det_in(1,1,i),N_int) enddo - call i8sort(bit_tmp,iorder,N_det) + call i8sort(bit_tmp,iorder,Ndet) !DIR$ IVDEP - do i=1,N_det + do i=1,Ndet do j=1,N_int - psi_det_sorted_bit(j,1,i) = psi_det(j,1,iorder(i)) - psi_det_sorted_bit(j,2,i) = psi_det(j,2,iorder(i)) + det_out(j,1,i) = det_in(j,1,iorder(i)) + det_out(j,2,i) = det_in(j,2,iorder(i)) enddo do k=1,N_states - psi_coef_sorted_bit(i,k) = psi_coef(iorder(i),k) + coef_out(i,k) = coef_in(iorder(i),k) enddo enddo deallocate(iorder, bit_tmp) -END_PROVIDER +end + subroutine int_of_3_highest_electrons( det_in, res, Nint ) implicit none diff --git a/src/Dets/psi_cas.irp.f b/src/Dets/psi_cas.irp.f new file mode 100644 index 00000000..299e5e8f --- /dev/null +++ b/src/Dets/psi_cas.irp.f @@ -0,0 +1,114 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_cas, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_cas ] + implicit none + BEGIN_DOC + ! CAS wave function, defined from the application of the CAS bitmask on the + ! determinants. idx_cas gives the indice of the CAS determinant in psi_det. + END_DOC + integer :: i, k, l + logical :: good + N_det_cas = 0 + do i=1,N_det + do l=1,n_cas_bitmask + good = .True. + do k=1,N_int + good = good .and. ( & + iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & + iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & + iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & + iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) + enddo + if (good) then + exit + endif + enddo + if (good) then + N_det_cas = N_det_cas+1 + do k=1,N_int + psi_cas(k,1,N_det_cas) = psi_det(k,1,i) + psi_cas(k,2,N_det_cas) = psi_det(k,2,i) + enddo + idx_cas(N_det_cas) = i + do k=1,N_states + psi_cas_coef(N_det_cas,k) = psi_coef(i,k) + enddo + endif + enddo + call write_int(output_dets,N_det_cas, 'Number of determinants in the CAS') + +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! CAS determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, & + psi_cas_sorted_bit, psi_cas_coef_sorted_bit) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_cas_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_non_cas, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_non_cas ] + implicit none + BEGIN_DOC + ! Set of determinants which are not part of the CAS, defined from the application + ! of the CAS bitmask on the determinants. + ! idx_non_cas gives the indice of the determinant in psi_det. + END_DOC + integer :: i_non_cas,j,k + integer :: degree + logical :: in_cas + i_non_cas =0 + do k=1,N_det + in_cas = .False. + do j=1,N_det_cas + call get_excitation_degree(psi_cas(1,1,j), psi_det(1,1,k), degree, N_int) + if (degree == 0) then + in_cas = .True. + exit + endif + enddo + if (.not.in_cas) then + double precision :: hij + i_non_cas += 1 + do j=1,N_int + psi_non_cas(j,1,i_non_cas) = psi_det(j,1,k) + psi_non_cas(j,2,i_non_cas) = psi_det(j,2,k) + enddo + do j=1,N_states + psi_non_cas_coef(i_non_cas,j) = psi_coef(k,j) + enddo + idx_non_cas(i_non_cas) = k + endif + enddo + N_det_non_cas = i_non_cas +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_cas_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! CAS determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, & + psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit) + +END_PROVIDER + + + + + diff --git a/src/Generators_CAS/generators.irp.f b/src/Generators_CAS/generators.irp.f index 5dcfc1b0..511877a0 100644 --- a/src/Generators_CAS/generators.irp.f +++ b/src/Generators_CAS/generators.irp.f @@ -62,7 +62,6 @@ END_PROVIDER psi_det_generators(k,2,m) = psi_det(k,2,i) enddo psi_coef_generators(m,:) = psi_coef(m,:) -! call debug_det(psi_det_generators(1,1,m),N_int) endif enddo diff --git a/src/MRCC/EZFIO.cfg b/src/MRCC/EZFIO.cfg new file mode 100644 index 00000000..ff586985 --- /dev/null +++ b/src/MRCC/EZFIO.cfg @@ -0,0 +1,4 @@ +[energy] +type: double precision +doc: Calculated MRCC energy +interface: output \ No newline at end of file diff --git a/src/MRCC/H_apply.irp.f b/src/MRCC/H_apply.irp.f index 4f284112..cadaca32 100644 --- a/src/MRCC/H_apply.irp.f +++ b/src/MRCC/H_apply.irp.f @@ -2,13 +2,13 @@ use bitmasks BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * -s = H_apply("mrcc") +s = H_apply("mrcc_simple") s.data["parameters"] = ", delta_ij_sd_, Ndet_sd" s.data["declarations"] += """ integer, intent(in) :: Ndet_sd double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) """ -s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" +s.data["keys_work"] = "call mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" s.data["params_post"] += ", delta_ij_sd_, Ndet_sd" s.data["params_main"] += "delta_ij_sd_, Ndet_sd" s.data["decls_main"] += """ @@ -18,8 +18,14 @@ s.data["decls_main"] += """ s.data["finalization"] = "" s.data["copy_buffer"] = "" s.data["generate_psi_guess"] = "" -s.data["size_max"] = "256" +s.data["size_max"] = "3072" +print s + + + +s.data["subroutine"] = "H_apply_mrcc" +s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" print s END_SHELL diff --git a/src/MRCC/NEEDED_MODULES b/src/MRCC/NEEDED_MODULES index 3711786f..dd64048c 100644 --- a/src/MRCC/NEEDED_MODULES +++ b/src/MRCC/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask CAS_SD_selected Dets Electrons Ezfio_files Generators_CAS Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils +AOs Bielec_integrals Bitmask CAS_SD Dets Electrons Ezfio_files Generators_CAS Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst index 1fbcccdb..3720884a 100644 --- a/src/MRCC/README.rst +++ b/src/MRCC/README.rst @@ -1,14 +1,82 @@ --4.142795384334731 +E(CAS+SD) = -3.93510180989509 + + -0.726935163332 + ++++---- + ++++---- + + 0.685482292841 + +++-+--- + +++-+--- + + 0.0409791961003 + ++-+-+-- + ++-+-+-- + +c2*c3/c1 = -0.038642391672 + + + -0.996413146877 + ++++---- + ++++---- + + -0.0598366139154 + ++-+-+-- + +++-+--- + + -0.0598366139154 + +++-+--- + ++-+-+-- + + +======= + + 0.706616635698 + +++-+--- + +++-+--- + + -0.692594178414 + ++++---- + ++++---- + + -0.0915716740265 + ++-+-+-- + +++-+--- + + -0.0915716740265 + +++-+--- + ++-+-+-- + + 0.0461418323107 + ++-+-+-- + ++-+-+-- + + -0.0458957789452 + ++--++-- + ++--++-- + + +---- + + 0.713102867186 + ++++---- + ++++---- + + -0.688067688244 + +++-+--- + +++-+--- + + 0.089262694488 + +++-+--- + ++-+-+-- + + 0.089262694488 + ++-+-+-- + +++-+--- + + -0.0459510603899 + ++-+-+-- + ++-+-+-- -4.695183071437694E-002 - Determinant 64 - --------------------------------------------- - 000000000000002E|000000000000002E - |-+++-+----------------------------------------------------------| - |-+++-+----------------------------------------------------------| -CAS-SD: -4.14214374069306 - : -4.14230904320551 -E0 = -11.5634986758976 diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f index 663b9f5a..2b879882 100644 --- a/src/MRCC/mrcc.irp.f +++ b/src/MRCC/mrcc.irp.f @@ -3,41 +3,65 @@ program mrcc read_wf = .True. TOUCH read_wf call run + call run_mrcc +! call run_mrcc_test end subroutine run implicit none - print *, N_det - print *, N_det_cas - print *, N_det_sd - -! call update_generators - integer :: i + integer :: i,j print *, 'CAS' print *, '===' do i=1,N_det_cas - print *, psi_cas_coefs(i,:) + print *, psi_cas_coef(i,:) call debug_det(psi_cas(1,1,i),N_int) enddo ! print *, 'SD' ! print *, '==' -! do i=1,N_det_sd -! print *, psi_sd_coefs(i,:) -! call debug_det(psi_sd(1,1,i),N_int) -! enddo -! - print *, 'MRCC' - print *, '====' - print *, ci_energy(:) - print *, h_matrix_all_dets(3,3), delta_ij(3,3,1) - print *, h_matrix_all_dets(3,3), delta_ij(3,3,1) - print *, ci_energy_dressed(:) -! print *, 'max', maxval(delta_ij_sd) -! print *, 'min', minval(delta_ij_sd) -! -! do i=1,N_det -! print '(10(F10.6,X))', delta_ij(i,1:N_det,1) +! do i=1,N_det_non_cas +! print *, psi_non_cas_coef(i,:) +! call debug_det(psi_non_cas(1,1,i),N_int) ! enddo + call write_double(6,ci_energy(1),"Initial CI energy") +end +subroutine run_mrcc_test + implicit none + integer :: i,j + double precision :: pt2 + pt2 = 0.d0 + do j=1,N_det + do i=1,N_det + pt2 += psi_coef(i,1)*psi_coef(j,1) * delta_ij(i,j,1) + enddo + enddo + print *, ci_energy(1) + print *, ci_energy(1)+pt2 +end +subroutine run_mrcc + implicit none + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + do while (delta_E > 1.d-8) + iteration += 1 + print *, '===========================' + print *, 'MRCC Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCC energy") + call diagonalize_ci_dressed + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + call ezfio_set_mrcc_energy(ci_energy_dressed(1)) +! call save_wavefunction + end diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index e5c8d233..083cea5a 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -1,15 +1,169 @@ -subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) +use omp_lib + +BEGIN_PROVIDER [ integer(omp_lock_kind), psi_cas_lock, (psi_det_size) ] + implicit none + BEGIN_DOC + ! Locks on CAS determinants to fill delta_ij + END_DOC + integer :: i + do i=1,psi_det_size + call omp_init_lock( psi_cas_lock(i) ) + enddo + +END_PROVIDER + +subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc - integer, intent(in) :: Ndet_sd - double precision, intent(inout) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) + integer, intent(in) :: Ndet + double precision, intent(inout) :: delta_ij_(Ndet,Ndet,*) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,l + integer :: degree_alpha(psi_det_size) + integer :: idx_alpha(0:psi_det_size) + logical :: good + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq, c_ref ,degree + integer :: connected_to_ref + + double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision, allocatable :: dIa_hla(:,:) + double precision :: haj, phase, phase2 + double precision :: f(N_states), ci_inv(N_states) + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + integer(bit_kind) :: tmp_det(Nint,2) + integer :: iint, ipos + integer :: i_state, k_sd, l_sd, i_I, i_alpha + + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) + + allocate (dIa_hla(N_states,Ndet)) + + ! |I> + + ! |alpha> + do i_alpha=1,N_tq + call get_excitation_degree_vector(psi_non_cas,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_cas,idx_alpha) + + ! |I> + do i_I=1,N_det_cas + ! Find triples and quadruple grand parents + call get_excitation_degree(tq(1,1,i_alpha),psi_cas(1,1,i_I),degree,Nint) + if (degree > 4) then + cycle + endif + + do i_state=1,N_states + dIa(i_state) = 0.d0 + enddo + + ! |alpha> + do k_sd=1,idx_alpha(0) + call get_excitation_degree(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),degree,Nint) + if (degree > 2) then + cycle + endif + ! + ! + call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),Nint,hIk) + do i_state=1,N_states + dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + enddo + ! |l> = Exc(k -> alpha) |I> + call get_excitation(psi_non_cas(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do k=1,N_int + tmp_det(k,1) = psi_cas(k,1,i_I) + tmp_det(k,2) = psi_cas(k,2,i_I) + enddo + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + if (degree_alpha(k_sd) == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + endif + + ! + do i_state=1,N_states + dka(i_state) = 0.d0 + enddo + do l_sd=k_sd+1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_cas(1,1,idx_alpha(l_sd)),degree,Nint) + if (degree == 0) then + call get_excitation(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) + call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hIl) + do i_state=1,N_states + dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + enddo + exit + endif + enddo + do i_state=1,N_states + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) + enddo + enddo + + do i_state=1,N_states + ci_inv(i_state) = 1.d0/psi_cas_coef(i_I,i_state) + enddo + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + call i_h_j(tq(1,1,i_alpha),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hla) + do i_state=1,N_states + dIa_hla(i_state,k_sd) = dIa(i_state) * hla + enddo + enddo + call omp_set_lock( psi_cas_lock(i_I) ) + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + do i_state=1,N_states + delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),i_state) += dIa_hla(i_state,k_sd) + delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),i_state) += dIa_hla(i_state,k_sd) + delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state) + enddo + enddo + call omp_unset_lock( psi_cas_lock(i_I) ) + enddo + enddo + deallocate (dIa_hla) +end + + + + + + + +subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,i_generator,n_selected,det_buffer,Nint,iproc) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet_non_cas + double precision, intent(inout) :: delta_ij_non_cas_(Ndet_non_cas,Ndet_non_cas,*) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,m integer :: new_size - logical :: is_in_wavefunction integer :: degree(psi_det_size) integer :: idx(0:psi_det_size) logical :: good @@ -18,6 +172,51 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin integer :: N_tq, c_ref integer :: connected_to_ref + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) + + ! Compute / (E0 - Haa) + double precision :: hka, haa + double precision :: haj + double precision :: f(N_states) + + do i=1,N_tq + call get_excitation_degree_vector(psi_non_cas,tq(1,1,i),degree,Nint,Ndet_non_cas,idx) + call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa) + do m=1,N_states + f(m) = 1.d0/(ci_electronic_energy(m)-haa) + enddo + do k=1,idx(0) + call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(k)),Nint,hka) + do j=k,idx(0) + call i_h_j(tq(1,1,i),psi_non_cas(1,1,idx(j)),Nint,haj) + do m=1,N_states + delta_ij_non_cas_(idx(k), idx(j),m) += haj*hka* f(m) + delta_ij_non_cas_(idx(j), idx(k),m) += haj*hka* f(m) + enddo + enddo + enddo + enddo +end + + +subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + logical :: is_in_wavefunction + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) + logical :: good + + integer(bit_kind), intent(out) :: tq(Nint,2,n_selected) + integer, intent(out) :: N_tq + integer :: c_ref + integer :: connected_to_ref + N_tq = 0 do i=1,N_selected c_ref = connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, & @@ -48,143 +247,11 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin endif enddo - ! Compute / (E0 - Haa) - double precision :: hka, haa - double precision :: haj - double precision :: f(N_states) - - do i=1,N_tq - call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx) - call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa) - do m=1,N_states - f(m) = 1.d0/(ci_electronic_energy(m)-haa) - enddo - do k=1,idx(0) - call i_h_j(tq(1,1,i),psi_sd(1,1,idx(k)),Nint,hka) - do j=k,idx(0) - call i_h_j(tq(1,1,i),psi_sd(1,1,idx(j)),Nint,haj) - do m=1,N_states - delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m) - delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m) - enddo - enddo - enddo - enddo end -BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ] - implicit none - BEGIN_DOC - ! Dressing matrix in SD basis - END_DOC - delta_ij_sd = 0.d0 - call H_apply_mrcc(delta_ij_sd,N_det_sd) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] - implicit none - BEGIN_DOC - ! Dressing matrix in N_det basis - END_DOC - integer :: i,j,m - delta_ij = 0.d0 - do m=1,N_states - do j=1,N_det_sd - do i=1,N_det_sd - delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m) - enddo - enddo - enddo -END_PROVIDER - -BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] - implicit none - BEGIN_DOC - ! Dressed H with Delta_ij - END_DOC - integer :: i, j - do j=1,N_det - do i=1,N_det - h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1) - if (i==j) then - print *, i, delta_ij(i,j,1), h_matrix_all_dets(i,j) - endif - enddo - enddo - -END_PROVIDER - BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! Eigenvectors/values of the CI matrix - END_DOC - integer :: i,j - - do j=1,N_states_diag - do i=1,N_det - CI_eigenvectors_dressed(i,j) = psi_coef(i,j) - enddo - enddo - - 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_Dets) - - else if (diag_algorithm == "Lapack") then - - double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) - allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) - allocate (eigenvalues(N_det)) - call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_dressed,size(H_matrix_dressed,1),N_det) - CI_electronic_energy_dressed(:) = 0.d0 - do i=1,N_det - CI_eigenvectors_dressed(i,1) = eigenvectors(i,1) - enddo - integer :: i_state - double precision :: s2 - i_state = 0 - do j=1,N_det - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - if(dabs(s2-expected_s2).le.0.3d0)then - i_state += 1 - do i=1,N_det - CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(i_state) = eigenvalues(j) - CI_eigenvectors_s2_dressed(i_state) = s2 - endif - if (i_state.ge.N_states_diag) then - exit - endif - enddo - deallocate(eigenvectors,eigenvalues) - endif - -END_PROVIDER -BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! N_states lowest eigenvalues of the dressed CI matrix - END_DOC - - integer :: j - character*(8) :: st - call write_time(output_Dets) - do j=1,N_states_diag - CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion - write(st,'(I4)') j - call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st)) - call write_double(output_Dets,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) - enddo -END_PROVIDER + diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 2c5dc811..eca1a706 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -1,150 +1,157 @@ - use bitmasks -BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] - implicit none - BEGIN_DOC - ! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference) - END_DOC - logical :: exists - integer :: i - PROVIDE ezfio_filename - - call ezfio_has_bitmasks_cas(exists) - if (exists) then - call ezfio_get_bitmasks_cas(cas_bitmask) - else - do i=1,N_cas_bitmask - cas_bitmask(:,:,i) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - enddo - endif - -END_PROVIDER - -BEGIN_PROVIDER [ integer, N_det_cas ] - implicit none - BEGIN_DOC - ! Number of generator detetrminants - END_DOC - integer :: i,k,l - logical :: good - call write_time(output_dets) - N_det_cas = 0 - do i=1,N_det - do l=1,n_cas_bitmask - good = .True. - do k=1,N_int - good = good .and. ( & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) - enddo - if (good) then - exit - endif - enddo - if (good) then - N_det_cas += 1 - endif - enddo - N_det_cas = max(N_det_cas, 1) - call write_int(output_dets,N_det_cas, 'Number of determinants in the CAS') -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_cas) ] -&BEGIN_PROVIDER [ double precision, psi_cas_coefs, (N_det_cas,n_states) ] -&BEGIN_PROVIDER [ integer, idx_cas, (N_det_cas) ] - implicit none - BEGIN_DOC - ! For Single reference wave functions, the generator is the - ! Hartree-Fock determinant - END_DOC - integer :: i, k, l, m - logical :: good - m=0 - do i=1,N_det - do l=1,n_cas_bitmask - good = .True. - do k=1,N_int - good = good .and. ( & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) - enddo - if (good) then - exit - endif - enddo - if (good) then - m = m+1 - do k=1,N_int - psi_cas(k,1,m) = psi_det(k,1,i) - psi_cas(k,2,m) = psi_det(k,2,i) - enddo - idx_cas(m) = i - do k=1,N_states - psi_cas_coefs(m,k) = psi_coef(i,k) - enddo - endif - enddo - -END_PROVIDER - - - - - BEGIN_PROVIDER [ integer(bit_kind), psi_sd, (N_int,2,N_det) ] -&BEGIN_PROVIDER [ double precision, psi_sd_coefs, (N_det,n_states) ] -&BEGIN_PROVIDER [ integer, idx_sd, (N_det) ] -&BEGIN_PROVIDER [ integer, N_det_sd] - implicit none - BEGIN_DOC - ! SD - END_DOC - integer :: i_sd,j,k - integer :: degree - logical :: in_cas - i_sd =0 - do k=1,N_det - in_cas = .False. - do j=1,N_det_cas - call get_excitation_degree(psi_cas(1,1,j), psi_det(1,1,k), degree, N_int) - if (degree == 0) then - in_cas = .True. - exit - endif - enddo - if (.not.in_cas) then - double precision :: hij - i_sd += 1 - psi_sd(1:N_int,1:2,i_sd) = psi_det(1:N_int,1:2,k) - psi_sd_coefs(i_sd,1:N_states) = psi_coef(k,1:N_states) - idx_sd(i_sd) = k - endif - enddo - N_det_sd = i_sd -END_PROVIDER - -BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] +BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] implicit none BEGIN_DOC ! cm/ END_DOC integer :: i,k double precision :: ihpsi(N_states) - do i=1,N_det_sd - call i_h_psi(psi_sd(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, & - size(psi_cas_coefs,1), n_states, ihpsi) + do i=1,N_det_non_cas + call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coef, N_int, N_det_cas, & + size(psi_cas_coef,1), n_states, ihpsi) double precision :: hij do k=1,N_states - if (dabs(ihpsi(k)) < 1.d-6) then - lambda_mrcc(i,k) = 0.d0 - else - lambda_mrcc(i,k) = psi_sd_coefs(i,k)/ihpsi(k) - lambda_mrcc(i,k) = min( lambda_mrcc (i,k),0.d0 ) - endif + lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k) + lambda_mrcc(k,i) = min( lambda_mrcc (k,i),0.d0 ) enddo enddo END_PROVIDER + + +BEGIN_PROVIDER [ character*(32), dressing_type ] + implicit none + BEGIN_DOC + ! [ Simple | MRCC ] + END_DOC + dressing_type = "MRCC" +END_PROVIDER + +BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ] + implicit none + BEGIN_DOC + ! Dressing matrix in SD basis + END_DOC + delta_ij_non_cas = 0.d0 + call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] + implicit none + BEGIN_DOC + ! Dressing matrix in N_det basis + END_DOC + integer :: i,j,m + delta_ij = 0.d0 + if (dressing_type == "MRCC") then + call H_apply_mrcc(delta_ij,N_det) + else if (dressing_type == "Simple") then + do m=1,N_states + do j=1,N_det_non_cas + do i=1,N_det_non_cas + delta_ij(idx_non_cas(i),idx_non_cas(j),m) = delta_ij_non_cas(i,j,m) + enddo + enddo + enddo + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] + implicit none + BEGIN_DOC + ! Dressed H with Delta_ij + END_DOC + integer :: i, j + do j=1,N_det + do i=1,N_det + h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] + implicit none + BEGIN_DOC + ! Eigenvectors/values of the CI matrix + END_DOC + integer :: i,j + + do j=1,N_states_diag + do i=1,N_det + CI_eigenvectors_dressed(i,j) = psi_coef(i,j) + enddo + enddo + + 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_Dets) + + else if (diag_algorithm == "Lapack") then + + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_dressed,size(H_matrix_dressed,1),N_det) + CI_electronic_energy_dressed(:) = 0.d0 + do i=1,N_det + CI_eigenvectors_dressed(i,1) = eigenvectors(i,1) + enddo + integer :: i_state + double precision :: s2 + i_state = 0 + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(s2-expected_s2).le.0.3d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(i_state) = eigenvalues(j) + CI_eigenvectors_s2_dressed(i_state) = s2 + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the dressed CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(output_Dets) + do j=1,N_states_diag + CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion + enddo + +END_PROVIDER + +subroutine diagonalize_CI_dressed + implicit none + BEGIN_DOC +! Replace the coefficients of the CI states by the coefficients of the +! eigenstates of the CI matrix + END_DOC + integer :: i,j + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + enddo + enddo + SOFT_TOUCH psi_coef + +end diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index a144c42c..efe8b8f8 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 Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected MRCC +AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets 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 diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 7d173246..a7462f94 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -451,3 +451,6 @@ subroutine set_zero_extra_diag(i1,i2,matrix,lda,m) end + + +