LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-22 05:02:15 +02:00
Code Issues Releases Wiki Activity

Merge branch 'master' of github.com:LCPQ/quantum_package

Browse Source 
Conflicts:
	src/BiInts/ao_bi_integrals.irp.f
This commit is contained in:
Anthony Scemama 2014-05-16 23:48:47 +02:00
commit c6d740bb89
68 changed files with 5702 additions and 373 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
scripts/generate_h_apply.py Executable file
View File

@ -0,0 +1,65 @@
#!/usr/bin/env python
import os
file = open(os.environ["QPACKAGE_ROOT"]+'/src/Dets/H_apply_template.f','r')
template = file.read()
file.close()
keywords = """
subroutine
parameters
initialization
declarations
keys_work
finalization
""".split()
def new_dict(openmp=True):
s ={}
for k in keywords:
s[k] = ""
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, &
!$OMP occ_particle,occ_hole,j_a,k_a,other_spin, &
!$OMP hole_save,ispin,jj,l_a,hij_elec,hij_tab, &
!$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, &
!$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,&
!$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, &
!$OMP N_elec_in_key_hole_2,ia_ja_pairs) &
!$OMP SHARED(key_in,N_int,elec_num_tab, &
!$OMP hole_1, particl_1, hole_2, particl_2, &
!$OMP lck,thresh,elec_alpha_num,E_ref)"""
s["omp_init_lock"] = "call omp_init_lock(lck)"
s["omp_set_lock"] = "call omp_set_lock(lck)"
s["omp_unset_lock"] = "call omp_unset_lock(lck)"
s["omp_test_lock"] = "omp_test_lock(lck)"
s["omp_destroy_lock"] = "call omp_destroy_lock(lck)"
s["omp_end_parallel"] = "!$OMP END PARALLEL"
s["omp_master"] = "!$OMP MASTER"
s["omp_end_master"] = "!$OMP END MASTER"
s["omp_barrier"] = "!$OMP BARRIER"
s["omp_do"] = "!$OMP DO SCHEDULE (static)"
s["omp_enddo"] = "!$OMP ENDDO NOWAIT"
if not openmp:
for k in s:
s[k] = ""
s["omp_test_lock"] = ".False."
s["size_max"] = str(1024*128)
s["set_i_H_j_threshold"] = """
thresh = H_apply_threshold
"""
return s
def create_h_apply(s):
buffer = template
for key in s:
buffer = buffer.replace('$'+key, s[key])
print buffer

110
scripts/update_README.py
View File

@ -8,14 +8,30 @@ __author__ = "Anthony Scemama <scemama@irsamc.ups-tlse.fr>"
README="README.rst"
Assum_key="Assumptions\n===========\n"
Needed_key="Needed Modules\n==============\n"
Doc_key="Documentation\n=============\n"
Sentinel="@@$%&@@"
URL="http://github.com/LCPQ/quantum_package/tree/master/src/"
import os
import subprocess
header = """
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
"""
try:
subprocess.check_output("git status".split())
has_git = True
except OSError:
has_git = False
def fetch_splitted_data():
"""Read the README.rst file and split it in 3 strings:
"""Read the README.rst file and split it in strings:
* The description
* The assumptions
* The documentation
* The needed modules
The result is given as a list of strings
"""
@ -26,6 +42,7 @@ def fetch_splitted_data():
# Place sentinels
data = data.replace(Assum_key,Sentinel+Assum_key)
data = data.replace(Doc_key,Sentinel+Doc_key)
data = data.replace(Needed_key,Sentinel+Needed_key)
# Now Split data using the sentinels
@ -46,11 +63,6 @@ def update_assumptions(data):
assumptions = file.read()
file.close()
header = """
.. Do not edit this section. It was auto-generated from the
.. ASSUMPTIONS.rst file.
"""
if assumptions.strip() != "":
assumptions = Assum_key + header + assumptions + '\n\n'
@ -74,11 +86,6 @@ def update_needed(data):
modules = file.read()
file.close()
header = """
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
"""
if modules.strip() != "":
modules = [ '* `%s <%s%s>`_'%(x,URL,x) for x in modules.split() ]
modules = "\n".join(modules)
@ -95,23 +102,92 @@ def update_needed(data):
return data
import subprocess
def update_documentation(data):
"""Reads the BEGIN_DOC ... END_DOC blocks and builds the documentation"""
# If the file does not exist, don't do anything
try:
file = open('tags','r')
except:
return
tags = file.readlines()
file.close()
def extract_doc(item):
"""Extracts the documentation contained in IRPF90_man file"""
file = open("IRPF90_man/%s.l"%(item),'r')
lines = file.readlines()
file.close()
result = []
inside = False
for line in lines:
if not inside:
inside = line.startswith(".SH Description")
else:
if line.startswith(".SH"):
return "".join(result)
result.append(" "+line.strip()+"\n")
items = []
command = "git ls-tree --full-tree --name-only HEAD:src/%s"
command = command%(os.path.basename(os.getcwd()))
try:
tracked_files = subprocess.check_output(command.split())
tracked_files = tracked_files.splitlines()
except:
tracked_files = []
for filename in tracked_files:
if filename.endswith('.irp.f'):
# Search for providers, subroutines and functions in each file using
# the tags file
search = "\t"+filename+"\t"
tmp = filter(lambda line: search in line, tags)
# Search for the documentation in the IRPF90_man directory
for item in tmp :
item, _, line = item.strip().split('\t')
doc = extract_doc(item)
items.append( (item, filename, doc, line) )
dirname = os.path.basename(os.getcwd())
# Write the documentation in the README
template = "`%(item)s <%(url)s%(dirname)s/%(filename)s#L%(line)s>`_\n%(doc)s\n"
documentation = Doc_key + header
url = URL
for item, filename, doc, line in items:
documentation += template%locals()
documentation += '\n\n'
has_doc = False
for i in range(len(data)):
if data[i].startswith(Doc_key):
has_doc = True
data[i] = documentation
if not has_doc:
data.append(documentation)
return data
def git_add():
"""Executes:
git add README.rst
if git is present on the machine."""
command = "git add "+README
try:
subprocess.call(command.split())
except OSError:
pass
os.system(command+" &> /dev/null")
def main():
if not has_git:
return
data = fetch_splitted_data()
data = update_assumptions(data)
data = update_documentation(data)
data = update_needed(data)
output = ''.join(data)

1
setup_environment.sh
View File

@ -22,6 +22,7 @@ fi
cat << EOF > quantum_package.rc
export IRPF90=${IRPF90}
export QPACKAGE_ROOT=${QPACKAGE_ROOT}
export PYTHONPATH+=:\${QPACKAGE_ROOT}/scripts
export PATH+=:\${QPACKAGE_ROOT}/scripts
export PATH+=:\${QPACKAGE_ROOT}/bin
export QPACKAGE_CACHE_URL="http://qmcchem.ups-tlse.fr/files/scemama/quantum_package/cache"

3
src/AOs/NEEDED_MODULES
View File

@ -1 +1,2 @@
Ezfio_files Nuclei Utils
Ezfio_files Nuclei Output Utils

49
src/AOs/README.rst
View File

@ -21,7 +21,7 @@ Assumptions
===========
.. Do not edit this section. It was auto-generated from the
.. ASSUMPTIONS.rst file.
.. NEEDED_MODULES file.
* The atomic orbitals are normalized:
@ -41,5 +41,52 @@ Needed Modules
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`ao_coef <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L21>`_
Coefficients, exponents and powers of x,y and z
ao_coef(i,j) = coefficient of the jth primitive on the ith ao
`ao_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L136>`_
Transposed ao_coef and ao_expo
`ao_expo <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L20>`_
Coefficients, exponents and powers of x,y and z
ao_coef(i,j) = coefficient of the jth primitive on the ith ao
`ao_expo_transp <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L137>`_
Transposed ao_coef and ao_expo
`ao_nucl <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L186>`_
Index of the nuclei on which the ao is centered
`ao_num <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L1>`_
Number of atomic orbitals
`ao_num_align <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L2>`_
Number of atomic orbitals
`ao_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L96>`_
matrix of the overlap for tha AOs
S(i,j) = overlap between the ith and the jth atomic basis function
`ao_power <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L19>`_
Coefficients, exponents and powers of x,y and z
ao_coef(i,j) = coefficient of the jth primitive on the ith ao
`ao_prim_num <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L154>`_
Number of primitives per atomic orbital
`ao_prim_num_max <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L176>`_
None
`ao_prim_num_max_align <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L177>`_
None

165
src/BiInts/README.rst
View File

@ -26,3 +26,168 @@ Needed Modules
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`ao_bielec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L1>`_
integral of the AO basis <ik|jl> or (ij|kl)
i(r1) j(r1) 1/r12 k(r2) l(r2)
`ao_bielec_integral_schwartz <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L326>`_
Needed to compuet Schwartz inequalities
`ao_bielec_integrals_in_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L189>`_
Map of Atomic integrals
i(r1) j(r2) 1/r12 k(r1) l(r2)
`compute_ao_bielec_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L148>`_
Compute AO 1/r12 integrals for all i and fixed j,k,l
`eri <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L487>`_
ATOMIC PRIMTIVE bielectronic integral between the 4 primitives ::
primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2)
primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2)
primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2)
primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2)
`general_primitive_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L352>`_
Computes the integral <pq|rs> where p,q,r,s are Gaussian primitives
`give_polynom_mult_center_x <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L632>`_
subroutine that returns the explicit polynom in term of the "t"
variable of the following polynomw :
I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q)
`i_x1_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L576>`_
recursive function involved in the bielectronic integral
`i_x1_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L695>`_
recursive function involved in the bielectronic integral
`i_x1_pol_mult_a1 <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L815>`_
recursive function involved in the bielectronic integral
`i_x1_pol_mult_a2 <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L869>`_
recursive function involved in the bielectronic integral
`i_x1_pol_mult_recurs <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L729>`_
recursive function involved in the bielectronic integral
`i_x2_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L599>`_
recursive function involved in the bielectronic integral
`i_x2_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L931>`_
recursive function involved in the bielectronic integral
`integrale_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L531>`_
calculate the integral of the polynom ::
I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q)
between ( 0 ; 1)
`n_pt_sup <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L618>`_
Returns the upper boundary of the degree of the polynom involved in the
bielctronic integral :
Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z)
`gauleg <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/gauss_legendre.irp.f#L20>`_
Gauss-Legendre
`gauleg_t2 <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/gauss_legendre.irp.f#L1>`_
t_w(i,1,k) = w(i)
t_w(i,2,k) = t(i)
`gauleg_w <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/gauss_legendre.irp.f#L2>`_
t_w(i,1,k) = w(i)
t_w(i,2,k) = t(i)
`ao_integrals_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L6>`_
AO integrals
`bielec_integrals_index <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L17>`_
None
`clear_ao_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L128>`_
Frees the memory of the AO map
`clear_mo_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L243>`_
Frees the memory of the MO map
`get_ao_bielec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L33>`_
Gets one AO bi-electronic integral from the AO map
`get_ao_bielec_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L51>`_
Gets multiple AO bi-electronic integral from the AO map .
All i are retrieved for j,k,l fixed.
`get_ao_bielec_integrals_non_zero <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L84>`_
Gets multiple AO bi-electronic integral from the AO map .
All non-zero i are retrieved for j,k,l fixed.
`get_ao_map_size <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L120>`_
Returns the number of elements in the AO map
`get_mo_bielec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L184>`_
Returns one integral <ij|kl> in the MO basis
`get_mo_bielec_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L213>`_
Returns multiple integrals <ij|kl> in the MO basis, all
i for j,k,l fixed.
`get_mo_map_size <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L235>`_
Return the number of elements in the MO map
`insert_into_ao_integrals_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L153>`_
Create new entry into AO map
`insert_into_mo_integrals_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L168>`_
Create new entry into MO map, or accumulate in an existing entry
`mo_bielec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L201>`_
Returns one integral <ij|kl> in the MO basis
`mo_integrals_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L142>`_
MO integrals
`add_integrals_to_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L40>`_
Adds integrals to tha MO map according to some bitmask
`mo_bielec_integral_jj <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L296>`_
Transform Bi-electronic integrals <ij|ij> and <ij|ji>
`mo_bielec_integral_jj_anti <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L298>`_
Transform Bi-electronic integrals <ij|ij> and <ij|ji>
`mo_bielec_integral_jj_exchange <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L297>`_
Transform Bi-electronic integrals <ij|ij> and <ij|ji>
`mo_bielec_integrals_in_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L21>`_
If True, the map of MO bielectronic integrals is provided
`mo_bielec_integrals_index <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L1>`_
Computes an unique index for i,j,k,l integrals
`ao_integrals_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L73>`_
If <pq|rs> < ao_integrals_threshold, <pq|rs> = 0
`do_direct_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L111>`_
If True, compute integrals on the fly
`mo_integrals_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L92>`_
If <ij|kl> < mo_integrals_threshold, <ij|kl> = 0
`read_ao_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L55>`_
If true, read AO integrals in EZFIO
`read_mo_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L37>`_
If true, read MO integrals in EZFIO
`write_ao_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L19>`_
If true, write AO integrals in EZFIO
`write_mo_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L1>`_
If true, write MO integrals in EZFIO

13
src/BiInts/ao_bi_integrals.irp.f
View File

@ -160,20 +160,17 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value)
thresh = ao_integrals_threshold
integer :: n_centers, i
integer*1 :: center_count(nucl_num)
PROVIDE gauleg_t2 ao_nucl all_utils
if (ao_overlap_abs(j,l) < thresh) then
buffer_value = 0.
buffer_value = 0._integral_kind
return
endif
center_count = 0
do i = 1, ao_num
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
buffer_value(i) = 0.
buffer_value(i) = 0._integral_kind
cycle
endif
!DIR$ FORCEINLINE
@ -193,7 +190,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
implicit none
use map_module
BEGIN_DOC
! Map of Atomic integrals :
! Map of Atomic integrals
! i(r1) j(r2) 1/r12 k(r1) l(r2)
END_DOC
@ -211,7 +208,6 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
integer :: n_integrals, n_centers
integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax
integer*1 :: center_count(nucl_num)
PROVIDE gauleg_t2 ao_integrals_map all_utils
integral = ao_bielec_integral(1,1,1,1)
@ -243,7 +239,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
call cpu_time(cpu_1)
!$OMP PARALLEL PRIVATE(i,j,k,l,kk, &
!$OMP integral,buffer_i,buffer_value,n_integrals, &
!$OMP cpu_2,wall_2,i1,j1,center_count) &
!$OMP cpu_2,wall_2,i1,j1) &
!$OMP DEFAULT(NONE) &
!$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, &
!$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, &
@ -252,7 +248,6 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
allocate(buffer_i(size_buffer))
allocate(buffer_value(size_buffer))
n_integrals = 0
center_count = 0
!$OMP DO SCHEDULE(dynamic)
do kk=1,lmax

1
src/BiInts/bi_integrals.ezfio_config
View File

@ -5,4 +5,5 @@ bielec_integrals
write_mo_integrals logical
threshold_ao double precision
threshold_mo double precision
direct logical

2
src/BiInts/map_integrals.irp.f
View File

@ -330,5 +330,5 @@ end
SUBST [ ao_integrals_map, ao_integrals, ao_num , get_ao_bielec_integral ]
ao_integrals_map ; ao_integrals ; ao_num ; get_ao_bielec_integral ;;
mo_integrals_map ; mo_integrals ; n_act ; get_mo_bielec_integral ;;
mo_integrals_map ; mo_integrals ; mo_tot_num ; get_mo_bielec_integral ;;
END_TEMPLATE

183
src/BiInts/mo_bi_integrals.irp.f
View File

@ -72,7 +72,6 @@ subroutine add_integrals_to_map(mask_ijkl)
PROVIDE N_int ao_bielec_integrals_in_map ao_integrals_map mo_coef mo_coef_transp
!Get list of MOs for i,j,k and l
!-------------------------------
@ -82,9 +81,6 @@ subroutine add_integrals_to_map(mask_ijkl)
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
l1_global = 0
size_buffer = min(ao_num*ao_num*ao_num,16000000)
write(output_BiInts,*) 'Providing the molecular integrals '
@ -294,63 +290,138 @@ subroutine add_integrals_to_map(mask_ijkl)
end
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num)]
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij>
! Transform Bi-electronic integrals <ij|ij> and <ij|ji>
END_DOC
integer :: i,j
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
do j= 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i= 1, mo_tot_num
mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map)
enddo
! Padding
do i= mo_tot_num+1, mo_tot_num_align
mo_bielec_integral_jj(i,j) = 0.d0
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij> - <ij|ji>
END_DOC
integer :: i,j
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
do j = 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i = 1,mo_tot_num
mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map)
enddo
! Padding
do i= mo_tot_num+1, mo_tot_num_align
mo_bielec_integral_jj_exchange(i,j) = 0.d0
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij> - <ij|ji>
END_DOC
integer :: i,j
PROVIDE mo_bielec_integrals_in_map
do j = 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i = 1,mo_tot_num_align
mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j)
integer :: i,j,p,q,r,s
double precision :: c
real(integral_kind) :: integral
integer :: n, pp
real(integral_kind), allocatable :: int_value(:)
integer, allocatable :: int_idx(:)
double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:)
PROVIDE ao_integrals_threshold
if (.not.do_direct_integrals) then
PROVIDE ao_bielec_integrals_in_map
endif
mo_bielec_integral_jj = 0.d0
mo_bielec_integral_jj_exchange = 0.d0
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, &
!$OMP iqrs, iqsr,iqri,iqis) &
!$OMP SHARED(mo_tot_num,mo_coef_transp,mo_tot_num_align,ao_num,&
!$OMP ao_integrals_threshold,do_direct_integrals) &
!$OMP REDUCTION(+:mo_bielec_integral_jj,mo_bielec_integral_jj_exchange)
allocate( int_value(ao_num), int_idx(ao_num), &
iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),&
iqsr(mo_tot_num_align,ao_num) )
!$OMP DO SCHEDULE (guided)
do s=1,ao_num
do q=1,ao_num
do j=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,j) = 0.d0
iqsr(i,j) = 0.d0
enddo
enddo
if (do_direct_integrals) then
double precision :: ao_bielec_integral
do r=1,ao_num
call compute_ao_bielec_integrals(q,r,s,ao_num,int_value)
do p=1,ao_num
integral = int_value(p)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
call compute_ao_bielec_integrals(q,s,r,ao_num,int_value)
do p=1,ao_num
integral = int_value(p)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqsr(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
enddo
else
do r=1,ao_num
call get_ao_bielec_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n)
do pp=1,n
p = int_idx(pp)
integral = int_value(pp)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
call get_ao_bielec_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n)
do pp=1,n
p = int_idx(pp)
integral = int_value(pp)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqsr(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
enddo
endif
iqis = 0.d0
iqri = 0.d0
do r=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqis(i) += mo_coef_transp(i,r) * iqrs(i,r)
iqri(i) += mo_coef_transp(i,r) * iqsr(i,r)
enddo
enddo
do i=1,mo_tot_num
!DIR$ VECTOR ALIGNED
do j=1,mo_tot_num
c = mo_coef_transp(j,q)*mo_coef_transp(j,s)
mo_bielec_integral_jj(j,i) += c * iqis(i)
mo_bielec_integral_jj_exchange(j,i) += c * iqri(i)
enddo
enddo
enddo
enddo
!$OMP END DO NOWAIT
deallocate(iqrs,iqsr,int_value,int_idx)
!$OMP END PARALLEL
mo_bielec_integral_jj_anti = mo_bielec_integral_jj - mo_bielec_integral_jj_exchange
END_PROVIDER

18
src/BiInts/options.irp.f
View File

@ -107,3 +107,21 @@ BEGIN_PROVIDER [ double precision, mo_integrals_threshold ]
END_PROVIDER
BEGIN_PROVIDER [ logical, do_direct_integrals ]
implicit none
BEGIN_DOC
! If True, compute integrals on the fly
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_direct(has)
if (has) then
call ezfio_get_bielec_integrals_direct(do_direct_integrals)
else
do_direct_integrals = .False.
call ezfio_set_bielec_integrals_direct(do_direct_integrals)
endif
END_PROVIDER

43
src/Bitmask/README.rst
View File

@ -22,7 +22,7 @@ Assumptions
===========
.. Do not edit this section. It was auto-generated from the
.. ASSUMPTIONS.rst file.
.. NEEDED_MODULES file.
``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent:
@ -48,3 +48,44 @@ Needed Modules
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`full_ijkl_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L12>`_
Bitmask to include all possible MOs
`hf_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L32>`_
Hartree Fock bit mask
`n_int <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L3>`_
Number of 64-bit integers needed to represent determinants as binary strings
`ref_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L50>`_
Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask
`bitstring_to_hexa <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L95>`_
Transform a bit string to a string in hexadecimal format for printing
`bitstring_to_list <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L1>`_
Gives the inidices(+1) of the bits set to 1 in the bit string
`bitstring_to_str <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L62>`_
Transform a bit string to a string for printing
`debug_det <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L117>`_
None
`list_to_bitstring <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L29>`_
return the physical string "string(N_int,2)" from the array of occupations "list(N_int*bit_kind_size,2)
list
<== ipos ==>
|
v
string :|------------------------|-------------------------|------------------------|
<==== bit_kind_size ====> <==== bit_kind_size ====> <==== bit_kind_size ====>
{ iint } { iint } { iint }

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

@ -55,3 +55,4 @@ BEGIN_PROVIDER [ integer(bit_kind), ref_bitmask, (N_int,2)]
ref_bitmask = HF_bitmask
END_PROVIDER

16
src/Bitmask/bitmasks_routines.irp.f
View File

@ -114,4 +114,18 @@ subroutine bitstring_to_hexa( output, string, Nint )
enddo
end
subroutine debug_det(string,Nint)
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: string(Nint,2)
character*(512) :: output(2)
call bitstring_to_hexa( output(1), string(1,1), Nint )
call bitstring_to_hexa( output(2), string(1,2), Nint )
print *, trim(output(1)) , '|', trim(output(2))
call bitstring_to_str( output(1), string(1,1), N_int )
call bitstring_to_str( output(2), string(1,2), N_int )
print *, trim(output(1))
print *, trim(output(2))
end

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

67
src/CISD/H_apply.irp.f Normal file
View File

@ -0,0 +1,67 @@
BEGIN_SHELL [ /bin/bash ]
./h_apply.py
END_SHELL
subroutine fill_H_apply_buffer_cisd(n_selected,det_buffer,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Fill the H_apply buffer with determiants for CISD
END_DOC
integer, intent(in) :: n_selected, Nint
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k
integer :: new_size
new_size = H_apply_buffer_N_det + n_selected
if (new_size > h_apply_buffer_size) then
call resize_h_apply_buffer_det(max(h_apply_buffer_size*2,new_size))
endif
do i=1,H_apply_buffer_N_det
ASSERT (sum(popcnt(h_apply_buffer_det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(h_apply_buffer_det(:,2,i))) == elec_beta_num)
enddo
do i=1,n_selected
do j=1,N_int
h_apply_buffer_det(j,1,i+H_apply_buffer_N_det) = det_buffer(j,1,i)
h_apply_buffer_det(j,2,i+H_apply_buffer_N_det) = det_buffer(j,2,i)
enddo
ASSERT (sum(popcnt(h_apply_buffer_det(:,1,i+H_apply_buffer_N_det)) )== elec_alpha_num)
ASSERT (sum(popcnt(h_apply_buffer_det(:,2,i+H_apply_buffer_N_det))) == elec_beta_num)
H_apply_buffer_coef(i,:) = 0.d0
enddo
H_apply_buffer_N_det = new_size
do i=1,H_apply_buffer_N_det
ASSERT (sum(popcnt(h_apply_buffer_det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(h_apply_buffer_det(:,2,i))) == elec_beta_num)
enddo
SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det
end
subroutine H_apply_cisd
implicit none
BEGIN_DOC
! Calls H_apply on the HF determinant and selects all connected single and double
! excitations (of the same symmetry).
END_DOC
integer(bit_kind) :: hole_mask(N_int,2)
integer(bit_kind) :: particle_mask(N_int,2)
hole_mask(:,1) = HF_bitmask(:,1)
hole_mask(:,2) = HF_bitmask(:,2)
particle_mask(:,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1))
particle_mask(:,2) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2))
call H_apply_cisd_monoexc(HF_bitmask, &
hole_mask, particle_mask)
call H_apply_cisd_diexc(HF_bitmask, &
hole_mask, particle_mask, &
hole_mask, particle_mask )
call copy_H_apply_buffer_to_wf
end

8
src/CISD/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/CISD/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils

36
src/CISD/README.rst Normal file
View File

@ -0,0 +1,36 @@
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`fill_h_apply_buffer_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L6>`_
Fill the H_apply buffer with determiants for CISD
`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L43>`_
Calls H_apply on the HF determinant and selects all connected single and double
excitations (of the same symmetry).
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd.irp.f#L1>`_
None

19
src/CISD/cisd.irp.f Normal file
View File

@ -0,0 +1,19 @@
program cisd
implicit none
integer :: i
call H_apply_cisd
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
allocate(eigvalues(n_det),eigvectors(n_det,n_det))
print *, 'N_det = ', N_det
call lapack_diag(eigvalues,eigvectors,H_matrix_all_dets,n_det,n_det)
! print *, H_matrix_all_dets
print *, '---'
print *, 'HF:', HF_energy
print *, '---'
do i = 1,3
print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion
enddo
! print *, eigvectors(:,1)
deallocate(eigvalues,eigvectors)
end

11
src/CISD/h_apply.py Executable file
View File

@ -0,0 +1,11 @@
#!/usr/bin/env python
import generate_h_apply
# H_apply
s = generate_h_apply.new_dict(openmp=True)
s["subroutine"] = "H_apply_cisd"
s["keys_work"] = "call fill_H_apply_buffer_cisd(key_idx,keys_out,N_int)"
generate_h_apply.create_h_apply(s)

33
src/CISD/tests/Makefile Normal file
View File

@ -0,0 +1,33 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90+= -I tests
REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f))
.PHONY: clean executables serial_tests parallel_tests
all: clean executables serial_tests parallel_tests
parallel_tests: $(REF_FILES)
@echo ; echo " ---- Running parallel tests ----" ; echo
@OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py
serial_tests: $(REF_FILES)
@echo ; echo " ---- Running serial tests ----" ; echo
@OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py
executables: $(wildcard *.irp.f) veryclean
$(MAKE) -C ..
%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables
$(QPACKAGE_ROOT)/scripts/create_test_ref.sh $*
clean:
$(MAKE) -C .. clean
veryclean:
$(MAKE) -C .. veryclean

31
src/DensityMatrix/README.rst Normal file
View File

@ -0,0 +1,31 @@
====================
DensityMatrix Module
====================
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_

7
src/Dets/ASSUMPTIONS.rst Normal file
View File

@ -0,0 +1,7 @@
* The MOs are orthonormal
* All the determinants have the same number of electrons
* The determinants are orthonormal
* The number of generator determinants <= the number of determinants
* All the determinants in the H_apply buffer are supposed to be different from the
wave function determinants
* All the determinants in the H_apply buffer are supposed to be unique

157
src/Dets/H_apply.irp.f Normal file
View File

@ -0,0 +1,157 @@
use bitmasks
BEGIN_PROVIDER [ double precision, H_apply_threshold ]
implicit none
BEGIN_DOC
! Theshold on | <Di|H|Dj> |
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_determinants_H_apply_threshold(has)
if (has) then
call ezfio_get_determinants_H_apply_threshold(H_apply_threshold)
else
H_apply_threshold = 1.d-10
call ezfio_set_determinants_H_apply_threshold(H_apply_threshold)
endif
call write_time(output_Dets)
call write_double(output_Dets, H_apply_threshold, &
'H_apply_threshold')
END_PROVIDER
BEGIN_PROVIDER [ integer*8, H_apply_buffer_size ]
implicit none
BEGIN_DOC
! Size of the H_apply buffer.
END_DOC
H_apply_buffer_size = 1000
END_PROVIDER
subroutine resize_H_apply_buffer_det(new_size)
implicit none
integer, intent(in) :: new_size
integer(bit_kind), allocatable :: buffer_det(:,:,:)
double precision, allocatable :: buffer_coef(:,:)
integer :: i,j,k
integer :: Ndet
ASSERT (new_size > 0)
allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states) )
do i=1,min(new_size,H_apply_buffer_N_det)
do k=1,N_int
buffer_det(k,1,i) = H_apply_buffer_det(k,1,i)
buffer_det(k,2,i) = H_apply_buffer_det(k,2,i)
enddo
ASSERT (sum(popcnt(H_apply_buffer_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer_det(:,2,i))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,min(new_size,H_apply_buffer_N_det)
buffer_coef(i,k) = H_apply_buffer_coef(i,k)
enddo
enddo
H_apply_buffer_size = new_size
Ndet = min(new_size,H_apply_buffer_N_det)
TOUCH H_apply_buffer_size
H_apply_buffer_N_det = Ndet
do i=1,H_apply_buffer_N_det
do k=1,N_int
H_apply_buffer_det(k,1,i) = buffer_det(k,1,i)
H_apply_buffer_det(k,2,i) = buffer_det(k,2,i)
enddo
ASSERT (sum(popcnt(H_apply_buffer_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer_det(:,2,i))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,H_apply_buffer_N_det
H_apply_buffer_coef(i,k) = buffer_coef(i,k)
enddo
enddo
deallocate (buffer_det, buffer_coef)
SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det
end
BEGIN_PROVIDER [ integer(bit_kind), H_apply_buffer_det,(N_int,2,H_apply_buffer_size) ]
&BEGIN_PROVIDER [ double precision, H_apply_buffer_coef,(H_apply_buffer_size,N_states) ]
&BEGIN_PROVIDER [ integer, H_apply_buffer_N_det ]
implicit none
BEGIN_DOC
! Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines.
END_DOC
H_apply_buffer_N_det = 0
END_PROVIDER
subroutine copy_H_apply_buffer_to_wf
implicit none
integer(bit_kind), allocatable :: buffer_det(:,:,:)
double precision, allocatable :: buffer_coef(:,:)
integer :: i,j,k
integer :: N_det_old
ASSERT (N_int > 0)
ASSERT (N_det > 0)
allocate ( buffer_det(N_int,2,N_det), buffer_coef(N_det,N_states) )
do i=1,N_det
do k=1,N_int
ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num)
buffer_det(k,1,i) = psi_det(k,1,i)
buffer_det(k,2,i) = psi_det(k,2,i)
enddo
enddo
do k=1,N_states
do i=1,N_det
buffer_coef(i,k) = psi_coef(i,k)
enddo
enddo
N_det_old = N_det
N_det = N_det + H_apply_buffer_N_det
TOUCH N_det
do i=1,N_det_old
do k=1,N_int
psi_det(k,1,i) = buffer_det(k,1,i)
psi_det(k,2,i) = buffer_det(k,2,i)
enddo
ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num )
enddo
do i=1,H_apply_buffer_N_det
do k=1,N_int
psi_det(k,1,i+N_det_old) = H_apply_buffer_det(k,1,i)
psi_det(k,2,i+N_det_old) = H_apply_buffer_det(k,2,i)
enddo
ASSERT (sum(popcnt(psi_det(:,1,i+N_det_old))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i+N_det_old))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,N_det_old
psi_coef(i,k) = buffer_coef(i,k)
enddo
do i=1,H_apply_buffer_N_det
psi_coef(i+N_det_old,k) = H_apply_buffer_coef(i,k)
enddo
enddo
SOFT_TOUCH psi_det psi_coef
end

356
src/Dets/H_apply_template.f Normal file
View File

@ -0,0 +1,356 @@
subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parameters )
use omp_lib
use bitmasks
implicit none
$declarations
integer(omp_lock_kind) :: lck
integer(bit_kind),intent(in) :: key_in(N_int,2)
integer(bit_kind),allocatable :: keys_out(:,:,:)
double precision, allocatable :: hij_tab(:)
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(bit_kind) :: hole_save(N_int,2)
integer(bit_kind) :: key(N_int,2),hole(N_int,2), particle(N_int,2)
integer(bit_kind) :: hole_tmp(N_int,2), particle_tmp(N_int,2)
integer :: ii,i,jj,j,k,ispin,l
integer :: occ_particle(N_int*bit_kind_size,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_particle_tmp(N_int*bit_kind_size,2)
integer :: occ_hole_tmp(N_int*bit_kind_size,2)
integer :: kk,pp,other_spin,key_idx
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
integer,parameter :: size_max = $size_max
double precision :: hij_elec, mo_bielec_integral, thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem, E_ref
PROVIDE mo_integrals_map
PROVIDE mo_bielec_integrals_in_map
$set_i_H_j_threshold
$omp_init_lock
E_ref = diag_H_mat_elem(key_in,N_int)
$initialization
$omp_parallel
allocate (keys_out(N_int,2,size_max),hij_tab(size_max))
!print*,'key_in !!'
!call print_key(key_in)
!print*,'hole_1, particl_1'
!call print_key(hole_1)
!call print_key(particl_1)
!print*,'hole_2, particl_2'
!call print_key(hole_2)
!call print_key(particl_2)
!!!! First couple hole particle
do j = 1, N_int
hole(j,1) = iand(hole_1(j,1),key_in(j,1))
hole(j,2) = iand(hole_1(j,2),key_in(j,2))
particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1))
particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2))
enddo
call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int)
call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
call bitstring_to_list(hole(1,1),occ_hole(1,1),N_elec_in_key_hole_1(1),N_int)
call bitstring_to_list(hole(1,2),occ_hole(1,2),N_elec_in_key_hole_1(2),N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2))
do ispin=1,2
i=0
do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole
i_a = occ_hole(ii,ispin)
ASSERT (i_a > 0)
ASSERT (i_a <= mo_tot_num)
do jj=1,N_elec_in_key_part_1(ispin) !particle
j_a = occ_particle(jj,ispin)
ASSERT (j_a > 0)
ASSERT (j_a <= mo_tot_num)
i += 1
ia_ja_pairs(1,i,ispin) = i_a
ia_ja_pairs(2,i,ispin) = j_a
enddo
enddo
ia_ja_pairs(1,0,ispin) = i
enddo
key_idx = 0
integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b
integer(bit_kind) :: test(N_int,2)
double precision :: accu
accu = 0.d0
hij_elec = 0.d0
do ispin=1,2
other_spin = iand(ispin,1)+1
$omp_do
do ii=1,ia_ja_pairs(1,0,ispin)
i_a = ia_ja_pairs(1,ii,ispin)
ASSERT (i_a > 0)
ASSERT (i_a <= mo_tot_num)
j_a = ia_ja_pairs(2,ii,ispin)
ASSERT (j_a > 0)
ASSERT (j_a <= mo_tot_num)
hole = key_in
k = ishft(i_a-1,-bit_kind_shift)+1
j = i_a-ishft(k-1,bit_kind_shift)-1
hole(k,ispin) = ibclr(hole(k,ispin),j)
k_a = ishft(j_a-1,-bit_kind_shift)+1
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
!!!! Second couple hole particle
do j = 1, N_int
hole_tmp(j,1) = iand(hole_2(j,1),hole(j,1))
hole_tmp(j,2) = iand(hole_2(j,2),hole(j,2))
particle_tmp(j,1) = iand(xor(particl_2(j,1),hole(j,1)),particl_2(j,1))
particle_tmp(j,2) = iand(xor(particl_2(j,2),hole(j,2)),particl_2(j,2))
enddo
call bitstring_to_list(particle_tmp(1,1),occ_particle_tmp(1,1),N_elec_in_key_part_2(1),N_int)
call bitstring_to_list(particle_tmp(1,2),occ_particle_tmp(1,2),N_elec_in_key_part_2(2),N_int)
call bitstring_to_list(hole_tmp (1,1),occ_hole_tmp (1,1),N_elec_in_key_hole_2(1),N_int)
call bitstring_to_list(hole_tmp (1,2),occ_hole_tmp (1,2),N_elec_in_key_hole_2(2),N_int)
! hole = a^(+)_j_a(ispin) a_i_a(ispin)|key_in> : mono exc :: orb(i_a,ispin) --> orb(j_a,ispin)
hole_save = hole
if (ispin == 1) then
integer :: jjj
do kk = 1,N_elec_in_key_hole_2(other_spin)
hole = hole_save
i_b = occ_hole_tmp(kk,other_spin)
ASSERT (i_b > 0)
ASSERT (i_b <= mo_tot_num)
k = ishft(i_b-1,-bit_kind_shift)+1
j = i_b-ishft(k-1,bit_kind_shift)-1
hole(k,other_spin) = ibclr(hole(k,other_spin),j)
do jjj=1,N_elec_in_key_part_2(other_spin) ! particule
j_b = occ_particle_tmp(jjj,other_spin)
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if(dabs( mo_bielec_integral(j_a,j_b,i_a,i_b))<thresh)cycle
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,other_spin) = ibset(key(k,other_spin),l)
call i_H_j(key,key_in,N_int,hij_elec)
if(dabs(hij_elec)>=thresh) then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$omp_set_lock
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
enddo
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
enddo
endif
! does all the mono excitations of the same spin
do kk = 1,N_elec_in_key_hole_2(ispin)
i_b = occ_hole_tmp(kk,ispin)
ASSERT (i_b > 0)
ASSERT (i_b <= mo_tot_num)
if (i_b <= i_a.or.i_b == j_a) cycle
hole = hole_save
k = ishft(i_b-1,-bit_kind_shift)+1
j = i_b-ishft(k-1,bit_kind_shift)-1
hole(k,ispin) = ibclr(hole(k,ispin),j)
do jjj=1,N_elec_in_key_part_2(ispin)
j_b = occ_particle_tmp(jjj,ispin)
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if (j_b <= j_a) cycle
if(dabs( mo_bielec_integral(j_a,j_b,i_b,i_a))<thresh)cycle
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,ispin) = ibset(key(k,ispin),l)
!! a^((+)_j_b(ispin) a_i_b(ispin) : mono exc :: orb(i_b,ispin) --> orb(j_b,ispin)
call i_H_j(key,key_in,N_int,hij_elec)
if(dabs(hij_elec)>=thresh) then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$omp_set_lock
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
enddo
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
enddo! kk
enddo ! ii
$omp_enddo
enddo ! ispin
$omp_set_lock
$keys_work
$omp_unset_lock
deallocate (keys_out,hij_tab,ia_ja_pairs)
$omp_end_parallel
$omp_destroy_lock
$finalization
end
subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
use omp_lib
use bitmasks
implicit none
$declarations
integer(omp_lock_kind) :: lck
integer(bit_kind),intent(in) :: key_in(N_int,2)
integer(bit_kind),allocatable :: keys_out(:,:,:)
double precision, allocatable :: hij_tab(:)
integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind) :: hole_2(N_int,2), particl_2(N_int,2)
integer(bit_kind) :: hole_save(N_int,2)
integer(bit_kind) :: key(N_int,2),hole(N_int,2), particle(N_int,2)
integer(bit_kind) :: hole_tmp(N_int,2), particle_tmp(N_int,2)
integer :: ii,i,jj,j,k,ispin,l
integer :: occ_particle(N_int*bit_kind_size,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_particle_tmp(N_int*bit_kind_size,2)
integer :: occ_hole_tmp(N_int*bit_kind_size,2)
integer :: kk,pp,other_spin,key_idx
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
integer,parameter :: size_max = $size_max
double precision :: hij_elec, thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem, E_ref
PROVIDE mo_integrals_map
PROVIDE mo_bielec_integrals_in_map
$set_i_H_j_threshold
$omp_init_lock
E_ref = diag_H_mat_elem(key_in,N_int)
$initialization
$omp_parallel
allocate (keys_out(N_int,2,size_max),hij_tab(size_max))
!!!! First couple hole particle
do j = 1, N_int
hole(j,1) = iand(hole_1(j,1),key_in(j,1))
hole(j,2) = iand(hole_1(j,2),key_in(j,2))
particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1))
particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2))
enddo
call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int)
call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
call bitstring_to_list(hole (1,1),occ_hole (1,1),N_elec_in_key_hole_1(1),N_int)
call bitstring_to_list(hole (1,2),occ_hole (1,2),N_elec_in_key_hole_1(2),N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2))
do ispin=1,2
i=0
do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole
i_a = occ_hole(ii,ispin)
do jj=1,N_elec_in_key_part_1(ispin) !particule
j_a = occ_particle(jj,ispin)
i += 1
ia_ja_pairs(1,i,ispin) = i_a
ia_ja_pairs(2,i,ispin) = j_a
enddo
enddo
ia_ja_pairs(1,0,ispin) = i
enddo
key_idx = 0
integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b
integer(bit_kind) :: test(N_int,2)
double precision :: accu
accu = 0.d0
hij_elec = 0.d0
do ispin=1,2
other_spin = iand(ispin,1)+1
$omp_do
do ii=1,ia_ja_pairs(1,0,ispin)
i_a = ia_ja_pairs(1,ii,ispin)
j_a = ia_ja_pairs(2,ii,ispin)
hole = key_in
k = ishft(i_a-1,-bit_kind_shift)+1
j = i_a-ishft(k-1,bit_kind_shift)-1
hole(k,ispin) = ibclr(hole(k,ispin),j)
k_a = ishft(j_a-1,-bit_kind_shift)+1
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
call i_H_j(hole,key_in,N_int,hij_elec)
if(dabs(hij_elec) .ge. thresh)then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)
keys_out(k,2,key_idx) = hole(k,2)
enddo
hij_tab(key_idx) = hij_elec
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
if (key_idx == size_max) then
$omp_set_lock
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
enddo ! ii
$omp_enddo
enddo ! ispin
$omp_set_lock
$keys_work
$omp_unset_lock
deallocate (keys_out,hij_tab,ia_ja_pairs)
$omp_end_parallel
$omp_destroy_lock
$finalization
end

1
src/Dets/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils

159
src/Dets/README.rst Normal file
View File

@ -0,0 +1,159 @@
===========
Dets Module
===========
This module contains the determinants of the CI wave function.
H is applied on the list of generator determinants. Selected determinants
are added into the *H_apply buffer*. Then the new wave function is
constructred as the concatenation of the odl wave function and
some determinants of the H_apply buffer. Generator determinants are built
as a subset of the determinants of the wave function.
Assumptions
===========
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* The MOs are orthonormal
* All the determinants have the same number of electrons
* The determinants are orthonormal
* The number of generator determinants <= the number of determinants
* All the determinants in the H_apply buffer are supposed to be different from the
wave function determinants
* All the determinants in the H_apply buffer are supposed to be unique
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L93>`_
None
`h_apply_buffer_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L82>`_
Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines.
`h_apply_buffer_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L81>`_
Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines.
`h_apply_buffer_n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L83>`_
Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines.
`h_apply_buffer_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L22>`_
Size of the H_apply buffer.
`h_apply_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L3>`_
Theshold on | <Di|H|Dj> |
`resize_h_apply_buffer_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L31>`_
None
`n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L11>`_
Number of determinants in the wave function
`n_det_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L47>`_
Number of generator determinants in the wave function
`n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L3>`_
Number of states to consider
`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
The wave function. Initialized with Hartree-Fock
`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L19>`_
The wave function. Initialized with Hartree-Fock
`psi_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L55>`_
Determinants on which H is applied
`double_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L40>`_
double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1
double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1
double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2
double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2
for a given couple of hole/particle excitations i.
`n_double_exc_bitmasks <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L31>`_
Number of double excitation bitmasks
`n_single_exc_bitmasks <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L8>`_
Number of single excitation bitmasks
`single_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L17>`_
single_exc_bitmask(:,1,i) is the bitmask for holes
single_exc_bitmask(:,2,i) is the bitmask for particles
for a given couple of hole/particle excitations i.
`get_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L1>`_
Returns <S^2>
`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L842>`_
Needed for diag_H_mat_elem
`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L887>`_
Needed for diag_H_mat_elem
`decode_exc <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L76>`_
Decodes the exc arrays returned by get_excitation.
h1,h2 : Holes
p1,p2 : Particles
s1,s2 : Spins (1:alpha, 2:beta)
degree : Degree of excitation
`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L779>`_
Computes <i|H|i>
`filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L602>`_
Filters out the determinants that are not connected by H
`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L687>`_
None
`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L140>`_
Returns the two excitation operators between two doubly excited determinants and the phase
`get_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L30>`_
Returns the excitation operators between two determinants and the phase
`get_excitation_degree <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L1>`_
Returns the excitation degree between two determinants
`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L518>`_
Applies get_excitation_degree to an array of determinants
`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L273>`_
Returns the excitation operator between two singly excited determinants and the phase
`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L935>`_
Returns a list of occupation numbers from a bitstring
`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L354>`_
Returns <i|H|j> where i and j are determinants
`i_h_psim <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L490>`_
None
`h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/utils.irp.f#L1>`_
H matrix on the basis of the slater deter;inants defined by psi_det

9
src/Dets/determinants.ezfio_config Normal file
View File

@ -0,0 +1,9 @@
determinants
N_int integer
bit_kind integer
n_dets integer
n_states integer
psi_coef double precision (determinants_n_dets,determinants_n_states)
psi_det integer (determinants_N_int*determinants_bit_kind/4,2,determinants_n_dets)
H_apply_threshold double precision

69
src/Dets/determinants.irp.f Normal file
View File

@ -0,0 +1,69 @@
use bitmasks
BEGIN_PROVIDER [ integer, N_states ]
implicit none
BEGIN_DOC
! Number of states to consider
END_DOC
N_states = 1
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det ]
implicit none
BEGIN_DOC
! Number of determinants in the wave function
END_DOC
N_det = max(1,N_states)
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,N_det) ]
&BEGIN_PROVIDER [ double precision, psi_coef, (N_det,N_states) ]
implicit none
BEGIN_DOC
! The wave function. Initialized with Hartree-Fock
END_DOC
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
psi_det = 0_bit_kind
psi_coef = 0.d0
integer :: i
do i=1,N_int
psi_det(i,1,1) = HF_bitmask(i,1)
psi_det(i,2,1) = HF_bitmask(i,2)
enddo
do i=1,N_states
psi_coef(i,i) = 1.d0
enddo
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_generators ]
implicit none
BEGIN_DOC
! Number of generator determinants in the wave function
END_DOC
N_det_generators = N_det
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,N_det) ]
implicit none
BEGIN_DOC
! Determinants on which H is applied
END_DOC
psi_generators = 0_bit_kind
integer :: i
do i=1,N_int
psi_generators(i,1,1) = psi_det(i,1,1)
psi_generators(i,2,1) = psi_det(i,1,1)
enddo
END_PROVIDER

57
src/Dets/determinants_bitmasks.irp.f Normal file
View File

@ -0,0 +1,57 @@
use bitmasks
integer, parameter :: hole_ = 1
integer, parameter :: particle_ = 2
integer, parameter :: hole2_ = 3
integer, parameter :: particle2_= 4
BEGIN_PROVIDER [ integer, N_single_exc_bitmasks ]
implicit none
BEGIN_DOC
! Number of single excitation bitmasks
END_DOC
N_single_exc_bitmasks = 1
!TODO : Read from input!
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), single_exc_bitmask, (N_int, 2, N_single_exc_bitmasks) ]
implicit none
BEGIN_DOC
! single_exc_bitmask(:,1,i) is the bitmask for holes
! single_exc_bitmask(:,2,i) is the bitmask for particles
! for a given couple of hole/particle excitations i.
END_DOC
single_exc_bitmask(:,hole_,1) = HF_bitmask(:,1)
single_exc_bitmask(:,particle_,1) = not(HF_bitmask(:,2))
!TODO : Read from input!
END_PROVIDER
BEGIN_PROVIDER [ integer, N_double_exc_bitmasks ]
implicit none
BEGIN_DOC
! Number of double excitation bitmasks
END_DOC
N_double_exc_bitmasks = 1
!TODO : Read from input!
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), double_exc_bitmask, (N_int, 4, N_double_exc_bitmasks) ]
implicit none
BEGIN_DOC
! double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1
! double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1
! double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2
! double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2
! for a given couple of hole/particle excitations i.
END_DOC
double_exc_bitmask(:,hole_,1) = HF_bitmask(:,1)
double_exc_bitmask(:,particle_,1) = not(HF_bitmask(:,2))
double_exc_bitmask(:,hole2_,1) = HF_bitmask(:,1)
double_exc_bitmask(:,particle2_,1) = not(HF_bitmask(:,2))
!TODO : Read from input!
END_PROVIDER

34
src/Dets/s2.irp.f Normal file
View File

@ -0,0 +1,34 @@
subroutine get_s2(key_i,key_j,phase,Nint)
implicit none
BEGIN_DOC
! Returns <S^2>
END_DOC
integer, intent(in) :: Nint
integer, intent(in) :: key_i(Nint,2)
integer, intent(in) :: key_j(Nint,2)
double precision, intent(out) :: phase
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase_spsm
integer :: nup, i
phase = 0.d0
!$FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint)
select case (degree)
case(2)
call get_double_excitation(key_i,key_j,exc,phase_spsm,Nint)
if (exc(0,1,1) == 1) then ! Mono alpha + mono-beta
if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then
phase = -phase_spsm
endif
endif
case(0)
nup = 0
do i=1,Nint
nup += popcnt(iand(xor(key_i(i,1),key_i(i,2)),key_i(i,1)))
enddo
phase = dble(nup)
end select
end

949
src/Dets/slater_rules.irp.f Normal file
View File

@ -0,0 +1,949 @@
subroutine get_excitation_degree(key1,key2,degree,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Returns the excitation degree between two determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key1(Nint,2)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: degree
integer :: l
ASSERT (Nint > 0)
degree = popcnt(xor( key1(1,1), key2(1,1))) + &
popcnt(xor( key1(1,2), key2(1,2)))
!DEC$ NOUNROLL
do l=2,Nint
degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) + &
popcnt(xor( key1(l,2), key2(l,2)))
enddo
ASSERT (degree >= 0)
degree = ishft(degree,-1)
end
subroutine get_excitation(det1,det2,exc,degree,phase,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Returns the excitation operators between two determinants and the phase
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2)
integer(bit_kind), intent(in) :: det2(Nint,2)
integer, intent(out) :: exc(0:2,2,2)
integer, intent(out) :: degree
double precision, intent(out) :: phase
! exc(number,hole/particle,spin)
! ex :
! exc(0,1,1) = number of holes alpha
! exc(0,2,1) = number of particle alpha
! exc(0,2,2) = number of particle beta
! exc(1,2,1) = first particle alpha
! exc(1,1,1) = first hole alpha
! exc(1,2,2) = first particle beta
! exc(1,1,2) = first hole beta
ASSERT (Nint > 0)
!DIR$ FORCEINLINE
call get_excitation_degree(det1,det2,degree,Nint)
select case (degree)
case (3:)
degree = -1
return
case (2)
call get_double_excitation(det1,det2,exc,phase,Nint)
return
case (1)
call get_mono_excitation(det1,det2,exc,phase,Nint)
return
case(0)
return
end select
end
subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
implicit none
BEGIN_DOC
! Decodes the exc arrays returned by get_excitation.
! h1,h2 : Holes
! p1,p2 : Particles
! s1,s2 : Spins (1:alpha, 2:beta)
! degree : Degree of excitation
END_DOC
integer, intent(in) :: exc(0:2,2,2),degree
integer, intent(out) :: h1,h2,p1,p2,s1,s2
ASSERT (degree > 0)
ASSERT (degree < 3)
select case(degree)
case(2)
if (exc(0,1,1) == 2) then
h1 = exc(1,1,1)
h2 = exc(2,1,1)
p1 = exc(1,2,1)
p2 = exc(2,2,1)
s1 = 1
s2 = 1
else if (exc(0,1,2) == 2) then
h1 = exc(1,1,2)
h2 = exc(2,1,2)
p1 = exc(1,2,2)
p2 = exc(2,2,2)
s1 = 2
s2 = 2
else
h1 = exc(1,1,1)
h2 = exc(1,1,2)
p1 = exc(1,2,1)
p2 = exc(1,2,2)
s1 = 1
s2 = 2
endif
case(1)
if (exc(0,1,1) == 1) then
h1 = exc(1,1,1)
h2 = 0
p1 = exc(1,2,1)
p2 = 0
s1 = 1
s2 = 0
else
h1 = exc(1,1,2)
h2 = 0
p1 = exc(1,2,2)
p2 = 0
s1 = 2
s2 = 0
endif
case(0)
h1 = 0
p1 = 0
h2 = 0
p2 = 0
s1 = 0
s2 = 0
end select
end
subroutine get_double_excitation(det1,det2,exc,phase,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Returns the two excitation operators between two doubly excited determinants and the phase
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2)
integer(bit_kind), intent(in) :: det2(Nint,2)
integer, intent(out) :: exc(0:2,2,2)
double precision, intent(out) :: phase
integer :: tz
integer :: l, ispin, idx_hole, idx_particle, ishift
integer :: nperm
integer :: i,j,k,m,n
integer :: high, low
integer :: a,b,c,d
integer(bit_kind) :: hole, particle, tmp
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
ASSERT (Nint > 0)
nperm = 0
exc(0,1,1) = 0
exc(0,2,1) = 0
exc(0,1,2) = 0
exc(0,2,2) = 0
do ispin = 1,2
idx_particle = 0
idx_hole = 0
ishift = 1-bit_kind_size
do l=1,Nint
ishift = ishift + bit_kind_size
if (det1(l,ispin) == det2(l,ispin)) then
cycle
endif
tmp = xor( det1(l,ispin), det2(l,ispin) )
particle = iand(tmp, det2(l,ispin))
hole = iand(tmp, det1(l,ispin))
do while (particle /= 0_bit_kind)
tz = trailz(particle)
idx_particle = idx_particle + 1
exc(0,2,ispin) = exc(0,2,ispin) + 1
exc(idx_particle,2,ispin) = tz+ishift
particle = iand(particle,particle-1_bit_kind)
enddo
if (iand(exc(0,1,ispin),exc(0,2,ispin))==2) then ! exc(0,1,ispin)==2 or exc(0,2,ispin)==2
exit
endif
do while (hole /= 0_bit_kind)
tz = trailz(hole)
idx_hole = idx_hole + 1
exc(0,1,ispin) = exc(0,1,ispin) + 1
exc(idx_hole,1,ispin) = tz+ishift
hole = iand(hole,hole-1_bit_kind)
enddo
if (iand(exc(0,1,ispin),exc(0,2,ispin))==2) then ! exc(0,1,ispin)==2 or exc(0,2,ispin)
exit
endif
enddo
! TODO : Voir si il faut sortir i,n,k,m du case.
select case (exc(0,1,ispin))
case(0)
cycle
case(1)
low = min(exc(1,1,ispin), exc(1,2,ispin))
high = max(exc(1,1,ispin), exc(1,2,ispin))
ASSERT (low > 0)
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
n = iand(low,bit_kind_size-1) ! mod(low,bit_kind_size)
ASSERT (high > 0)
k = ishft(high-1,-bit_kind_shift)+1
m = iand(high,bit_kind_size-1)
if (j==k) then
nperm = nperm + popcnt(iand(det1(j,ispin), &
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
else
nperm = nperm + popcnt(iand(det1(k,ispin), &
ibset(0_bit_kind,m-1)-1_bit_kind)) + &
popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind))
do i=j+1,k-1
nperm = nperm + popcnt(det1(i,ispin))
end do
endif
case (2)
do i=1,2
low = min(exc(i,1,ispin), exc(i,2,ispin))
high = max(exc(i,1,ispin), exc(i,2,ispin))
ASSERT (low > 0)
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
n = iand(low,bit_kind_size-1) ! mod(low,bit_kind_size)
ASSERT (high > 0)
k = ishft(high-1,-bit_kind_shift)+1
m = iand(high,bit_kind_size-1)
if (j==k) then
nperm = nperm + popcnt(iand(det1(j,ispin), &
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
else
nperm = nperm + popcnt(iand(det1(k,ispin), &
ibset(0_bit_kind,m-1)-1_bit_kind)) + &
popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind))
do l=j+1,k-1
nperm = nperm + popcnt(det1(l,ispin))
end do
endif
enddo
a = min(exc(1,1,ispin), exc(1,2,ispin))
b = max(exc(1,1,ispin), exc(1,2,ispin))
c = min(exc(2,1,ispin), exc(2,2,ispin))
d = max(exc(2,1,ispin), exc(2,2,ispin))
if (c>a .and. c<b .and. d>b) then
nperm = nperm + 1
endif
exit
end select
enddo
phase = phase_dble(iand(nperm,1))
end
subroutine get_mono_excitation(det1,det2,exc,phase,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Returns the excitation operator between two singly excited determinants and the phase
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2)
integer(bit_kind), intent(in) :: det2(Nint,2)
integer, intent(out) :: exc(0:2,2,2)
double precision, intent(out) :: phase
integer :: tz
integer :: l, ispin, idx_hole, idx_particle, ishift
integer :: nperm
integer :: i,j,k,m,n
integer :: high, low
integer :: a,b,c,d
integer(bit_kind) :: hole, particle, tmp
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
ASSERT (Nint > 0)
nperm = 0
exc(0,1,1) = 0
exc(0,2,1) = 0
exc(0,1,2) = 0
exc(0,2,2) = 0
do ispin = 1,2
ishift = 1-bit_kind_size
do l=1,Nint
ishift = ishift + bit_kind_size
if (det1(l,ispin) == det2(l,ispin)) then
cycle
endif
tmp = xor( det1(l,ispin), det2(l,ispin) )
particle = iand(tmp, det2(l,ispin))
hole = iand(tmp, det1(l,ispin))
if (particle /= 0_bit_kind) then
tz = trailz(particle)
exc(0,2,ispin) = 1
exc(1,2,ispin) = tz+ishift
endif
if (hole /= 0_bit_kind) then
tz = trailz(hole)
exc(0,1,ispin) = 1
exc(1,1,ispin) = tz+ishift
endif
if ( iand(exc(0,1,ispin),exc(0,2,ispin)) /= 1) then ! exc(0,1,ispin)/=1 and exc(0,2,ispin) /= 1
cycle
endif
low = min(exc(1,1,ispin),exc(1,2,ispin))
high = max(exc(1,1,ispin),exc(1,2,ispin))
ASSERT (low > 0)
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
n = iand(low,bit_kind_size-1) ! mod(low,bit_kind_size)
ASSERT (high > 0)
k = ishft(high-1,-bit_kind_shift)+1
m = iand(high,bit_kind_size-1)
if (j==k) then
nperm = popcnt(iand(det1(j,ispin), &
iand(ibset(0_bit_kind,m-1)-1_bit_kind,ibclr(-1_bit_kind,n)+1_bit_kind)))
else
nperm = nperm + popcnt(iand(det1(k,ispin),ibset(0_bit_kind,m-1)-1_bit_kind)) +&
popcnt(iand(det1(j,ispin),ibclr(-1_bit_kind,n)+1_bit_kind))
do i=j+1,k-1
nperm = nperm + popcnt(det1(i,ispin))
end do
endif
phase = phase_dble(iand(nperm,1))
return
enddo
enddo
end
subroutine i_H_j(key_i,key_j,Nint,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|H|j> where i and j are determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij
integer :: exc(0:2,2,2)
integer :: degree
double precision :: get_mo_bielec_integral
integer :: m,n,p,q
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem, phase,phase_2
integer :: n_occ_alpha, n_occ_beta
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
hij = 0.d0
!DEC$ FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint)
select case (degree)
case (2)
call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral( &
exc(1,1,1), &
exc(1,1,2), &
exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map)
else if (exc(0,1,1) == 2) then
! Double alpha
hij = phase*(get_mo_bielec_integral( &
exc(1,1,1), &
exc(2,1,1), &
exc(1,2,1), &
exc(2,2,1) ,mo_integrals_map) - &
get_mo_bielec_integral( &
exc(1,1,1), &
exc(2,1,1), &
exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) )
else if (exc(0,1,2) == 2) then
! Double beta
hij = phase*(get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,2), &
exc(1,2,2), &
exc(2,2,2) ,mo_integrals_map) - &
get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,2), &
exc(2,2,2), &
exc(1,2,2) ,mo_integrals_map) )
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
m = exc(1,1,1)
p = exc(1,2,1)
do k = 1, elec_alpha_num
i = occ(k,1)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_beta_num
i = occ(k,2)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_alpha_num
hij = hij + mipi(occ(k,1)) - miip(occ(k,1))
enddo
do k = 1, elec_beta_num
hij = hij + mipi(occ(k,2))
enddo
else
! Mono beta
m = exc(1,1,2)
p = exc(1,2,2)
do k = 1, elec_beta_num
i = occ(k,2)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_alpha_num
i = occ(k,1)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True.
endif
enddo
do k = 1, elec_alpha_num
hij = hij + mipi(occ(k,1))
enddo
do k = 1, elec_beta_num
hij = hij + mipi(occ(k,2)) - miip(occ(k,2))
enddo
endif
hij = phase*(hij + mo_mono_elec_integral(m,p))
case (0)
hij = diag_H_mat_elem(key_i,Nint)
end select
end
subroutine i_H_psim(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
implicit none
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate
integer, intent(in) :: keys(Nint,2,Ndet_max)
double precision, intent(in) :: coef(Ndet_max,Nstate)
integer, intent(in) :: key(Nint,2)
double precision, intent(out) :: i_H_psi_array(Nstate)
integer :: i, ii,j
double precision :: phase
integer :: exc(0:2,2,2)
double precision :: hij
integer :: idx(0:Ndet)
i_H_psi_array = 0.d0
call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
do ii=1,idx(0)
i = idx(ii)
!DEC$ FORCEINLINE
call i_H_j(keys(1,1,i),key,Nint,hij)
do j = 1, Nstate
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
enddo
enddo
end
subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,sze_max,idx)
use bitmasks
implicit none
BEGIN_DOC
! Applies get_excitation_degree to an array of determinants
END_DOC
integer, intent(in) :: Nint, sze,sze_max
integer(bit_kind), intent(in) :: key1(Nint,2,sze_max)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: degree(sze)
integer, intent(out) :: idx(0:sze)
integer :: i,l
ASSERT (Nint > 0)
ASSERT (sze > 0)
ASSERT (sze_max >= sze)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree(l) = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))),-1)
if (degree(l) < 3) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree(l) = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))),-1)
if (degree(l) < 3) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree(l) = ishft( popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2))),-1)
if (degree(l) < 3) then
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree(l) = 0
!DEC$ LOOP COUNT MIN(4)
do l=1,Nint
degree(l) = degree(l)+ popcnt(xor( key1(l,1,i), key2(l,1))) +&
popcnt(xor( key1(l,2,i), key2(l,2)))
enddo
degree(l) = ishft(degree(l),-1)
if (degree(l) < 3) then
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
end
subroutine filter_connected(key1,key2,Nint,sze,sze_max,idx)
use bitmasks
implicit none
BEGIN_DOC
! Filters out the determinants that are not connected by H
END_DOC
integer, intent(in) :: Nint, sze,sze_max
integer(bit_kind), intent(in) :: key1(Nint,2,sze_max)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,j,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze > 0)
ASSERT (sze_max >= sze)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do j=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +&
popcnt(xor( key1(j,2,i), key2(j,2)))
if (degree_x2 > 4) then
exit
endif
enddo
if (degree_x2 <= 5) then
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
end
subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,sze_max,idx)
use bitmasks
implicit none
integer, intent(in) :: Nint, sze,sze_max
integer(bit_kind), intent(in) :: key1(Nint,2,sze_max)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze > 0)
ASSERT (sze_max >= sze)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do l=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +&
popcnt(xor( key1(l,2,i), key2(l,2)))
if (degree_x2 > 4) then
exit
endif
enddo
if (degree_x2 <= 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
endif
idx(0) = l-1
end
double precision function diag_H_mat_elem(det_in,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes <i|H|i>
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: det_in(Nint,2)
integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: particle(Nint,2)
integer :: i, nexc(2), ispin
integer :: occ_particle(Nint*bit_kind_size,2)
integer :: occ_hole(Nint*bit_kind_size,2)
integer(bit_kind) :: det_tmp(Nint,2)
integer :: na, nb
ASSERT (Nint > 0)
ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num)
nexc(1) = 0
nexc(2) = 0
do i=1,Nint
hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1))
hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2))
particle(i,1) = iand(hole(i,1),det_in(i,1))
particle(i,2) = iand(hole(i,2),det_in(i,2))
hole(i,1) = iand(hole(i,1),ref_bitmask(i,1))
hole(i,2) = iand(hole(i,2),ref_bitmask(i,2))
nexc(1) += popcnt(hole(i,1))
nexc(2) += popcnt(hole(i,2))
enddo
diag_H_mat_elem = ref_bitmask_energy
if (nexc(1)+nexc(2) == 0) then
return
endif
!call debug_det(det_in,Nint)
integer :: tmp
call bitstring_to_list(particle(1,1), occ_particle(1,1), tmp, Nint)
ASSERT (tmp == nexc(1))
call bitstring_to_list(particle(1,2), occ_particle(1,2), tmp, Nint)
ASSERT (tmp == nexc(2))
call bitstring_to_list(hole(1,1), occ_hole(1,1), tmp, Nint)
ASSERT (tmp == nexc(1))
call bitstring_to_list(hole(1,2), occ_hole(1,2), tmp, Nint)
ASSERT (tmp == nexc(2))
det_tmp = ref_bitmask
do ispin=1,2
na = elec_num_tab(ispin)
nb = elec_num_tab(iand(ispin,1)+1)
do i=1,nexc(ispin)
!DIR$ FORCEINLINE
call ac_operator( occ_particle(i,ispin), ispin, det_tmp, diag_H_mat_elem, Nint,na,nb)
!DIR$ FORCEINLINE
call a_operator ( occ_hole (i,ispin), ispin, det_tmp, diag_H_mat_elem, Nint,na,nb)
enddo
enddo
end
subroutine a_operator(iorb,ispin,key,hjj,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Needed for diag_H_mat_elem
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer, intent(inout) :: na, nb
integer(bit_kind), intent(inout) :: key(Nint,2)
double precision, intent(inout) :: hjj
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k > 0)
l = iorb - ishft(k-1,bit_kind_shift)-1
key(k,ispin) = ibclr(key(k,ispin),l)
other_spin = iand(ispin,1)+1
!DIR$ FORCEINLINE
call get_occ_from_key(key,occ,Nint)
na -= 1
hjj -= mo_mono_elec_integral(iorb,iorb)
! Same spin
do i=1,na
hjj -= mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
hjj -= mo_bielec_integral_jj(occ(i,other_spin),iorb)
enddo
end
subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Needed for diag_H_mat_elem
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer, intent(inout) :: na, nb
integer(bit_kind), intent(inout) :: key(Nint,2)
double precision, intent(inout) :: hjj
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
integer :: tmp
!DIR$ FORCEINLINE
call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint)
ASSERT (tmp == elec_alpha_num)
!DIR$ FORCEINLINE
call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint)
ASSERT (tmp == elec_beta_num)
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k > 0)
l = iorb - ishft(k-1,bit_kind_shift)-1
key(k,ispin) = ibset(key(k,ispin),l)
other_spin = iand(ispin,1)+1
hjj += mo_mono_elec_integral(iorb,iorb)
! Same spin
do i=1,na
hjj += mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
hjj += mo_bielec_integral_jj(occ(i,other_spin),iorb)
enddo
na += 1
end
subroutine get_occ_from_key(key,occ,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Returns a list of occupation numbers from a bitstring
END_DOC
integer(bit_kind), intent(in) :: key(Nint,2)
integer , intent(in) :: Nint
integer , intent(out) :: occ(Nint*bit_kind_size,2)
integer :: tmp
call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint)
call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint)
end

15
src/Dets/utils.irp.f Normal file
View File

@ -0,0 +1,15 @@
BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(n_det,n_det) ]
implicit none
BEGIN_DOC
! H matrix on the basis of the slater deter;inants defined by psi_det
END_DOC
integer :: i,j
double precision :: hij
do i =1,N_det
do j =i,N_det
call i_H_j(psi_det(1,1,i),psi_det(1,1,j),N_int,hij)
H_matrix_all_dets(i,j) = hij
H_matrix_all_dets(j,i) = hij
enddo
enddo
END_PROVIDER

22
src/Electrons/README.rst
View File

@ -10,7 +10,7 @@ Assumptions
===========
.. Do not edit this section. It was auto-generated from the
.. ASSUMPTIONS.rst file.
.. NEEDED_MODULES file.
* ``elec_num`` >= 0
* ``elec_alpha_num`` >= 0
@ -28,3 +28,23 @@ Needed Modules
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`elec_alpha_num <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons/electrons.irp.f#L1>`_
Numbers of alpha ("up") , beta ("down") and total electrons
`elec_beta_num <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons/electrons.irp.f#L2>`_
Numbers of alpha ("up") , beta ("down") and total electrons
`elec_num <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons/electrons.irp.f#L3>`_
Numbers of alpha ("up") , beta ("down") and total electrons
`elec_num_tab <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons/electrons.irp.f#L4>`_
Numbers of alpha ("up") , beta ("down") and total electrons

25
src/Ezfio_files/README.rst
View File

@ -5,3 +5,28 @@ Ezfio_files Module
This modules essentially contains the name of the EZFIO directory in the
``ezfio_filename`` variable. This is read as the first argument of the
command-line, or as the ``QPACKAGE_INPUT`` environment variable.
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`ezfio_filename <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/ezfio.irp.f#L1>`_
Name of EZFIO file. It is obtained from the QPACKAGE_INPUT environment
variable if it is set, or as the 1st argument of the command line.
`getunitandopen <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/get_unit_and_open.irp.f#L1>`_
:f:
file name
.br
:mode:
'R' : READ, UNFORMATTED
'W' : WRITE, UNFORMATTED
'r' : READ, FORMATTED
'w' : WRITE, FORMATTED
'a' : APPEND, FORMATTED
'x' : READ/WRITE, FORMATTED
.br

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

219
src/Hartree_Fock/Fock_matrix.irp.f Normal file
View File

@ -0,0 +1,219 @@
BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo, (mo_tot_num)]
implicit none
BEGIN_DOC
! Fock matrix on the MO basis.
! For open shells, the ROHF Fock Matrix is
!
! | F-K | F + K/2 | F |
! |---------------------------------|
! | F + K/2 | F | F - K/2 |
! |---------------------------------|
! | F | F - K/2 | F + K |
!
! F = 1/2 (Fa + Fb)
!
! K = Fb - Fa
!
END_DOC
integer :: i,j,n
if (elec_alpha_num == elec_beta_num) then
Fock_matrix_mo = Fock_matrix_alpha_mo
else
do j=1,elec_beta_num
! F-K
do i=1,elec_beta_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
- (Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
enddo
! F+K/2
do i=elec_beta_num+1,elec_alpha_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
+ 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
enddo
! F
do i=elec_alpha_num+1, mo_tot_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))
enddo
enddo
do j=elec_beta_num+1,elec_alpha_num
! F+K/2
do i=1,elec_beta_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
+ 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
enddo
! F
do i=elec_beta_num+1,elec_alpha_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))
enddo
! F-K/2
do i=elec_alpha_num+1, mo_tot_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
- 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
enddo
enddo
do j=elec_alpha_num+1, mo_tot_num
! F
do i=1,elec_beta_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))
enddo
! F-K/2
do i=elec_beta_num+1,elec_alpha_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
- 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
enddo
! F+K
do i=elec_alpha_num+1,mo_tot_num
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j)) &
+ (Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
enddo
enddo
endif
do i = 1, mo_tot_num
Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_ao, (ao_num_align, ao_num) ]
&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_ao, (ao_num_align, ao_num) ]
implicit none
BEGIN_DOC
! Alpha Fock matrix in AO basis set
END_DOC
integer :: i,j,k,l,k1,kmax
double precision, allocatable :: ao_ints_val(:)
integer, allocatable :: ao_ints_idx(:)
double precision :: integral
double precision :: ao_bielec_integral
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_ints_idx, ao_ints_val
if (do_direct_integrals) then
PROVIDE all_utils ao_overlap_abs ao_integrals_threshold
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral) &
!$OMP SHARED(ao_num,Fock_matrix_alpha_ao,ao_mono_elec_integral,&
!$OMP ao_num_align,Fock_matrix_beta_ao,HF_density_matrix_ao_alpha, &
!$OMP HF_density_matrix_ao_beta)
!$OMP DO SCHEDULE(GUIDED)
do j=1,ao_num
do i=1,j
Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j)
Fock_matrix_beta_ao (i,j) = ao_mono_elec_integral(i,j)
do l=1,ao_num
do k=1,ao_num
PROVIDE HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
if ((abs(HF_density_matrix_ao_alpha(k,l)) > 1.d-9).or. &
(abs(HF_density_matrix_ao_beta (k,l)) > 1.d-9)) then
integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta (k,l)) * ao_bielec_integral(k,l,i,j)
Fock_matrix_alpha_ao(i,j) += integral
Fock_matrix_beta_ao (i,j) += integral
integral = ao_bielec_integral(k,j,i,l)
Fock_matrix_alpha_ao(i,j) -= HF_density_matrix_ao_alpha(k,l)*integral
Fock_matrix_beta_ao (i,j) -= HF_density_matrix_ao_beta (k,l)*integral
endif
enddo
enddo
Fock_matrix_alpha_ao(j,i) = Fock_matrix_alpha_ao(i,j)
Fock_matrix_beta_ao (j,i) = Fock_matrix_beta_ao (i,j)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP END PARALLEL
else
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral,ao_ints_val,ao_ints_idx,kmax) &
!$OMP SHARED(ao_num,Fock_matrix_alpha_ao,ao_mono_elec_integral,&
!$OMP ao_num_align,Fock_matrix_beta_ao,HF_density_matrix_ao_alpha, &
!$OMP HF_density_matrix_ao_beta)
allocate(ao_ints_idx(ao_num_align),ao_ints_val(ao_num_align))
!$OMP DO SCHEDULE(GUIDED)
do j=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,ao_num
Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j)
Fock_matrix_beta_ao (i,j) = ao_mono_elec_integral(i,j)
enddo
do l=1,ao_num
do i=1,ao_num
call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax)
!DIR$ VECTOR ALIGNED
do k1=1,kmax
k = ao_ints_idx(k1)
integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * ao_ints_val(k1)
Fock_matrix_alpha_ao(i,j) += integral
Fock_matrix_beta_ao (i,j) += integral
integral = ao_ints_val(k1)
Fock_matrix_alpha_ao(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral
Fock_matrix_beta_ao (l,j) -= HF_density_matrix_ao_beta (k,i) * integral
enddo
enddo
enddo
enddo
!$OMP END DO
deallocate(ao_ints_val,ao_ints_idx)
!$OMP END PARALLEL
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Fock matrix on the MO basis
END_DOC
double precision, allocatable :: T(:,:)
allocate ( T(ao_num_align,mo_tot_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, Fock_matrix_alpha_ao,size(Fock_matrix_alpha_ao,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, ao_num_align)
call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, size(T,1), &
0.d0, Fock_matrix_alpha_mo, mo_tot_num_align)
deallocate(T)
END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Fock matrix on the MO basis
END_DOC
double precision, allocatable :: T(:,:)
allocate ( T(ao_num_align,mo_tot_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, Fock_matrix_beta_ao,size(Fock_matrix_beta_ao,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, ao_num_align)
call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, size(T,1), &
0.d0, Fock_matrix_beta_mo, mo_tot_num_align)
deallocate(T)
END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_energy ]
implicit none
BEGIN_DOC
! Hartree-Fock energy
END_DOC
HF_energy = nuclear_repulsion + ref_bitmask_energy
END_PROVIDER

87
src/Hartree_Fock/HF_density_matrix_ao.irp.f Normal file
View File

@ -0,0 +1,87 @@
BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Alpha and Beta density matrix in the AO basis
END_DOC
integer :: i,j,k,l1,l2
integer, allocatable :: mo_occ(:,:)
allocate ( mo_occ(elec_alpha_num,2) )
call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
ASSERT ( j==elec_alpha_num )
call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
ASSERT ( j==elec_beta_num )
do j=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,ao_num_align
HF_density_matrix_ao_alpha(i,j) = 0.d0
HF_density_matrix_ao_beta (i,j) = 0.d0
enddo
do k=1,elec_beta_num
l1 = mo_occ(k,1)
l2 = mo_occ(k,2)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
mo_coef(i,l1) * mo_coef(j,l1)
HF_density_matrix_ao_beta (i,j) = HF_density_matrix_ao_beta (i,j) +&
mo_coef(i,l2) * mo_coef(j,l2)
enddo
enddo
do k=elec_beta_num+1,elec_alpha_num
l1 = mo_occ(k,1)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
mo_coef(i,l1) * mo_coef(j,l1)
enddo
enddo
enddo
deallocate(mo_occ)
END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Density matrix in the AO basis
END_DOC
integer :: i,j,k,l1,l2
integer, allocatable :: mo_occ(:,:)
allocate ( mo_occ(elec_alpha_num,2) )
call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
ASSERT ( j==elec_alpha_num )
call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
ASSERT ( j==elec_beta_num )
do j=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,ao_num_align
HF_density_matrix_ao(i,j) = 0.d0
enddo
do k=1,elec_beta_num
l1 = mo_occ(k,1)
l2 = mo_occ(k,2)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + &
mo_coef(i,l1) * mo_coef(j,l1) + &
mo_coef(i,l2) * mo_coef(j,l2)
enddo
enddo
do k=elec_beta_num+1,elec_alpha_num
l1 = mo_occ(k,1)
!DIR$ VECTOR ALIGNED
do i=1,ao_num
HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + &
mo_coef(i,l1) * mo_coef(j,l1)
enddo
enddo
enddo
deallocate(mo_occ)
END_PROVIDER

8
src/Hartree_Fock/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/Hartree_Fock/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Electrons Ezfio_files MonoInts MOs Nuclei Output Utils

94
src/Hartree_Fock/README.rst
View File

@ -20,3 +20,97 @@ Needed Modules
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`fock_matrix_alpha_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L83>`_
Alpha Fock matrix in AO basis set
`fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L172>`_
Fock matrix on the MO basis
`fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L84>`_
Alpha Fock matrix in AO basis set
`fock_matrix_beta_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L192>`_
Fock matrix on the MO basis
`fock_matrix_diag_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L2>`_
Fock matrix on the MO basis.
For open shells, the ROHF Fock Matrix is
.br
| F-K | F + K/2 | F |
|---------------------------------|
| F + K/2 | F | F - K/2 |
|---------------------------------|
| F | F - K/2 | F + K |
.br
F = 1/2 (Fa + Fb)
.br
K = Fb - Fa
.br
`fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L1>`_
Fock matrix on the MO basis.
For open shells, the ROHF Fock Matrix is
.br
| F-K | F + K/2 | F |
|---------------------------------|
| F + K/2 | F | F - K/2 |
|---------------------------------|
| F | F - K/2 | F + K |
.br
F = 1/2 (Fa + Fb)
.br
K = Fb - Fa
.br
`hf_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L211>`_
Hartree-Fock energy
`hf_density_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L46>`_
Density matrix in the AO basis
`hf_density_matrix_ao_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L1>`_
Alpha and Beta density matrix in the AO basis
`hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L2>`_
Alpha and Beta density matrix in the AO basis
`diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L1>`_
Diagonal Fock matrix in the MO basis
`eigenvectors_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L2>`_
Diagonal Fock matrix in the MO basis
`scf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L1>`_
None
`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L41>`_
If True, compute integrals on the fly
`n_it_scf_max <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L22>`_
Maximum number of SCF iterations
`thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L1>`_
Threshold on the convergence of the Hartree Fock energy
`bi_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L5>`_
Energy of the reference bitmask used in Slater rules
`kinetic_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L3>`_
Energy of the reference bitmask used in Slater rules
`mono_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L2>`_
Energy of the reference bitmask used in Slater rules
`nucl_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L4>`_
Energy of the reference bitmask used in Slater rules
`ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L1>`_
Energy of the reference bitmask used in Slater rules

20
src/Hartree_Fock/diagonalize_fock.irp.f Normal file
View File

@ -0,0 +1,20 @@
BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Diagonal Fock matrix in the MO basis
END_DOC
double precision, allocatable :: R(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: R
allocate(R(mo_tot_num_align,mo_tot_num))
call lapack_diag(diagonal_Fock_matrix_mo,R,Fock_matrix_mo,size(Fock_matrix_mo,1),&
mo_tot_num)
call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num,1.d0,mo_coef,size(mo_coef,1),&
R,size(R,1),0.d0,eigenvectors_Fock_matrix_mo,size(eigenvectors_Fock_matrix_mo,1))
deallocate(R)
END_PROVIDER

5
src/Hartree_Fock/hartree_fock.ezfio_config Normal file
View File

@ -0,0 +1,5 @@
hartree_fock
thresh_scf double precision
n_it_scf_max integer
diis logical

37
src/Hartree_Fock/mo_SCF_iterations.irp.f Normal file
View File

@ -0,0 +1,37 @@
program scf_iteration
use bitmasks
implicit none
double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral
double precision :: E0
integer :: i_it
E0 = HF_energy
i_it = 0
n_it_scf_max = 100
SCF_energy_before = huge(1.d0)
SCF_energy_after = E0
print *, E0
mo_label = "Canonical"
thresh_SCF = 1.d-10
do while (i_it < 10 .or. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF)
SCF_energy_before = SCF_energy_after
mo_coef = eigenvectors_Fock_matrix_mo
TOUCH mo_coef mo_label
call clear_mo_map
SCF_energy_after = HF_energy
print*,SCF_energy_after
i_it +=1
if(i_it > n_it_scf_max)exit
enddo
if (i_it == n_it_scf_max) then
stop 'Failed'
endif
if (SCF_energy_after - E0 > thresh_SCF) then
stop 'Failed'
endif
mo_label = "Canonical"
TOUCH mo_label mo_coef
call save_mos
end

58
src/Hartree_Fock/options.irp.f Normal file
View File

@ -0,0 +1,58 @@
BEGIN_PROVIDER [ double precision,thresh_SCF ]
implicit none
BEGIN_DOC
! Threshold on the convergence of the Hartree Fock energy
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_Hartree_Fock_thresh_SCF(has)
if (has) then
call ezfio_get_Hartree_Fock_thresh_SCF(thresh_SCF)
else
thresh_SCF = 1.d-10
call ezfio_set_Hartree_Fock_thresh_SCF(thresh_SCF)
endif
call write_time(output_Hartree_Fock)
call write_double(output_Hartree_Fock, thresh_SCF, &
'thresh_SCF')
END_PROVIDER
BEGIN_PROVIDER [ integer, n_it_scf_max]
implicit none
BEGIN_DOC
! Maximum number of SCF iterations
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_Hartree_Fock_n_it_scf_max (has)
if (has) then
call ezfio_get_Hartree_Fock_n_it_scf_max(n_it_scf_max)
else
n_it_scf_max = 30
call ezfio_set_Hartree_Fock_n_it_scf_max(n_it_scf_max)
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, do_DIIS ]
implicit none
BEGIN_DOC
! If True, compute integrals on the fly
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_Hartree_Fock_DIIS(has)
if (has) then
call ezfio_get_Hartree_Fock_DIIS(do_DIIS)
else
do_DIIS = .False.
call ezfio_set_Hartree_Fock_DIIS(do_DIIS)
endif
END_PROVIDER

57
src/Hartree_Fock/ref_bitmask.irp.f Normal file
View File

@ -0,0 +1,57 @@
BEGIN_PROVIDER [ double precision, ref_bitmask_energy ]
&BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ]
&BEGIN_PROVIDER [ double precision, kinetic_ref_bitmask_energy ]
&BEGIN_PROVIDER [ double precision, nucl_elec_ref_bitmask_energy ]
&BEGIN_PROVIDER [ double precision, bi_elec_ref_bitmask_energy ]
use bitmasks
implicit none
BEGIN_DOC
! Energy of the reference bitmask used in Slater rules
END_DOC
integer :: occ(N_int*bit_kind_size,2)
integer :: i,j
call bitstring_to_list(ref_bitmask(1,1), occ(1,1), i, N_int)
call bitstring_to_list(ref_bitmask(1,2), occ(1,2), i, N_int)
ref_bitmask_energy = 0.d0
mono_elec_ref_bitmask_energy = 0.d0
kinetic_ref_bitmask_energy = 0.d0
nucl_elec_ref_bitmask_energy = 0.d0
bi_elec_ref_bitmask_energy = 0.d0
do i = 1, elec_beta_num
ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1)) + mo_mono_elec_integral(occ(i,2),occ(i,2))
kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1)) + mo_kinetic_integral(occ(i,2),occ(i,2))
nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1)) + mo_nucl_elec_integral(occ(i,2),occ(i,2))
enddo
do i = elec_beta_num+1,elec_alpha_num
ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1))
kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1))
nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1))
enddo
do j= 1, elec_alpha_num
do i = j+1, elec_alpha_num
bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1))
ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1))
enddo
enddo
do j= 1, elec_beta_num
do i = j+1, elec_beta_num
bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2))
ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2))
enddo
do i= 1, elec_alpha_num
bi_elec_ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2))
ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2))
enddo
enddo
mono_elec_ref_bitmask_energy = kinetic_ref_bitmask_energy + nucl_elec_ref_bitmask_energy
END_PROVIDER

33
src/Hartree_Fock/tests/Makefile Normal file
View File

@ -0,0 +1,33 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90+= -I tests
REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f))
.PHONY: clean executables serial_tests parallel_tests
all: clean executables serial_tests parallel_tests
parallel_tests: $(REF_FILES)
@echo ; echo " ---- Running parallel tests ----" ; echo
@OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py
serial_tests: $(REF_FILES)
@echo ; echo " ---- Running serial tests ----" ; echo
@OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py
executables: $(wildcard *.irp.f) veryclean
$(MAKE) -C ..
%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables
$(QPACKAGE_ROOT)/scripts/create_test_ref.sh $*
clean:
$(MAKE) -C .. clean
veryclean:
$(MAKE) -C .. veryclean

9
src/Hartree_Fock/tests/hf_energy.irp.f Normal file
View File

@ -0,0 +1,9 @@
program test
implicit none
print *, 'HF energy :', HF_energy
print *, 'ref_bitmask_energy :', ref_bitmask_energy
print *, 'mono_elec_ref_bitmask_energy :', mono_elec_ref_bitmask_energy
print *, 'kinetic_ref_bitmask_energy :', kinetic_ref_bitmask_energy
print *, 'nucl_elec_ref_bitmask_energy :', nucl_elec_ref_bitmask_energy
print *, 'bi_elec_ref_bitmask_energy :', bi_elec_ref_bitmask_energy
end

1690
src/Hartree_Fock/tests/hf_energy.ref Normal file
View File

File diff suppressed because it is too large Load Diff

30
src/Hartree_Fock_AOs/README.rst Normal file
View File

@ -0,0 +1,30 @@
===================
Hartree-Fock Module
===================
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.

35
src/MOGuess/README.rst
View File

@ -2,3 +2,38 @@
MOGuess Module
==============
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L1>`_
None
`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L2>`_
matrix of the coefficients of the mos generated by the
orthonormalization by the S^{-1/2} canonical transformation of the aos
ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital
`ao_ortho_lowdin_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L26>`_
overlap matrix of the ao_ortho_lowdin
supposed to be the Identity
`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L1>`_
None

114
src/MOGuess/localize.irp.f
View File

@ -1,114 +0,0 @@
!TODO Ecrire un cholesky avec bitmask
subroutine localize_mos(mask, nint)
implicit none
use bitmasks
integer, intent(in) :: nint
integer(bit_kind), intent(in) :: mask(nint)
integer :: i,j,k,l
double precision, allocatable :: DM(:,:)
double precision, allocatable :: mo_coef_new(:,:), R(:,:)
integer :: n
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: DM, mo_coef_new, R
integer :: rank
integer, parameter :: n_core = 2
allocate(R(mo_tot_num,mo_tot_num))
allocate(DM(ao_num_align,ao_num))
allocate(mo_coef_new(ao_num_align,mo_tot_num))
n = ao_num
mo_coef_new = mo_coef
BEGIN_TEMPLATE
DM = 0.d0
if ($START < $END) then
do k=$START, $END
do j=1,n
!DEC$ VECTOR ALIGNED
do i=1,n
DM(i,j) += mo_coef_new(i,k)*mo_coef_new(j,k)
enddo
enddo
enddo
call cholesky_mo(n,$END-$START+1,DM,mo_coef_new(1,$START),size(mo_coef_new,1),-1.d0,rank)
endif
SUBST [ START, END ]
1 ; n_core ;;
END_TEMPLATE
deallocate(DM)
call find_rotation(mo_coef,ao_num_align,mo_coef_new,ao_num,R,mo_tot_num)
mo_coef = mo_coef_new
deallocate(mo_coef_new)
double precision,allocatable :: mo_energy_new(:)
integer, allocatable :: iorder(:)
allocate(mo_energy_new(mo_tot_num),iorder(mo_tot_num))
do i=1,mo_tot_num
iorder(i) = i
mo_energy_new(i) = 0.d0
do k=1,mo_tot_num
mo_energy_new(i) += R(k,i)*R(k,i)*mo_energy(k)
enddo
enddo
mo_energy = mo_energy_new
call dsort(mo_energy(1),iorder(1),n_core)
allocate (mo_coef_new(ao_num_align,mo_tot_num))
mo_coef_new = mo_coef
do j=1,mo_tot_num
do i=1,ao_num
mo_coef(i,j) = mo_coef_new(i,iorder(j))
enddo
enddo
deallocate (mo_coef_new,R)
deallocate(mo_energy_new,iorder)
mo_label = 'localized'
SOFT_TOUCH mo_coef mo_energy mo_label
end
subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank)
implicit none
integer, intent(in) :: n,m, LDC
double precision, intent(in) :: P(LDC,n)
double precision, intent(out) :: C(LDC,m)
double precision, intent(in) :: tol_in
integer, intent(out) :: rank
integer :: info
integer :: i,k
integer :: ipiv(n)
double precision:: tol
double precision, allocatable :: W(:,:), work(:)
!DEC$ ATTRIBUTES ALIGN: 32 :: W
!DEC$ ATTRIBUTES ALIGN: 32 :: work
!DEC$ ATTRIBUTES ALIGN: 32 :: ipiv
allocate(W(LDC,n),work(2*n))
tol=tol_in
info = 0
do i=1,n
do k=1,i
W(i,k) = P(i,k)
enddo
do k=i+1,n
W(i,k) = 0.
enddo
enddo
call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info )
do i=1,n
do k=1,min(m,rank)
C(ipiv(i),k) = W(i,k)
enddo
enddo
deallocate(W,work)
end

2
src/MOs/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Ezfio_files Nuclei Utils
AOs Ezfio_files Nuclei Output Utils

32
src/MOs/README.rst
View File

@ -25,5 +25,37 @@ Needed Modules
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`mo_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L22>`_
Molecular orbital coefficients on AO basis set
mo_coef(i,j) = coefficient of the ith ao on the jth mo
mo_label : Label characterizing the MOS (local, canonical, natural, etc)
`mo_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L61>`_
Molecular orbital coefficients on AO basis set
`mo_label <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L23>`_
Molecular orbital coefficients on AO basis set
mo_coef(i,j) = coefficient of the ith ao on the jth mo
mo_label : Label characterizing the MOS (local, canonical, natural, etc)
`mo_tot_num <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L1>`_
Total number of molecular orbitals and the size of the keys corresponding
`mo_tot_num_align <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L12>`_
Aligned variable for dimensioning of arrays
`mo_as_eigvectors_of_mo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L21>`_
None
`save_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L1>`_
None

2
src/MOs/utils.irp.f
View File

@ -3,7 +3,7 @@ subroutine save_mos
double precision, allocatable :: buffer(:,:)
integer :: i,j
call system('save_current_mos.sh '//trim(ezfio_filename))
call system('$QPACKAGE_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename))
call ezfio_set_mo_basis_mo_label(mo_label)
allocate ( buffer(ao_num,mo_tot_num) )

4
src/Makefile.common
View File

@ -113,9 +113,9 @@ clean_links:
endif
LIB+=$(EZFIO) $(MKL)
IRPF90+=$(patsubst %, -I %, $(INCLUDE_DIRS))
IRPF90+=$(patsubst %, -I %, $(INCLUDE_DIRS)) $(IRPF90_FLAGS)
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO) NEEDED_MODULES
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO) NEEDED_MODULES $(wildcard *.py)
$(IRPF90)
Makefile.depend: Makefile

99
src/MonoInts/README.rst
View File

@ -11,3 +11,102 @@ Needed Modules
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`ao_mono_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L122>`_
array of the mono electronic hamiltonian on the AOs basis
: sum of the kinetic and nuclear electronic potential
`ao_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L1>`_
Overlap between atomic basis functions:
:math:`\int \chi_i(r) \chi_j(r) dr)`
`ao_overlap_abs <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L65>`_
Overlap between absolute value of atomic basis functions:
:math:`\int |\chi_i(r)| |\chi_j(r)| dr)`
`ao_overlap_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L2>`_
Overlap between atomic basis functions:
:math:`\int \chi_i(r) \chi_j(r) dr)`
`ao_overlap_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L3>`_
Overlap between atomic basis functions:
:math:`\int \chi_i(r) \chi_j(r) dr)`
`ao_overlap_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L4>`_
Overlap between atomic basis functions:
:math:`\int \chi_i(r) \chi_j(r) dr)`
`check_ortho <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/check_orthonormality.irp.f#L1>`_
None
`do_print <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/check_orthonormality.irp.f#L11>`_
None
`n_pt_max_i_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L2>`_
None
`n_pt_max_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L1>`_
None
`ao_deriv2_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L1>`_
second derivatives matrix elements in the ao basis
.. math::
.br
{\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
`ao_deriv2_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L2>`_
second derivatives matrix elements in the ao basis
.. math::
.br
{\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
`ao_deriv2_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L3>`_
second derivatives matrix elements in the ao basis
.. math::
.br
{\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
`ao_kinetic_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L125>`_
array of the priminitve basis kinetic integrals
\langle \chi_i |\hat{T}| \chi_j \rangle
`mo_kinetic_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_mo_ints.irp.f#L1>`_
None
`mo_mono_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/mo_mono_ints.irp.f#L35>`_
array of the mono electronic hamiltonian on the MOs basis
: sum of the kinetic and nuclear electronic potential
`mo_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/mo_mono_ints.irp.f#L1>`_
None
`orthonormalize_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/orthonormalize.irp.f#L1>`_
None
`ao_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L1>`_
interaction nuclear electron
`give_polynom_mult_center_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L157>`_
None
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L285>`_
None
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L357>`_
None
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L428>`_
None
`nai_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L82>`_
None
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L409>`_
None
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L473>`_
None
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L457>`_
None
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L486>`_
None
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L502>`_
None
`mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_mo_ints.irp.f#L1>`_
None
`save_ortho_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/save_ortho_mos.irp.f#L1>`_
None

14
src/MonoInts/ao_mono_ints.irp.f
View File

@ -119,3 +119,17 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
implicit none
integer :: i,j,n,l
BEGIN_DOC
! array of the mono electronic hamiltonian on the AOs basis
! : sum of the kinetic and nuclear electronic potential
END_DOC
do j = 1, ao_num
do i = 1, ao_num
ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j)
enddo
enddo
END_PROVIDER

4
src/MonoInts/mo_mono_ints.irp.f
View File

@ -39,8 +39,8 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to
! array of the mono electronic hamiltonian on the MOs basis
! : sum of the kinetic and nuclear electronic potential
END_DOC
do i = 1, mo_tot_num
do j = 1, mo_tot_num
do j = 1, mo_tot_num
do i = 1, mo_tot_num
mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j)
enddo
enddo

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD

54
src/Nuclei/README.rst
View File

@ -16,3 +16,57 @@ Needed Modules
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`nucl_charge <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L23>`_
Nuclear charges
`nucl_coord <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L55>`_
Nuclear coordinates in the format (:, {x,y,z})
`nucl_coord_transp <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L110>`_
Transposed array of nucl_coord
`nucl_dist <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L129>`_
nucl_dist : Nucleus-nucleus distances
nucl_dist_2 : Nucleus-nucleus distances squared
nucl_dist_vec : Nucleus-nucleus distances vectors
`nucl_dist_2 <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L125>`_
nucl_dist : Nucleus-nucleus distances
nucl_dist_2 : Nucleus-nucleus distances squared
nucl_dist_vec : Nucleus-nucleus distances vectors
`nucl_dist_vec_x <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L126>`_
nucl_dist : Nucleus-nucleus distances
nucl_dist_2 : Nucleus-nucleus distances squared
nucl_dist_vec : Nucleus-nucleus distances vectors
`nucl_dist_vec_y <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L127>`_
nucl_dist : Nucleus-nucleus distances
nucl_dist_2 : Nucleus-nucleus distances squared
nucl_dist_vec : Nucleus-nucleus distances vectors
`nucl_dist_vec_z <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L128>`_
nucl_dist : Nucleus-nucleus distances
nucl_dist_2 : Nucleus-nucleus distances squared
nucl_dist_vec : Nucleus-nucleus distances vectors
`nucl_label <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L41>`_
Nuclear labels
`nucl_num <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L1>`_
Number of nuclei
`nucl_num_aligned <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L2>`_
Number of nuclei
`nuclear_repulsion <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei/nuclei.irp.f#L171>`_
Nuclear repulsion energy

330
src/Nuclei/nuclei.irp.f
View File

@ -1,174 +1,198 @@
BEGIN_PROVIDER [ integer, nucl_num ]
&BEGIN_PROVIDER [ integer, nucl_num_aligned ]
implicit none
BEGIN_DOC
! Number of nuclei
END_DOC
PROVIDE ezfio_filename
nucl_num = 0
call ezfio_get_nuclei_nucl_num(nucl_num)
ASSERT (nucl_num > 0)
integer :: align_double
nucl_num_aligned = align_double(nucl_num)
implicit none
BEGIN_DOC
! Number of nuclei
END_DOC
PROVIDE ezfio_filename
nucl_num = 0
logical :: has
call ezfio_has_nuclei_nucl_num(has)
if (has) then
call ezfio_get_nuclei_nucl_num(nucl_num)
else
print *, irp_here
stop 1
endif
ASSERT (nucl_num > 0)
integer :: align_double
nucl_num_aligned = align_double(nucl_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_charge, (nucl_num) ]
implicit none
BEGIN_DOC
! Nuclear charges
END_DOC
PROVIDE ezfio_filename
nucl_charge = -1.d0
call ezfio_get_nuclei_nucl_charge(nucl_charge)
ASSERT (nucl_charge(:) >= 0.d0)
implicit none
BEGIN_DOC
! Nuclear charges
END_DOC
PROVIDE ezfio_filename
nucl_charge = -1.d0
logical :: has
call ezfio_has_nuclei_nucl_charge(has)
if (has) then
call ezfio_get_nuclei_nucl_charge(nucl_charge)
else
print *, irp_here
stop 1
endif
ASSERT (minval(nucl_charge) >= 0.d0)
END_PROVIDER
BEGIN_PROVIDER [ character*(32), nucl_label, (nucl_num) ]
implicit none
BEGIN_DOC
! Nuclear labels
END_DOC
PROVIDE ezfio_filename
nucl_label = ""
call ezfio_get_nuclei_nucl_label(nucl_label)
implicit none
BEGIN_DOC
! Nuclear labels
END_DOC
PROVIDE ezfio_filename
nucl_label = ""
logical :: has
call ezfio_has_nuclei_nucl_label(has)
if (has) then
call ezfio_get_nuclei_nucl_label(nucl_label)
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_coord, (nucl_num_aligned,3) ]
implicit none
BEGIN_DOC
! Nuclear coordinates in the format (:, {x,y,z})
END_DOC
PROVIDE ezfio_filename
double precision, allocatable :: buffer(:,:)
nucl_coord = 0.d0
allocate (buffer(nucl_num,3))
buffer = 0.d0
call ezfio_get_nuclei_nucl_coord(buffer)
integer :: i,j
do i=1,3
do j=1,nucl_num
nucl_coord(j,i) = buffer(j,i)
implicit none
BEGIN_DOC
! Nuclear coordinates in the format (:, {x,y,z})
END_DOC
PROVIDE ezfio_filename
double precision, allocatable :: buffer(:,:)
nucl_coord = 0.d0
allocate (buffer(nucl_num,3))
buffer = 0.d0
logical :: has
call ezfio_has_nuclei_nucl_coord(has)
if (.not.has) then
print *, irp_here
stop 1
endif
call ezfio_get_nuclei_nucl_coord(buffer)
integer :: i,j
do i=1,3
do j=1,nucl_num
nucl_coord(j,i) = buffer(j,i)
enddo
enddo
enddo
deallocate(buffer)
deallocate(buffer)
character*(64), parameter :: f = '(A16, 4(X,F12.6))'
character*(64), parameter :: ft= '(A16, 4(X,A12 ))'
double precision, parameter :: a0= 0.529177249d0
call write_time(output_Nuclei)
write(output_Nuclei,'(A)') ''
write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)'
write(output_Nuclei,'(A)') '==============================='
write(output_Nuclei,'(A)') ''
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,*) &
' Atom Charge X Y Z '
write(output_Nuclei,ft) &
'================','============','============','============','============'
do i=1,nucl_num
write(output_Nuclei,f) nucl_label(i), nucl_charge(i), &
nucl_coord(i,1)*a0, &
nucl_coord(i,2)*a0, &
nucl_coord(i,3)*a0
enddo
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,'(A)') ''
END_PROVIDER
character*(64), parameter :: f = '(A16, 4(X,F12.6))'
character*(64), parameter :: ft= '(A16, 4(X,A12 ))'
double precision, parameter :: a0= 0.529177249d0
call write_time(output_Nuclei)
write(output_Nuclei,'(A)') ''
write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)'
write(output_Nuclei,'(A)') '==============================='
write(output_Nuclei,'(A)') ''
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,*) &
' Atom Charge X Y Z '
write(output_Nuclei,ft) &
'================','============','============','============','============'
do i=1,nucl_num
write(output_Nuclei,f) nucl_label(i), nucl_charge(i), &
nucl_coord(i,1)*a0, &
nucl_coord(i,2)*a0, &
nucl_coord(i,3)*a0
enddo
write(output_Nuclei,ft) &
'================','============','============','============','============'
write(output_Nuclei,'(A)') ''
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_coord_transp, (3,nucl_num) ]
implicit none
BEGIN_DOC
! Transposed array of nucl_coord
END_DOC
integer :: i, k
nucl_coord_transp = 0.d0
do i=1,nucl_num
nucl_coord_transp(1,i) = nucl_coord(i,1)
nucl_coord_transp(2,i) = nucl_coord(i,2)
nucl_coord_transp(3,i) = nucl_coord(i,3)
enddo
implicit none
BEGIN_DOC
! Transposed array of nucl_coord
END_DOC
integer :: i, k
nucl_coord_transp = 0.d0
do i=1,nucl_num
nucl_coord_transp(1,i) = nucl_coord(i,1)
nucl_coord_transp(2,i) = nucl_coord(i,2)
nucl_coord_transp(3,i) = nucl_coord(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, nucl_dist_2, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_x, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_y, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist_vec_z, (nucl_num_aligned,nucl_num) ]
&BEGIN_PROVIDER [ double precision, nucl_dist, (nucl_num_aligned,nucl_num) ]
implicit none
BEGIN_DOC
! nucl_dist : Nucleus-nucleus distances
! nucl_dist_2 : Nucleus-nucleus distances squared
! nucl_dist_vec : Nucleus-nucleus distances vectors
END_DOC
integer :: ie1, ie2, l
integer,save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
nucl_dist = 0.d0
nucl_dist_2 = 0.d0
nucl_dist_vec_x = 0.d0
nucl_dist_vec_y = 0.d0
nucl_dist_vec_z = 0.d0
endif
implicit none
BEGIN_DOC
! nucl_dist : Nucleus-nucleus distances
! nucl_dist_2 : Nucleus-nucleus distances squared
! nucl_dist_vec : Nucleus-nucleus distances vectors
END_DOC
integer :: ie1, ie2, l
integer,save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
nucl_dist = 0.d0
nucl_dist_2 = 0.d0
nucl_dist_vec_x = 0.d0
nucl_dist_vec_y = 0.d0
nucl_dist_vec_z = 0.d0
endif
do ie2 = 1,nucl_num
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do ie1 = 1,nucl_num_aligned
nucl_dist_vec_x(ie1,ie2) = nucl_coord(ie1,1) - nucl_coord(ie2,1)
nucl_dist_vec_y(ie1,ie2) = nucl_coord(ie1,2) - nucl_coord(ie2,2)
nucl_dist_vec_z(ie1,ie2) = nucl_coord(ie1,3) - nucl_coord(ie2,3)
enddo
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do ie1 = 1,nucl_num_aligned
nucl_dist_2(ie1,ie2) = nucl_dist_vec_x(ie1,ie2)*nucl_dist_vec_x(ie1,ie2) +&
nucl_dist_vec_y(ie1,ie2)*nucl_dist_vec_y(ie1,ie2) + &
nucl_dist_vec_z(ie1,ie2)*nucl_dist_vec_z(ie1,ie2)
nucl_dist(ie1,ie2) = sqrt(nucl_dist_2(ie1,ie2))
ASSERT (nucl_dist(ie1,ie2) > 0.d0)
enddo
enddo
END_PROVIDER
do ie2 = 1,nucl_num
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do ie1 = 1,nucl_num_aligned
nucl_dist_vec_x(ie1,ie2) = nucl_coord(ie1,1) - nucl_coord(ie2,1)
nucl_dist_vec_y(ie1,ie2) = nucl_coord(ie1,2) - nucl_coord(ie2,2)
nucl_dist_vec_z(ie1,ie2) = nucl_coord(ie1,3) - nucl_coord(ie2,3)
enddo
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do ie1 = 1,nucl_num_aligned
nucl_dist_2(ie1,ie2) = nucl_dist_vec_x(ie1,ie2)*nucl_dist_vec_x(ie1,ie2) +&
nucl_dist_vec_y(ie1,ie2)*nucl_dist_vec_y(ie1,ie2) + &
nucl_dist_vec_z(ie1,ie2)*nucl_dist_vec_z(ie1,ie2)
nucl_dist(ie1,ie2) = sqrt(nucl_dist_2(ie1,ie2))
ASSERT (nucl_dist(ie1,ie2) > 0.d0)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
implicit none
BEGIN_DOC
! Nuclear repulsion energy
END_DOC
integer :: k,l
double precision :: Z12, r2, x(3)
nuclear_repulsion = 0.d0
do l = 1, nucl_num
do k = 1, nucl_num
if(k /= l) then
Z12 = nucl_charge(k)*nucl_charge(l)
x(1) = nucl_coord(k,1) - nucl_coord(l,1)
x(2) = nucl_coord(k,2) - nucl_coord(l,2)
x(3) = nucl_coord(k,3) - nucl_coord(l,3)
r2 = x(1)*x(1) + x(2)*x(2) + x(3)*x(3)
nuclear_repulsion += Z12/dsqrt(r2)
endif
implicit none
BEGIN_DOC
! Nuclear repulsion energy
END_DOC
integer :: k,l
double precision :: Z12, r2, x(3)
nuclear_repulsion = 0.d0
do l = 1, nucl_num
do k = 1, nucl_num
if(k /= l) then
Z12 = nucl_charge(k)*nucl_charge(l)
x(1) = nucl_coord(k,1) - nucl_coord(l,1)
x(2) = nucl_coord(k,2) - nucl_coord(l,2)
x(3) = nucl_coord(k,3) - nucl_coord(l,3)
r2 = x(1)*x(1) + x(2)*x(2) + x(3)*x(3)
nuclear_repulsion += Z12/dsqrt(r2)
endif
enddo
enddo
enddo
nuclear_repulsion *= 0.5d0
call write_time(output_Nuclei)
call write_double(output_Nuclei,nuclear_repulsion, &
'Nuclear repulsion energy')
nuclear_repulsion *= 0.5d0
call write_time(output_Nuclei)
call write_double(output_Nuclei,nuclear_repulsion, &
'Nuclear repulsion energy')
END_PROVIDER

31
src/Output/README.rst
View File

@ -25,3 +25,34 @@ All output should be printed using routines present in this module.
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`output_cpu_time_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Output/output.irp.f#L2>`_
Initial CPU and wall times when printing in the output files
`output_wall_time_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Output/output.irp.f#L1>`_
Initial CPU and wall times when printing in the output files
`write_double <http://github.com/LCPQ/quantum_package/tree/master/src/Output/output.irp.f#L49>`_
Write a double precision value in output
`write_int <http://github.com/LCPQ/quantum_package/tree/master/src/Output/output.irp.f#L64>`_
Write an integer value in output
`write_time <http://github.com/LCPQ/quantum_package/tree/master/src/Output/output.irp.f#L33>`_
Write a time stamp in the output for chronological reconstruction

8
src/README.rst
View File

@ -122,5 +122,11 @@ Needed Modules
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `Hartree-fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree-fock>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOGuess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `DensityMatrix <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix>`_
* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_

153
src/Utils/README.rst
View File

@ -4,3 +4,156 @@ Utils Module
Contains general purpose utilities.
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`apply_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L146>`_
Apply the rotation found by find_rotation
`find_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L127>`_
Find A.C = B
`get_pseudo_inverse <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L73>`_
Find C = A^-1
`lapack_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L158>`_
Diagonalize matrix H
`ortho_lowdin <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L1>`_
Compute U.S^-1/2 canonical orthogonalization
`add_poly <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L243>`_
Add two polynomials
D(t) =! D(t) +( B(t)+C(t))
`add_poly_multiply <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L271>`_
Add a polynomial multiplied by a constant
D(t) =! D(t) +( cst * B(t))
`f_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L345>`_
function that calculates the following integral
\int_{\-infty}^{+\infty} x^n \exp(-p x^2) dx
`gaussian_product <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L121>`_
Gaussian product in 1D.
e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
`gaussian_product_x <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L163>`_
Gaussian product in 1D.
e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
`give_explicit_poly_and_gaussian <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L46>`_
Transforms the product of
(x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
into
fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
* [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
* [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
`give_explicit_poly_and_gaussian_x <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L1>`_
Transform the product of
(x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
into
fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2)
`hermite <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L468>`_
Hermite polynomial
`multiply_poly <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L201>`_
Multiply two polynomials
D(t) =! D(t) +( B(t)*C(t))
`recentered_poly2 <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L300>`_
Recenter two polynomials
`rint <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L373>`_
.. math::
.br
\int_0^1 dx \exp(-p x^2) x^n
.br
`rint1 <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L524>`_
Standard version of rint
`rint_large_n <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L493>`_
Version of rint for large values of n
`rint_sum <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L417>`_
Needed for the calculation of two-electron integrals.
`overlap_gaussian_x <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/one_e_integration.irp.f#L1>`_
.. math::
.br
\sum_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha(x-A_x)^2) exp(-beta(x-B_X)^2) dx
.br
`overlap_gaussian_xyz <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/one_e_integration.irp.f#L37>`_
.. math::
.br
S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx \\
S = S_x S_y S_z
.br
`overlap_x_abs <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/one_e_integration.irp.f#L99>`_
.. math ::
.br
\int_{-infty}^{+infty} (x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) dx
.br
`align_double <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L65>`_
Compute 1st dimension such that it is aligned for vectorization.
`all_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L1>`_
Dummy provider to provide all utils
`binom <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L47>`_
Binomial coefficients
`binom_func <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L16>`_
.. math ::
.br
\frac{i!}{j!(i-j)!}
.br
`binom_transp <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L48>`_
Binomial coefficients
`dble_fact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L124>`_
n!!
`fact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L80>`_
n!
`fact_inv <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L112>`_
1/n!
`inv_int <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L171>`_
1/i
`normalize <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L272>`_
Normalizes vector u
u is expected to be aligned in memory.
`nproc <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L197>`_
Number of current OpenMP threads
`u_dot_u <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L244>`_
Compute <u|u>
u is expected to be aligned in memory.
`u_dot_v <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L213>`_
Compute <u|v>
u and v are expected to be aligned in memory.
`wall_time <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L182>`_
The equivalent of cpu_time, but for the wall time.
`write_git_log <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L157>`_
Write the last git commit in file iunit.

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

@ -382,7 +382,8 @@ BEGIN_TEMPLATE
return
endif
ASSERT (iradix > 0)
ASSERT (iradix >= 0)
if (isize < 48) then
call insertion_$Xsort$big(x,iorder,isize)
return

104
src/Utils/util.irp.f
View File

@ -1,16 +1,16 @@
BEGIN_PROVIDER [ logical, all_utils ]
implicit none
BEGIN_DOC
! Dummy provider to provide all utils
END_DOC
! Do not move this : it greps itself
BEGIN_SHELL [ /bin/bash ]
implicit none
BEGIN_DOC
! Dummy provider to provide all utils
END_DOC
! Do not move this : it greps itself
BEGIN_SHELL [ /bin/bash ]
for i in $(grep "BEGIN_PROVIDER" $QPACKAGE_ROOT/src/Utils/*.irp.f | cut -d ',' -f 2 | cut -d ']' -f 1 | tail --lines=+3 )
do
echo PROVIDE $i
echo PROVIDE $i
done
END_SHELL
END_SHELL
END_PROVIDER
double precision function binom_func(i,j)
@ -209,3 +209,89 @@ BEGIN_PROVIDER [ integer, nproc ]
!$OMP END PARALLEL
END_PROVIDER
double precision function u_dot_v(u,v,sze)
implicit none
BEGIN_DOC
! Compute <u|v>
! u and v are expected to be aligned in memory.
END_DOC
integer, intent(in) :: sze
double precision, intent(in) :: u(sze),v(sze)
integer :: i,t1, t2, t3, t4
ASSERT (sze > 0)
t1 = 0
t2 = sze/4
t3 = t2+t2
t4 = t3+t2
u_dot_v = 0.d0
!DIR$ VECTOR ALWAYS
!DIR$ VECTOR ALIGNED
do i=1,t2
u_dot_v = u_dot_v + u(t1+i)*v(t1+i) + u(t2+i)*v(t2+i) + &
u(t3+i)*v(t3+i) + u(t4+i)*v(t4+i)
enddo
!DIR$ VECTOR ALWAYS
!DIR$ VECTOR ALIGNED
do i=t4+t2+1,sze
u_dot_v = u_dot_v + u(i)*v(i)
enddo
end
double precision function u_dot_u(u,sze)
implicit none
BEGIN_DOC
! Compute <u|u>
! u is expected to be aligned in memory.
END_DOC
integer, intent(in) :: sze
double precision, intent(in) :: u(sze)
integer :: i
integer :: t1, t2, t3, t4
ASSERT (sze > 0)
t1 = 0
t2 = sze/4
t3 = t2+t2
t4 = t3+t2
u_dot_u = 0.d0
do i=1,t2
u_dot_u = u_dot_u + u(t1+i)*u(t1+i) + u(t2+i)*u(t2+i) + &
u(t3+i)*u(t3+i) + u(t4+i)*u(t4+i)
enddo
do i=t4+t2+1,sze
u_dot_u = u_dot_u+u(i)*u(i)
enddo
end
subroutine normalize(u,sze)
implicit none
BEGIN_DOC
! Normalizes vector u
! u is expected to be aligned in memory.
END_DOC
integer, intent(in) :: sze
double precision, intent(inout):: u(sze)
double precision :: d
double precision, external :: u_dot_u
integer :: i
!DIR$ FORCEINLINE
d = 1.d0/dsqrt( u_dot_u(u,sze) )
if (d /= 1.d0) then
do i=1,sze
u(i) = d*u(i)
enddo
endif
end