mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-10-10 01:43:06 +02:00
Merge branch 'dev-stable-tc-scf' of https://github.com/AbdAmmar/qp2 into dev-stable-tc-scf
This commit is contained in:
AbdAmmar
commit
45f7d69e70
2
configure
vendored
2
configure
vendored
@ -215,7 +215,6 @@ EOF
|
||||
cd trexio-${VERSION}
|
||||
./configure --prefix=\${QP_ROOT} --without-hdf5
|
||||
make -j 8 && make -j 8 check && make -j 8 install
|
||||
cp ${QP_ROOT}/include/trexio_f.f90 ${QP_ROOT}/src/ezfio_files
|
||||
tar -zxvf "\${QP_ROOT}"/external/qp2-dependencies/${ARCHITECTURE}/ninja.tar.gz
|
||||
mv ninja "\${QP_ROOT}"/bin/
|
||||
EOF
|
||||
@ -229,7 +228,6 @@ EOF
|
||||
cd trexio-${VERSION}
|
||||
./configure --prefix=\${QP_ROOT}
|
||||
make -j 8 && make -j 8 check && make -j 8 install
|
||||
cp ${QP_ROOT}/include/trexio_f.f90 ${QP_ROOT}/src/ezfio_files
|
||||
EOF
|
||||
|
||||
|
||||
|
||||
@ -44,8 +44,12 @@ end = struct
|
||||
let get_default = Qpackage.get_ezfio_default "ao_basis";;
|
||||
|
||||
let read_ao_basis () =
|
||||
Ezfio.get_ao_basis_ao_basis ()
|
||||
|> AO_basis_name.of_string
|
||||
let result =
|
||||
Ezfio.get_ao_basis_ao_basis ()
|
||||
in
|
||||
if result <> "None" then
|
||||
AO_basis_name.of_string result
|
||||
else failwith "No basis"
|
||||
;;
|
||||
|
||||
let read_ao_num () =
|
||||
@ -192,7 +196,7 @@ end = struct
|
||||
ao_expo ;
|
||||
ao_cartesian ;
|
||||
ao_normalized ;
|
||||
primitives_normalized ;
|
||||
primitives_normalized ;
|
||||
} = b
|
||||
in
|
||||
write_md5 b ;
|
||||
@ -207,7 +211,7 @@ end = struct
|
||||
Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
|
||||
~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
|
||||
|
||||
let ao_nucl =
|
||||
let ao_nucl =
|
||||
Array.to_list ao_nucl
|
||||
|> list_map Nucl_number.to_int
|
||||
in
|
||||
@ -215,7 +219,7 @@ end = struct
|
||||
~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ;
|
||||
|
||||
let ao_power =
|
||||
let l = Array.to_list ao_power in
|
||||
let l = Array.to_list ao_power in
|
||||
List.concat [
|
||||
(list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.x) l) ;
|
||||
(list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.y) l) ;
|
||||
@ -227,7 +231,7 @@ end = struct
|
||||
Ezfio.set_ao_basis_ao_cartesian(ao_cartesian);
|
||||
Ezfio.set_ao_basis_ao_normalized(ao_normalized);
|
||||
Ezfio.set_ao_basis_primitives_normalized(primitives_normalized);
|
||||
|
||||
|
||||
let ao_coef =
|
||||
Array.to_list ao_coef
|
||||
|> list_map AO_coef.to_float
|
||||
@ -267,7 +271,10 @@ end = struct
|
||||
|> Ezfio.set_ao_basis_ao_md5 ;
|
||||
Some result
|
||||
with
|
||||
| _ -> (Ezfio.set_ao_basis_ao_md5 "None" ; None)
|
||||
| _ -> ( "None"
|
||||
|> Digest.string
|
||||
|> Digest.to_hex
|
||||
|> Ezfio.set_ao_basis_ao_md5 ; None)
|
||||
;;
|
||||
|
||||
|
||||
@ -276,7 +283,7 @@ end = struct
|
||||
to_basis b
|
||||
|> Long_basis.of_basis
|
||||
|> Array.of_list
|
||||
and unordered_basis =
|
||||
and unordered_basis =
|
||||
to_long_basis b
|
||||
|> Array.of_list
|
||||
in
|
||||
@ -289,15 +296,15 @@ end = struct
|
||||
(a.(i) <- None ; i)
|
||||
else
|
||||
find x a (i+1)
|
||||
and find2 (s,g,n) a i =
|
||||
and find2 (s,g,n) a i =
|
||||
if i = Array.length a then -1
|
||||
else
|
||||
match a.(i) with
|
||||
match a.(i) with
|
||||
| None -> find2 (s,g,n) a (i+1)
|
||||
| Some (s', g', n') ->
|
||||
if s <> s' || n <> n' then find2 (s,g,n) a (i+1)
|
||||
else
|
||||
let lc = list_map (fun (prim, _) -> prim) g.Gto.lc
|
||||
let lc = list_map (fun (prim, _) -> prim) g.Gto.lc
|
||||
and lc' = list_map (fun (prim, _) -> prim) g'.Gto.lc
|
||||
in
|
||||
if lc <> lc' then find2 (s,g,n) a (i+1) else (a.(i) <- None ; i)
|
||||
@ -313,13 +320,13 @@ end = struct
|
||||
let ao_num = List.length long_basis |> AO_number.of_int in
|
||||
let ao_prim_num =
|
||||
list_map (fun (_,g,_) -> List.length g.Gto.lc
|
||||
|> AO_prim_number.of_int ) long_basis
|
||||
|> AO_prim_number.of_int ) long_basis
|
||||
|> Array.of_list
|
||||
and ao_nucl =
|
||||
list_map (fun (_,_,n) -> n) long_basis
|
||||
list_map (fun (_,_,n) -> n) long_basis
|
||||
|> Array.of_list
|
||||
and ao_power =
|
||||
list_map (fun (x,_,_) -> x) long_basis
|
||||
list_map (fun (x,_,_) -> x) long_basis
|
||||
|> Array.of_list
|
||||
in
|
||||
let ao_prim_num_max = Array.fold_left (fun s x ->
|
||||
@ -329,16 +336,16 @@ end = struct
|
||||
in
|
||||
|
||||
let gtos =
|
||||
list_map (fun (_,x,_) -> x) long_basis
|
||||
list_map (fun (_,x,_) -> x) long_basis
|
||||
in
|
||||
let create_expo_coef ec =
|
||||
let coefs =
|
||||
begin match ec with
|
||||
| `Coefs -> list_map (fun x->
|
||||
list_map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos
|
||||
list_map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos
|
||||
| `Expos -> list_map (fun x->
|
||||
list_map (fun (prim,_) -> AO_expo.to_float
|
||||
prim.GaussianPrimitive.expo) x.Gto.lc ) gtos
|
||||
prim.GaussianPrimitive.expo) x.Gto.lc ) gtos
|
||||
end
|
||||
in
|
||||
let rec get_n n accu = function
|
||||
@ -360,7 +367,7 @@ end = struct
|
||||
let ao_coef = create_expo_coef `Coefs
|
||||
|> Array.of_list
|
||||
|> Array.map AO_coef.of_float
|
||||
and ao_expo = create_expo_coef `Expos
|
||||
and ao_expo = create_expo_coef `Expos
|
||||
|> Array.of_list
|
||||
|> Array.map AO_expo.of_float
|
||||
in
|
||||
@ -372,7 +379,7 @@ end = struct
|
||||
}
|
||||
;;
|
||||
|
||||
let reorder b =
|
||||
let reorder b =
|
||||
let order = ordering b in
|
||||
let f a = Array.init (Array.length a) (fun i -> a.(order.(i))) in
|
||||
let ao_prim_num_max = AO_prim_number.to_int b.ao_prim_num_max
|
||||
@ -464,7 +471,7 @@ Basis set (read-only) ::
|
||||
| line :: tail ->
|
||||
let line = String.trim line in
|
||||
if line = "Basis set (read-only) ::" then
|
||||
String.concat "\n" tail
|
||||
String.concat "\n" tail
|
||||
else
|
||||
extract_basis tail
|
||||
in
|
||||
|
||||
@ -56,7 +56,10 @@ end = struct
|
||||
let read_ao_md5 () =
|
||||
let ao_md5 =
|
||||
match (Input_ao_basis.Ao_basis.read ()) with
|
||||
| None -> failwith "Unable to read AO basis"
|
||||
| None -> ("None"
|
||||
|> Digest.string
|
||||
|> Digest.to_hex
|
||||
|> MD5.of_string)
|
||||
| Some result -> Input_ao_basis.Ao_basis.to_md5 result
|
||||
in
|
||||
let result =
|
||||
|
||||
@ -38,7 +38,8 @@ let run slave ?prefix exe ezfio_file =
|
||||
| Unix.Unix_error _ -> try_new_port (port_number+100)
|
||||
in
|
||||
let result =
|
||||
try_new_port 41279
|
||||
let port = 10*(Unix.getpid () mod 2823) + 32_769 in
|
||||
try_new_port port
|
||||
in
|
||||
Zmq.Socket.close dummy_socket;
|
||||
Zmq.Context.terminate zmq_context;
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
#!/usr/bin/env python3
|
||||
"""
|
||||
Save the .o from a .f90
|
||||
and is the .o is asked a second time, retur it
|
||||
and is the .o is asked a second time, return it
|
||||
Take in argv command like:
|
||||
ifort -g -openmp -I IRPF90_temp/Ezfio_files/ -c IRPF90_temp/Integrals_Monoelec/kin_ao_ints.irp.module.F90 -o IRPF90_temp/Integrals_Monoelec/kin_ao_ints.irp.module.o
|
||||
"""
|
||||
|
||||
@ -13,11 +13,17 @@ Options:
|
||||
|
||||
import sys
|
||||
import os
|
||||
import trexio
|
||||
import numpy as np
|
||||
from functools import reduce
|
||||
from ezfio import ezfio
|
||||
from docopt import docopt
|
||||
import qp_bitmasks
|
||||
|
||||
try:
|
||||
import trexio
|
||||
except ImportError:
|
||||
print("Error: trexio python module is not found. Try python3 -m pip install trexio")
|
||||
sys.exit(1)
|
||||
|
||||
|
||||
try:
|
||||
@ -90,14 +96,15 @@ def write_ezfio(trexio_filename, filename):
|
||||
p = re.compile(r'(\d*)$')
|
||||
label = [p.sub("", x).capitalize() for x in label]
|
||||
ezfio.set_nuclei_nucl_label(label)
|
||||
print("OK")
|
||||
|
||||
else:
|
||||
ezfio.set_nuclei_nucl_num(1)
|
||||
ezfio.set_nuclei_nucl_charge([0.])
|
||||
ezfio.set_nuclei_nucl_coord([0.,0.,0.])
|
||||
ezfio.set_nuclei_nucl_label(["X"])
|
||||
print("None")
|
||||
|
||||
print("OK")
|
||||
|
||||
|
||||
print("Electrons\t...\t", end=' ')
|
||||
@ -105,12 +112,12 @@ def write_ezfio(trexio_filename, filename):
|
||||
try:
|
||||
num_beta = trexio.read_electron_dn_num(trexio_file)
|
||||
except:
|
||||
num_beta = sum(charge)//2
|
||||
num_beta = int(sum(charge))//2
|
||||
|
||||
try:
|
||||
num_alpha = trexio.read_electron_up_num(trexio_file)
|
||||
except:
|
||||
num_alpha = sum(charge) - num_beta
|
||||
num_alpha = int(sum(charge)) - num_beta
|
||||
|
||||
if num_alpha == 0:
|
||||
print("\n\nError: There are zero electrons in the TREXIO file.\n\n")
|
||||
@ -118,7 +125,7 @@ def write_ezfio(trexio_filename, filename):
|
||||
ezfio.set_electrons_elec_alpha_num(num_alpha)
|
||||
ezfio.set_electrons_elec_beta_num(num_beta)
|
||||
|
||||
print("OK")
|
||||
print(f"{num_alpha} {num_beta}")
|
||||
|
||||
print("Basis\t\t...\t", end=' ')
|
||||
|
||||
@ -126,60 +133,113 @@ def write_ezfio(trexio_filename, filename):
|
||||
try:
|
||||
basis_type = trexio.read_basis_type(trexio_file)
|
||||
|
||||
if basis_type.lower() not in ["gaussian", "slater"]:
|
||||
raise TypeError
|
||||
if basis_type.lower() in ["gaussian", "slater"]:
|
||||
shell_num = trexio.read_basis_shell_num(trexio_file)
|
||||
prim_num = trexio.read_basis_prim_num(trexio_file)
|
||||
ang_mom = trexio.read_basis_shell_ang_mom(trexio_file)
|
||||
nucl_index = trexio.read_basis_nucleus_index(trexio_file)
|
||||
exponent = trexio.read_basis_exponent(trexio_file)
|
||||
coefficient = trexio.read_basis_coefficient(trexio_file)
|
||||
shell_index = trexio.read_basis_shell_index(trexio_file)
|
||||
ao_shell = trexio.read_ao_shell(trexio_file)
|
||||
|
||||
shell_num = trexio.read_basis_shell_num(trexio_file)
|
||||
prim_num = trexio.read_basis_prim_num(trexio_file)
|
||||
ang_mom = trexio.read_basis_shell_ang_mom(trexio_file)
|
||||
nucl_index = trexio.read_basis_nucleus_index(trexio_file)
|
||||
exponent = trexio.read_basis_exponent(trexio_file)
|
||||
coefficient = trexio.read_basis_coefficient(trexio_file)
|
||||
shell_index = trexio.read_basis_shell_index(trexio_file)
|
||||
ao_shell = trexio.read_ao_shell(trexio_file)
|
||||
ezfio.set_basis_basis("Read from TREXIO")
|
||||
ezfio.set_ao_basis_ao_basis("Read from TREXIO")
|
||||
ezfio.set_basis_shell_num(shell_num)
|
||||
ezfio.set_basis_prim_num(prim_num)
|
||||
ezfio.set_basis_shell_ang_mom(ang_mom)
|
||||
ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
|
||||
ezfio.set_basis_prim_expo(exponent)
|
||||
ezfio.set_basis_prim_coef(coefficient)
|
||||
|
||||
ezfio.set_basis_basis("Read from TREXIO")
|
||||
ezfio.set_basis_shell_num(shell_num)
|
||||
ezfio.set_basis_prim_num(prim_num)
|
||||
ezfio.set_basis_shell_ang_mom(ang_mom)
|
||||
ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
|
||||
ezfio.set_basis_prim_expo(exponent)
|
||||
ezfio.set_basis_prim_coef(coefficient)
|
||||
nucl_shell_num = []
|
||||
prev = None
|
||||
m = 0
|
||||
for i in ao_shell:
|
||||
if i != prev:
|
||||
m += 1
|
||||
if prev is None or nucl_index[i] != nucl_index[prev]:
|
||||
nucl_shell_num.append(m)
|
||||
m = 0
|
||||
prev = i
|
||||
assert (len(nucl_shell_num) == nucl_num)
|
||||
|
||||
nucl_shell_num = []
|
||||
prev = None
|
||||
m = 0
|
||||
for i in ao_shell:
|
||||
if i != prev:
|
||||
m += 1
|
||||
if prev is None or nucl_index[i] != nucl_index[prev]:
|
||||
nucl_shell_num.append(m)
|
||||
m = 0
|
||||
prev = i
|
||||
assert (len(nucl_shell_num) == nucl_num)
|
||||
shell_prim_num = []
|
||||
prev = shell_index[0]
|
||||
count = 0
|
||||
for i in shell_index:
|
||||
if i != prev:
|
||||
shell_prim_num.append(count)
|
||||
count = 0
|
||||
count += 1
|
||||
prev = i
|
||||
shell_prim_num.append(count)
|
||||
|
||||
shell_prim_num = []
|
||||
prev = shell_index[0]
|
||||
count = 0
|
||||
for i in shell_index:
|
||||
if i != prev:
|
||||
shell_prim_num.append(count)
|
||||
count = 0
|
||||
count += 1
|
||||
prev = i
|
||||
shell_prim_num.append(count)
|
||||
assert (len(shell_prim_num) == shell_num)
|
||||
|
||||
assert (len(shell_prim_num) == shell_num)
|
||||
|
||||
ezfio.set_basis_shell_prim_num(shell_prim_num)
|
||||
ezfio.set_basis_shell_index([x+1 for x in shell_index])
|
||||
ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
|
||||
ezfio.set_basis_shell_prim_num(shell_prim_num)
|
||||
ezfio.set_basis_shell_index([x+1 for x in shell_index])
|
||||
ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
|
||||
|
||||
|
||||
shell_factor = trexio.read_basis_shell_factor(trexio_file)
|
||||
prim_factor = trexio.read_basis_prim_factor(trexio_file)
|
||||
shell_factor = trexio.read_basis_shell_factor(trexio_file)
|
||||
prim_factor = trexio.read_basis_prim_factor(trexio_file)
|
||||
|
||||
print("OK")
|
||||
elif basis_type.lower() == "numerical":
|
||||
|
||||
shell_num = trexio.read_basis_shell_num(trexio_file)
|
||||
prim_num = shell_num
|
||||
ang_mom = trexio.read_basis_shell_ang_mom(trexio_file)
|
||||
nucl_index = trexio.read_basis_nucleus_index(trexio_file)
|
||||
exponent = [1.]*prim_num
|
||||
coefficient = [1.]*prim_num
|
||||
shell_index = [i for i in range(shell_num)]
|
||||
ao_shell = trexio.read_ao_shell(trexio_file)
|
||||
|
||||
ezfio.set_basis_basis("None")
|
||||
ezfio.set_ao_basis_ao_basis("None")
|
||||
ezfio.set_basis_shell_num(shell_num)
|
||||
ezfio.set_basis_prim_num(prim_num)
|
||||
ezfio.set_basis_shell_ang_mom(ang_mom)
|
||||
ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
|
||||
ezfio.set_basis_prim_expo(exponent)
|
||||
ezfio.set_basis_prim_coef(coefficient)
|
||||
|
||||
nucl_shell_num = []
|
||||
prev = None
|
||||
m = 0
|
||||
for i in ao_shell:
|
||||
if i != prev:
|
||||
m += 1
|
||||
if prev is None or nucl_index[i] != nucl_index[prev]:
|
||||
nucl_shell_num.append(m)
|
||||
m = 0
|
||||
prev = i
|
||||
assert (len(nucl_shell_num) == nucl_num)
|
||||
|
||||
shell_prim_num = []
|
||||
prev = shell_index[0]
|
||||
count = 0
|
||||
for i in shell_index:
|
||||
if i != prev:
|
||||
shell_prim_num.append(count)
|
||||
count = 0
|
||||
count += 1
|
||||
prev = i
|
||||
shell_prim_num.append(count)
|
||||
|
||||
assert (len(shell_prim_num) == shell_num)
|
||||
|
||||
ezfio.set_basis_shell_prim_num(shell_prim_num)
|
||||
ezfio.set_basis_shell_index([x+1 for x in shell_index])
|
||||
ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
|
||||
|
||||
shell_factor = trexio.read_basis_shell_factor(trexio_file)
|
||||
prim_factor = [1.]*prim_num
|
||||
else:
|
||||
raise TypeError
|
||||
|
||||
print(basis_type)
|
||||
except:
|
||||
print("None")
|
||||
ezfio.set_ao_basis_ao_cartesian(True)
|
||||
@ -256,9 +316,11 @@ def write_ezfio(trexio_filename, filename):
|
||||
# ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
|
||||
ezfio.set_ao_basis_ao_coef(coef)
|
||||
ezfio.set_ao_basis_ao_expo(expo)
|
||||
ezfio.set_ao_basis_ao_basis("Read from TREXIO")
|
||||
|
||||
print("OK")
|
||||
print("OK")
|
||||
|
||||
else:
|
||||
print("None")
|
||||
|
||||
|
||||
# _
|
||||
@ -279,6 +341,7 @@ def write_ezfio(trexio_filename, filename):
|
||||
except:
|
||||
label = "None"
|
||||
ezfio.set_mo_basis_mo_label(label)
|
||||
ezfio.set_determinants_mo_label(label)
|
||||
|
||||
try:
|
||||
clss = trexio.read_mo_class(trexio_file)
|
||||
@ -303,10 +366,10 @@ def write_ezfio(trexio_filename, filename):
|
||||
for i in range(num_beta):
|
||||
mo_occ[i] += 1.
|
||||
ezfio.set_mo_basis_mo_occ(mo_occ)
|
||||
print("OK")
|
||||
except:
|
||||
pass
|
||||
print("None")
|
||||
|
||||
print("OK")
|
||||
|
||||
|
||||
print("Pseudos\t\t...\t", end=' ')
|
||||
@ -386,8 +449,23 @@ def write_ezfio(trexio_filename, filename):
|
||||
ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl)
|
||||
ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl)
|
||||
ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl)
|
||||
print("OK")
|
||||
|
||||
else:
|
||||
print("None")
|
||||
|
||||
print("Determinant\t\t...\t", end=' ')
|
||||
alpha = [ i for i in range(num_alpha) ]
|
||||
beta = [ i for i in range(num_beta) ]
|
||||
if trexio.has_mo_spin(trexio_file):
|
||||
spin = trexio.read_mo_spin(trexio_file)
|
||||
beta = [ i for i in range(mo_num) if spin[i] == 1 ]
|
||||
beta = [ beta[i] for i in range(num_beta) ]
|
||||
|
||||
alpha = qp_bitmasks.BitMask(alpha)
|
||||
beta = qp_bitmasks.BitMask(beta )
|
||||
print(alpha)
|
||||
print(beta)
|
||||
print("OK")
|
||||
|
||||
|
||||
@ -22,7 +22,7 @@ def int_to_string(s):
|
||||
assert s>=0
|
||||
AssertionError
|
||||
"""
|
||||
assert type(s) in (int, long)
|
||||
assert type(s) == int
|
||||
assert s>=0
|
||||
return '{s:0b}'.format(s=s)
|
||||
|
||||
@ -62,7 +62,7 @@ def int_to_bitmask(s,bit_kind_size=BIT_KIND_SIZE):
|
||||
['1111111111111111111111111111111111111111111111111111111111110110']
|
||||
>>>
|
||||
"""
|
||||
assert type(s) in (int, long)
|
||||
assert type(s) == int
|
||||
if s < 0:
|
||||
s = s + (1 << bit_kind_size)
|
||||
return ['{s:0{width}b}'.format(s=s,width=bit_kind_size)]
|
||||
@ -104,7 +104,7 @@ class BitMask(object):
|
||||
return self._data_int[i]
|
||||
|
||||
def __setitem__(self,i,value):
|
||||
if type(value) in (int,long):
|
||||
if type(value) == int :
|
||||
self._data_int[i] = value
|
||||
elif type(value) == str:
|
||||
s = string_to_bitmask(value,bit_kind_size=self.bit_kind_size)[0]
|
||||
|
||||
@ -67,3 +67,15 @@ doc: Use normalized primitive functions
|
||||
interface: ezfio, provider
|
||||
default: true
|
||||
|
||||
[ao_expoim_cosgtos]
|
||||
type: double precision
|
||||
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
|
||||
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
|
||||
interface: ezfio, provider
|
||||
|
||||
[use_cosgtos]
|
||||
type: logical
|
||||
doc: If true, use cosgtos for AO integrals
|
||||
interface: ezfio
|
||||
default: False
|
||||
|
||||
|
||||
34
src/ao_basis/cosgtos.irp.f
Normal file
34
src/ao_basis/cosgtos.irp.f
Normal file
@ -0,0 +1,34 @@
|
||||
BEGIN_PROVIDER [ logical, use_cosgtos ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! If true, use cosgtos for AO integrals
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
use_cosgtos = .False.
|
||||
if (mpi_master) then
|
||||
call ezfio_has_ao_basis_use_cosgtos(has)
|
||||
if (has) then
|
||||
! write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..'
|
||||
call ezfio_get_ao_basis_use_cosgtos(use_cosgtos)
|
||||
else
|
||||
call ezfio_set_ao_basis_use_cosgtos(use_cosgtos)
|
||||
endif
|
||||
endif
|
||||
IRP_IF MPI_DEBUG
|
||||
print *, irp_here, mpi_rank
|
||||
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
|
||||
IRP_ENDIF
|
||||
IRP_IF MPI
|
||||
include 'mpif.h'
|
||||
integer :: ierr
|
||||
call MPI_BCAST( use_cosgtos, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
|
||||
if (ierr /= MPI_SUCCESS) then
|
||||
stop 'Unable to read use_cosgtos with MPI'
|
||||
endif
|
||||
IRP_ENDIF
|
||||
|
||||
! call write_time(6)
|
||||
|
||||
END_PROVIDER
|
||||
@ -62,6 +62,7 @@ END_PROVIDER
|
||||
double precision :: tmp_cent_x, tmp_cent_y, tmp_cent_z
|
||||
|
||||
provide j1b_pen
|
||||
provide j1b_pen_coef
|
||||
|
||||
List_all_comb_b2_coef = 0.d0
|
||||
List_all_comb_b2_expo = 0.d0
|
||||
@ -127,8 +128,8 @@ END_PROVIDER
|
||||
List_all_comb_b2_expo( 1) = 0.d0
|
||||
List_all_comb_b2_cent(1:3,1) = 0.d0
|
||||
do i = 1, nucl_num
|
||||
List_all_comb_b2_coef( i+1) = -1.d0
|
||||
List_all_comb_b2_expo( i+1) = j1b_pen( i)
|
||||
List_all_comb_b2_coef( i+1) = -1.d0 * j1b_pen_coef(i)
|
||||
List_all_comb_b2_expo( i+1) = j1b_pen(i)
|
||||
List_all_comb_b2_cent(1,i+1) = nucl_coord(i,1)
|
||||
List_all_comb_b2_cent(2,i+1) = nucl_coord(i,2)
|
||||
List_all_comb_b2_cent(3,i+1) = nucl_coord(i,3)
|
||||
@ -225,6 +226,7 @@ END_PROVIDER
|
||||
double precision :: dx, dy, dz, r2
|
||||
|
||||
provide j1b_pen
|
||||
provide j1b_pen_coef
|
||||
|
||||
List_all_comb_b3_coef = 0.d0
|
||||
List_all_comb_b3_expo = 0.d0
|
||||
@ -296,8 +298,8 @@ END_PROVIDER
|
||||
|
||||
do i = 1, nucl_num
|
||||
ii = ii + 1
|
||||
List_all_comb_b3_coef( ii) = -2.d0
|
||||
List_all_comb_b3_expo( ii) = j1b_pen( i)
|
||||
List_all_comb_b3_coef( ii) = -2.d0 * j1b_pen_coef(i)
|
||||
List_all_comb_b3_expo( ii) = j1b_pen(i)
|
||||
List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
|
||||
List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
|
||||
List_all_comb_b3_cent(3,ii) = nucl_coord(i,3)
|
||||
@ -305,7 +307,7 @@ END_PROVIDER
|
||||
|
||||
do i = 1, nucl_num
|
||||
ii = ii + 1
|
||||
List_all_comb_b3_coef( ii) = 1.d0
|
||||
List_all_comb_b3_coef( ii) = 1.d0 * j1b_pen_coef(i) * j1b_pen_coef(i)
|
||||
List_all_comb_b3_expo( ii) = 2.d0 * j1b_pen(i)
|
||||
List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
|
||||
List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
|
||||
@ -337,7 +339,7 @@ END_PROVIDER
|
||||
|
||||
ii = ii + 1
|
||||
! x 2 to avoid doing integrals twice
|
||||
List_all_comb_b3_coef( ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2)
|
||||
List_all_comb_b3_coef( ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2) * j1b_pen_coef(i) * j1b_pen_coef(j)
|
||||
List_all_comb_b3_expo( ii) = tmp3
|
||||
List_all_comb_b3_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj)
|
||||
List_all_comb_b3_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj)
|
||||
|
||||
@ -1,3 +1,2 @@
|
||||
ao_basis
|
||||
pseudo
|
||||
cosgtos_ao_int
|
||||
|
||||
@ -104,6 +104,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
|
||||
IF(do_pseudo) THEN
|
||||
ao_integrals_n_e += ao_pseudo_integrals
|
||||
ENDIF
|
||||
IF(point_charges) THEN
|
||||
ao_integrals_n_e += ao_integrals_pt_chrg
|
||||
ENDIF
|
||||
|
||||
endif
|
||||
|
||||
@ -455,10 +458,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(c)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
ny=0
|
||||
call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in)
|
||||
call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
! call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,R1x,d,nd)
|
||||
else
|
||||
do ix=0,n_pt_in
|
||||
X(ix) = 0.d0
|
||||
@ -469,7 +474,8 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(a-1)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
|
||||
nx = nd
|
||||
do ix=0,n_pt_in
|
||||
@ -479,10 +485,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(c)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
ny=0
|
||||
call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in)
|
||||
call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
! call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,R1x,d,nd)
|
||||
endif
|
||||
end
|
||||
|
||||
@ -519,7 +527,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(c-1)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
ny = 0
|
||||
do ix=0,dim
|
||||
Y(ix) = 0.d0
|
||||
@ -527,7 +536,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
|
||||
|
||||
call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim)
|
||||
if(ny.ge.0)then
|
||||
call multiply_poly(Y,ny,R1xp,2,d,nd)
|
||||
! call multiply_poly(Y,ny,R1xp,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,R1xp,d,nd)
|
||||
endif
|
||||
endif
|
||||
end
|
||||
|
||||
@ -4,6 +4,19 @@ doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: None
|
||||
|
||||
[ao_integrals_threshold]
|
||||
type: Threshold
|
||||
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-15
|
||||
ezfio_name: threshold_ao
|
||||
|
||||
[ao_cholesky_threshold]
|
||||
type: Threshold
|
||||
doc: If | (ii|jj) | < `ao_cholesky_threshold` then (ii|jj) is zero
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-12
|
||||
|
||||
[do_direct_integrals]
|
||||
type: logical
|
||||
doc: Compute integrals on the fly (very slow, only for debugging)
|
||||
|
||||
@ -4,29 +4,7 @@ BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
|
||||
! Number of Cholesky vectors in AO basis
|
||||
END_DOC
|
||||
|
||||
integer :: i,j,k,l
|
||||
double precision :: xnorm0, x, integral
|
||||
double precision, external :: ao_two_e_integral
|
||||
|
||||
cholesky_ao_num_guess = 0
|
||||
xnorm0 = 0.d0
|
||||
x = 0.d0
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
integral = ao_two_e_integral(i,i,j,j)
|
||||
if (integral > ao_integrals_threshold) then
|
||||
cholesky_ao_num_guess += 1
|
||||
else
|
||||
x += integral
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
print *, 'Cholesky decomposition of AO integrals'
|
||||
print *, '--------------------------------------'
|
||||
print *, ''
|
||||
print *, 'Estimated Error: ', x
|
||||
print *, 'Guess size: ', cholesky_ao_num_guess, '(', 100.d0*dble(cholesky_ao_num_guess)/dble(ao_num*ao_num), ' %)'
|
||||
|
||||
cholesky_ao_num_guess = ao_num*ao_num
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, cholesky_ao_num ]
|
||||
@ -39,7 +17,7 @@ END_PROVIDER
|
||||
END_DOC
|
||||
|
||||
type(c_ptr) :: ptr
|
||||
integer :: fd, i,j,k,l, rank
|
||||
integer :: fd, i,j,k,l,m,rank
|
||||
double precision, pointer :: ao_integrals(:,:,:,:)
|
||||
double precision, external :: ao_two_e_integral
|
||||
|
||||
@ -49,28 +27,83 @@ END_PROVIDER
|
||||
8, fd, .False., ptr)
|
||||
call c_f_pointer(ptr, ao_integrals, (/ao_num, ao_num, ao_num, ao_num/))
|
||||
|
||||
double precision :: integral
|
||||
print*, 'Providing the AO integrals (Cholesky)'
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
|
||||
ao_integrals = 0.d0
|
||||
|
||||
double precision :: integral, cpu_1, cpu_2, wall_1, wall_2
|
||||
logical, external :: ao_two_e_integral_zero
|
||||
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l, integral) SCHEDULE(dynamic)
|
||||
do l=1,ao_num
|
||||
do j=1,l
|
||||
do k=1,ao_num
|
||||
do i=1,k
|
||||
if (ao_two_e_integral_zero(i,j,k,l)) cycle
|
||||
integral = ao_two_e_integral(i,k,j,l)
|
||||
ao_integrals(i,k,j,l) = integral
|
||||
ao_integrals(k,i,j,l) = integral
|
||||
ao_integrals(i,k,l,j) = integral
|
||||
ao_integrals(k,i,l,j) = integral
|
||||
enddo
|
||||
double precision, external :: get_ao_two_e_integral
|
||||
|
||||
if (read_ao_two_e_integrals) then
|
||||
PROVIDE ao_two_e_integrals_in_map
|
||||
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l, integral, wall_2)
|
||||
do m=0,9
|
||||
do l=1+m,ao_num,10
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do j=1,ao_num
|
||||
do k=1,ao_num
|
||||
do i=1,ao_num
|
||||
if (ao_two_e_integral_zero(i,j,k,l)) cycle
|
||||
integral = get_ao_two_e_integral(i,j,k,l, ao_integrals_map)
|
||||
ao_integrals(i,k,j,l) = integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
enddo
|
||||
!$OMP MASTER
|
||||
call wall_time(wall_2)
|
||||
print '(I10,'' % in'', 4X, F10.2, '' s.'')', (m+1) * 10, wall_2-wall_1
|
||||
!$OMP END MASTER
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
else
|
||||
|
||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l, integral, wall_2)
|
||||
do m=0,9
|
||||
do l=1+m,ao_num,10
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do j=1,l
|
||||
do k=1,ao_num
|
||||
do i=1,min(k,j)
|
||||
if (ao_two_e_integral_zero(i,j,k,l)) cycle
|
||||
integral = ao_two_e_integral(i,k,j,l)
|
||||
ao_integrals(i,k,j,l) = integral
|
||||
ao_integrals(k,i,j,l) = integral
|
||||
ao_integrals(i,k,l,j) = integral
|
||||
ao_integrals(k,i,l,j) = integral
|
||||
ao_integrals(j,l,i,k) = integral
|
||||
ao_integrals(j,l,k,i) = integral
|
||||
ao_integrals(l,j,i,k) = integral
|
||||
ao_integrals(l,j,k,i) = integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
enddo
|
||||
!$OMP MASTER
|
||||
call wall_time(wall_2)
|
||||
print '(I10,'' % in'', 4X, F10.2, '' s.'')', (m+1) * 10, wall_2-wall_1
|
||||
!$OMP END MASTER
|
||||
enddo
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall_2)
|
||||
call cpu_time(cpu_2)
|
||||
print*, 'AO integrals provided:'
|
||||
print*, ' cpu time :',cpu_2 - cpu_1, 's'
|
||||
print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
|
||||
|
||||
endif
|
||||
|
||||
! Call Lapack
|
||||
cholesky_ao_num = cholesky_ao_num_guess
|
||||
call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_integrals_threshold, ao_num*ao_num, cholesky_ao)
|
||||
call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_cholesky_threshold, ao_num*ao_num, cholesky_ao)
|
||||
print *, 'Rank: ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
|
||||
|
||||
! Remove mmap
|
||||
|
||||
@ -29,14 +29,14 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
|
||||
complex*16 :: integral5, integral6, integral7, integral8
|
||||
complex*16 :: integral_tot
|
||||
|
||||
double precision :: ao_two_e_integral_cosgtos_schwartz_accel
|
||||
double precision :: ao_2e_cosgtos_schwartz_accel
|
||||
complex*16 :: ERI_cosgtos
|
||||
complex*16 :: general_primitive_integral_cosgtos
|
||||
|
||||
if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then
|
||||
|
||||
!print *, ' with shwartz acc '
|
||||
ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos_schwartz_accel(i, j, k, l)
|
||||
ao_two_e_integral_cosgtos = ao_2e_cosgtos_schwartz_accel(i, j, k, l)
|
||||
|
||||
else
|
||||
!print *, ' without shwartz acc '
|
||||
@ -294,7 +294,7 @@ end function ao_two_e_integral_cosgtos
|
||||
|
||||
! ---
|
||||
|
||||
double precision function ao_two_e_integral_cosgtos_schwartz_accel(i, j, k, l)
|
||||
double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
|
||||
|
||||
BEGIN_DOC
|
||||
! integral of the AO basis <ik|jl> or (ij|kl)
|
||||
@ -329,7 +329,7 @@ double precision function ao_two_e_integral_cosgtos_schwartz_accel(i, j, k, l)
|
||||
complex*16 :: ERI_cosgtos
|
||||
complex*16 :: general_primitive_integral_cosgtos
|
||||
|
||||
ao_two_e_integral_cosgtos_schwartz_accel = 0.d0
|
||||
ao_2e_cosgtos_schwartz_accel = 0.d0
|
||||
|
||||
dim1 = n_pt_max_integrals
|
||||
|
||||
@ -519,8 +519,7 @@ double precision function ao_two_e_integral_cosgtos_schwartz_accel(i, j, k, l)
|
||||
|
||||
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
|
||||
|
||||
ao_two_e_integral_cosgtos_schwartz_accel = ao_two_e_integral_cosgtos_schwartz_accel &
|
||||
+ coef4 * 2.d0 * real(integral_tot)
|
||||
ao_2e_cosgtos_schwartz_accel = ao_2e_cosgtos_schwartz_accel + coef4 * 2.d0 * real(integral_tot)
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
@ -698,8 +697,7 @@ double precision function ao_two_e_integral_cosgtos_schwartz_accel(i, j, k, l)
|
||||
|
||||
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
|
||||
|
||||
ao_two_e_integral_cosgtos_schwartz_accel = ao_two_e_integral_cosgtos_schwartz_accel &
|
||||
+ coef4 * 2.d0 * real(integral_tot)
|
||||
ao_2e_cosgtos_schwartz_accel = ao_2e_cosgtos_schwartz_accel + coef4 * 2.d0 * real(integral_tot)
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
@ -709,11 +707,11 @@ double precision function ao_two_e_integral_cosgtos_schwartz_accel(i, j, k, l)
|
||||
|
||||
deallocate(schwartz_kl)
|
||||
|
||||
end function ao_two_e_integral_cosgtos_schwartz_accel
|
||||
end function ao_2e_cosgtos_schwartz_accel
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_two_e_integral_cosgtos_schwartz, (ao_num,ao_num) ]
|
||||
BEGIN_PROVIDER [ double precision, ao_2e_cosgtos_schwartz, (ao_num,ao_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
! Needed to compute Schwartz inequalities
|
||||
@ -723,16 +721,16 @@ BEGIN_PROVIDER [ double precision, ao_two_e_integral_cosgtos_schwartz, (ao_num,a
|
||||
integer :: i, k
|
||||
double precision :: ao_two_e_integral_cosgtos
|
||||
|
||||
ao_two_e_integral_cosgtos_schwartz(1,1) = ao_two_e_integral_cosgtos(1, 1, 1, 1)
|
||||
ao_2e_cosgtos_schwartz(1,1) = ao_two_e_integral_cosgtos(1, 1, 1, 1)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,k) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP SHARED(ao_num, ao_two_e_integral_cosgtos_schwartz) &
|
||||
!$OMP PARALLEL DO PRIVATE(i,k) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP SHARED(ao_num, ao_2e_cosgtos_schwartz) &
|
||||
!$OMP SCHEDULE(dynamic)
|
||||
do i = 1, ao_num
|
||||
do k = 1, i
|
||||
ao_two_e_integral_cosgtos_schwartz(i,k) = dsqrt(ao_two_e_integral_cosgtos(i, i, k, k))
|
||||
ao_two_e_integral_cosgtos_schwartz(k,i) = ao_two_e_integral_cosgtos_schwartz(i,k)
|
||||
ao_2e_cosgtos_schwartz(i,k) = dsqrt(ao_two_e_integral_cosgtos(i, i, k, k))
|
||||
ao_2e_cosgtos_schwartz(k,i) = ao_2e_cosgtos_schwartz(i,k)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
@ -590,8 +590,20 @@ double precision function general_primitive_integral(dim, &
|
||||
d_poly(i)=0.d0
|
||||
enddo
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
|
||||
! call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
|
||||
integer :: ib, ic
|
||||
if (ior(n_Ix,n_Iy) >= 0) then
|
||||
do ib=0,n_Ix
|
||||
do ic = 0,n_Iy
|
||||
d_poly(ib+ic) = d_poly(ib+ic) + Iy_pol(ic) * Ix_pol(ib)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do n_pt_tmp = n_Ix+n_Iy, 0, -1
|
||||
if (d_poly(n_pt_tmp) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
|
||||
if (n_pt_tmp == -1) then
|
||||
return
|
||||
endif
|
||||
@ -600,8 +612,21 @@ double precision function general_primitive_integral(dim, &
|
||||
d1(i)=0.d0
|
||||
enddo
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
|
||||
! call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
|
||||
if (ior(n_pt_tmp,n_Iz) >= 0) then
|
||||
! Bottleneck here
|
||||
do ib=0,n_pt_tmp
|
||||
do ic = 0,n_Iz
|
||||
d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do n_pt_out = n_pt_tmp+n_Iz, 0, -1
|
||||
if (d1(n_pt_out) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
|
||||
|
||||
double precision :: rint_sum
|
||||
accu = accu + rint_sum(n_pt_out,const,d1)
|
||||
|
||||
@ -948,8 +973,9 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
|
||||
X(ix) *= dble(a-1)
|
||||
enddo
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(X,nx,B_10,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_10,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,B_10,d,nd)
|
||||
|
||||
nx = nd
|
||||
!DIR$ LOOP COUNT(8)
|
||||
@ -970,8 +996,9 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
|
||||
X(ix) *= c
|
||||
enddo
|
||||
endif
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,B_00,d,nd)
|
||||
endif
|
||||
|
||||
ny=0
|
||||
@ -988,9 +1015,9 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
|
||||
call I_x1_pol_mult_recurs(a-1,c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in)
|
||||
endif
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,C_00,d,nd)
|
||||
end
|
||||
|
||||
recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
@ -1028,8 +1055,9 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
enddo
|
||||
endif
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,B_00,d,nd)
|
||||
|
||||
ny=0
|
||||
|
||||
@ -1039,8 +1067,9 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
enddo
|
||||
call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,C_00,d,nd)
|
||||
|
||||
end
|
||||
|
||||
@ -1067,8 +1096,9 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
nx = 0
|
||||
call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(X,nx,B_10,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_10,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,B_10,d,nd)
|
||||
|
||||
nx = nd
|
||||
!DIR$ LOOP COUNT(8)
|
||||
@ -1086,8 +1116,9 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
enddo
|
||||
endif
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,B_00,d,nd)
|
||||
|
||||
ny=0
|
||||
!DIR$ LOOP COUNT(8)
|
||||
@ -1097,9 +1128,9 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
!DIR$ FORCEINLINE
|
||||
call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,C_00,d,nd)
|
||||
end
|
||||
|
||||
recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
@ -1146,8 +1177,10 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
Y(1) = D_00(1)
|
||||
Y(2) = D_00(2)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(Y,ny,D_00,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,D_00,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,D_00,d,nd)
|
||||
|
||||
return
|
||||
|
||||
case default
|
||||
@ -1164,8 +1197,9 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
X(ix) *= dble(c-1)
|
||||
enddo
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(X,nx,B_01,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_01,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,B_01,d,nd)
|
||||
|
||||
ny = 0
|
||||
!DIR$ LOOP COUNT(6)
|
||||
@ -1174,8 +1208,9 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
enddo
|
||||
call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,Y,ny,dim)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_poly(Y,ny,D_00,2,d,nd)
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,D_00,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,D_00,d,nd)
|
||||
|
||||
end select
|
||||
end
|
||||
@ -1233,3 +1268,34 @@ subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine multiply_poly_local(b,nb,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nb, nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:nb), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:nb+nc)
|
||||
|
||||
integer :: ndtmp
|
||||
integer :: ib, ic, id, k
|
||||
if(ior(nc,nb) < 0) return !False if nc>=0 and nb>=0
|
||||
|
||||
do ib=0,nb
|
||||
do ic = 0,nc
|
||||
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do nd = nb+nc,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||