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

Merge pull request #6 from TApplencourt/master

Browse Source 
Merge thomas
This commit is contained in:
Anthony Scemama 2015-04-10 19:07:50 +02:00
commit b4f7bc3724
36 changed files with 1364 additions and 755 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
INSTALL.md
View File

@ -6,8 +6,9 @@
* m4
* GNU make
* Fortran compiler (ifort or gfortran are tested)
* Python >= 2.7
* Python >= 2.6
* Bash
* Patch (for opam)
## Standard installation

66
data/ezfio_defaults/WILL_BE_DELETED.ezfio_default Normal file
View File

@ -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

6
data/ezfio_defaults/q_package.ezfio_default
View File

@ -6,6 +6,12 @@ cis_dressed
standard_doubles True
en_2_2 false
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

177
ocaml/test_input.ml
View File

@ -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 ();

333
scripts/ezfio_interface.py Executable file
View File

@ -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()

60
scripts/ezfio_interface/ei_handler.py
View File

@ -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,7 +799,11 @@ if __name__ == "__main__":
# O c a m l #
# ~#~#~#~#~#~#
if do_all or arguments["--ocaml"]:
try:
str_ocaml_input = create_ocaml_input(dict_ezfio_cfg, module_lower)
except ValueError:
pass
else:
save_ocaml_input(module_lower, str_ocaml_input)
# ~#~#~#~#~#~#~#~#~#~#~#~#~ #

24
scripts/ezfio_with_default.py
View File

@ -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 = open(filename,'r')
lines = file.readlines()
file.close()
k=-1
# Search directory
for k,line in enumerate(lines):
if line[0] != ' ':
@ -120,18 +115,21 @@ END_PROVIDER
break
v = buffer[1]
name = self.name
true = True
false= False
try:
true = True
false = False
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

2
scripts/generate_h_apply.py
View File

@ -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"

193
scripts/qp_convert.py
View File

@ -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]
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])
@ -235,34 +176,10 @@ def write_ezfioFile(res,filename):
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_tot_num(mo_tot_num)
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_coef(MoMatrix)
del MoMatrix

6
src/Bielec_integrals/map_integrals.irp.f
View File

@ -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

2
src/Bitmask/bitmasks.irp.f
View File

@ -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

0
src/CAS_SD_selected/ASSUMPTIONS.rst → src/CAS_SD/ASSUMPTIONS.rst
View File

36
src/CAS_SD/EZFIO.cfg Normal file
View File

@ -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

7
src/CAS_SD_selected/H_apply.irp.f → src/CAS_SD/H_apply.irp.f
View File

@ -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

0
src/CAS_SD_selected/Makefile → src/CAS_SD/Makefile
View File

0
src/CAS_SD_selected/NEEDED_MODULES → src/CAS_SD/NEEDED_MODULES
View File

9
src/CAS_SD/README.rst Normal file
View File

@ -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

14
src/CAS_SD_selected/cas_sd.irp.f → src/CAS_SD/cas_sd.irp.f
View File

@ -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)

87
src/CAS_SD/cas_sd_selected.irp.f Normal file
View File

@ -0,0 +1,87 @@
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
PROVIDE N_det_cas
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_cas
do i=1,min(N_det_cas,10)
do k=i,N_det_cas
call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
exc_max = max(exc_max,degree)
enddo
call debug_det(psi_cas(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_cas
call get_excitation_degree(psi_cas(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

5
src/CAS_SD_selected/README.rst
View File

@ -1,5 +0,0 @@
======================
CAS_SD_selected Module
======================
Selected CAS + SD module

32
src/CAS_SD_selected/options.irp.f
View File

@ -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

4
src/CISD/cisd.irp.f
View File

@ -14,4 +14,8 @@ program cisd
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)
! enddo
end

25
src/Dets/H_apply_template.f
View File

@ -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)

54
src/Dets/connected_to_ref.irp.f
View File

@ -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

37
src/Dets/determinants.irp.f
View File

@ -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

114
src/Dets/psi_cas.irp.f Normal file
View File

@ -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

1
src/Generators_CAS/generators.irp.f
View File

@ -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

4
src/MRCC/EZFIO.cfg Normal file
View File

@ -0,0 +1,4 @@
[energy]
type: double precision
doc: Calculated MRCC energy
interface: output

12
src/MRCC/H_apply.irp.f
View File

@ -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

2
src/MRCC/NEEDED_MODULES
View File

@ -1,2 +1,2 @@
AOs Bielec_integrals 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 Dets Electrons Ezfio_files Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils

88
src/MRCC/README.rst
View File

@ -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

73
src/MRCC/mrcc.irp.f
View File

@ -3,40 +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
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)
! print *, 'SD'
! print *, '=='
! 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

339
src/MRCC/mrcc_dress.irp.f
View File

@ -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
! <I| <> |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
! <I| /k\ |alpha>
! <I|H|k>
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
! <I| \l/ |alpha>
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 <k|H|a><a|H|j> / (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 <k|H|a><a|H|j> / (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

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

@ -1,150 +1,161 @@
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/<Psi_0|H|D_m>
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
if (dabs(ihpsi(k)) > 1.d-5) then
lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k)
lambda_mrcc(k,i) = min( lambda_mrcc (k,i),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 )
lambda_mrcc(k,i) = 0.d0
endif
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

2
src/NEEDED_MODULES
View File

@ -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

3
src/Utils/LinearAlgebra.irp.f
View File

@ -451,3 +451,6 @@ subroutine set_zero_extra_diag(i1,i2,matrix,lda,m)
end