diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 9c0441b8..58df65d9 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -13,11 +13,9 @@ Options: import sys import os -from functools import reduce -from ezfio import ezfio +import tempfile from docopt import docopt - try: QP_ROOT = os.environ["QP_ROOT"] QP_EZFIO = os.environ["QP_EZFIO"] @@ -30,401 +28,41 @@ else: QP_ROOT + "/install", QP_ROOT + "/scripts"] + sys.path -from resultsFile import * +from ezfio import ezfio +from qp_import_trexio import write_ezfio + try: - from resultsFile import * + from trexio_tools import trexio_run except: - print("Error: resultsFile Python library not installed") + print("Error: trexio-tools python package is not installed.") + print("Use pip install trexio-tools to install it.") sys.exit(1) - - -def write_ezfio(res, filename): - - res.clean_uncontractions() - ezfio.set_file(filename) - ezfio.set_ezfio_files_ezfio_convention(20250211) - - # _ - # |_ | _ _ _|_ ._ _ ._ _ - # |_ | (/_ (_ |_ | (_) | | _> - # - print("Electrons\t...\t", end=' ') - ezfio.set_electrons_elec_alpha_num(res.num_alpha) - ezfio.set_electrons_elec_beta_num(res.num_beta) - print("OK") - - # - # |\ | _ | _ o - # | \| |_| (_ | (/_ | - # - - print("Nuclei\t\t...\t", end=' ') - # ~#~#~#~ # - # I n i t # - # ~#~#~#~ # - - charge = [] - coord_x = [] - coord_y = [] - coord_z = [] - - # ~#~#~#~#~#~#~ # - # P a r s i n g # - # ~#~#~#~#~#~#~ # - - for a in res.geometry: - charge.append(a.charge) - if res.units == 'BOHR': - coord_x.append(a.coord[0]) - coord_y.append(a.coord[1]) - coord_z.append(a.coord[2]) - else: - coord_x.append(a.coord[0] / a0) - coord_y.append(a.coord[1] / a0) - coord_z.append(a.coord[2] / a0) - - - # ~#~#~#~#~ # - # W r i t e # - # ~#~#~#~#~ # - - ezfio.set_nuclei_nucl_num(len(res.geometry)) - ezfio.set_nuclei_nucl_charge(charge) - - # Transformt H1 into H - import re - p = re.compile(r'(\d*)$') - label = [p.sub("", x.name).capitalize() for x in res.geometry] - ezfio.set_nuclei_nucl_label(label) - - ezfio.set_nuclei_nucl_coord(coord_x + coord_y + coord_z) - print("OK") - - # _ - # /\ _ _ |_) _. _ o _ - # /--\ (_) _> |_) (_| _> | _> - # - - print("AOS\t\t...\t", end=' ') - # ~#~#~#~ # - # I n i t # - # ~#~#~#~ # - - at = [] - num_prim = [] - power_x = [] - power_y = [] - power_z = [] - coefficient = [] - exponent = [] - - res.convert_to_cartesian() - - # ~#~#~#~#~#~#~ # - # P a r s i n g # - # ~#~#~#~#~#~#~ # - - for b in res.basis: - c = b.center - for i, atom in enumerate(res.geometry): - if atom.coord == c: - at.append(i + 1) - num_prim.append(len(b.prim)) - s = b.sym - power_x.append(str.count(s, "x")) - power_y.append(str.count(s, "y")) - power_z.append(str.count(s, "z")) - coefficient.append(b.coef) - exponent.append([p.expo for p in b.prim]) - - # ~#~#~#~#~ # - # W r i t e # - # ~#~#~#~#~ # - - ezfio.set_ao_basis_ao_num(len(res.basis)) - ezfio.set_ao_basis_ao_nucl(at) - ezfio.set_ao_basis_ao_prim_num(num_prim) - ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) - - # ~#~#~#~#~#~#~ # - # P a r s i n g # - # ~#~#~#~#~#~#~ # - - prim_num_max = ezfio.get_ao_basis_ao_prim_num_max() - - for i in range(len(res.basis)): - coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)] - exponent[i] += [0. for j in range(len(exponent[i]), prim_num_max)] - - coefficient = reduce(lambda x, y: x + y, coefficient, []) - exponent = reduce(lambda x, y: x + y, exponent, []) - - coef = [] - expo = [] - for i in range(prim_num_max): - for j in range(i, len(coefficient), prim_num_max): - coef.append(coefficient[j]) - expo.append(exponent[j]) - - # ~#~#~#~#~ # - # W r i t e # - # ~#~#~#~#~ # - - ezfio.set_ao_basis_ao_coef(coef) - ezfio.set_basis_ao_normalized(True) - ezfio.set_ao_basis_ao_expo(expo) - ezfio.set_ao_basis_ao_basis("Read by resultsFile") - - print("OK") - - # _ - # |_) _. _ o _ - # |_) (_| _> | _> - # - - print("Basis\t\t...\t", end=' ') - # ~#~#~#~ # - # I n i t # - # ~#~#~#~ # - - coefficient = [] - exponent = [] - - # ~#~#~#~#~#~#~ # - # P a r s i n g # - # ~#~#~#~#~#~#~ # - - inucl = {} - for i, a in enumerate(res.geometry): - inucl[a.coord] = i - - nbasis = 0 - nucl_index = [] - curr_center = -1 - nucl_shell_num = [] - ang_mom = [] - nshell = 0 - nshell_tot = 0 - shell_index = [] - shell_prim_num = [] - for b in res.basis: - s = b.sym - if str.count(s, "y") + str.count(s, "x") == 0: - c = inucl[b.center] - nshell += 1 - nshell_tot += 1 - if c != curr_center: - curr_center = c - nucl_shell_num.append(nshell) - nshell = 0 - nbasis += 1 - nucl_index.append(c+1) - coefficient += b.coef[:len(b.prim)] - exponent += [p.expo for p in b.prim] - ang_mom.append(str.count(s, "z")) - shell_prim_num.append(len(b.prim)) - shell_index += [nshell_tot] * len(b.prim) - - shell_num = len(ang_mom) - assert(shell_index[0] == 1) - assert(shell_index[-1] == shell_num) - - # ~#~#~#~#~ # - # W r i t e # - # ~#~#~#~#~ # - - ezfio.set_basis_basis("Read from ResultsFile") - ezfio.set_basis_shell_num(shell_num) - ezfio.set_basis_basis_nucleus_index(nucl_index) - ezfio.set_basis_prim_num(len(coefficient)) - - ezfio.set_basis_nucleus_shell_num(nucl_shell_num) - ezfio.set_basis_prim_coef(coefficient) - ezfio.set_basis_prim_expo(exponent) - ezfio.set_basis_shell_ang_mom(ang_mom) - ezfio.set_basis_shell_prim_num(shell_prim_num) - ezfio.set_basis_shell_index(shell_index) - - print("OK") - - # _ - # |\/| _ _ |_) _. _ o _ - # | | (_) _> |_) (_| _> | _> - # - - print("MOS\t\t...\t", end=' ') - # ~#~#~#~ # - # I n i t # - # ~#~#~#~ # - - MoTag = res.determinants_mo_type - ezfio.set_mo_basis_mo_label('Orthonormalized') - ezfio.set_determinants_mo_label('Orthonormalized') - MO_type = MoTag - allMOs = res.mo_sets[MO_type] - - # ~#~#~#~#~#~#~ # - # P a r s i n g # - # ~#~#~#~#~#~#~ # - - try: - closed = [(allMOs[i].eigenvalue, i) for i in res.closed_mos] - active = [(allMOs[i].eigenvalue, i) for i in res.active_mos] - virtual = [(allMOs[i].eigenvalue, i) for i in res.virtual_mos] - except: - closed = [] - virtual = [] - active = [(allMOs[i].eigenvalue, i) for i in range(len(allMOs))] - - closed = [x[1] for x in closed] - active = [x[1] for x in active] - virtual = [x[1] for x in virtual] - MOindices = closed + active + virtual - - MOs = [] - for i in MOindices: - MOs.append(allMOs[i]) - - mo_num = len(MOs) - while len(MOindices) < mo_num: - MOindices.append(len(MOindices)) - - MOmap = list(MOindices) - for i in range(len(MOindices)): - MOmap[i] = MOindices.index(i) - - energies = [] - for i in range(mo_num): - energies.append(MOs[i].eigenvalue) - - OccNum = [] - if res.occ_num is not None: - for i in MOindices: - OccNum.append(res.occ_num[MO_type][i]) - else: - for i in range(res.num_beta): - OccNum.append(2.) - for i in range(res.num_beta,res.num_alpha): - OccNum.append(1.) - - while len(OccNum) < mo_num: - OccNum.append(0.) - - MoMatrix = [] - sym0 = [i.sym for i in res.mo_sets[MO_type]] - sym = [i.sym for i in res.mo_sets[MO_type]] - for i in range(len(sym)): - sym[MOmap[i]] = sym0[i] - - irrep = {} - for i in sym: - irrep[i] = 0 - - for i, j in enumerate(irrep.keys()): - irrep[j] = i+1 - - sym = [ irrep[k] for k in sym ] - - MoMatrix = [] - for i in range(len(MOs)): - m = MOs[i] - for coef in m.vector: - MoMatrix.append(coef) - - while len(MoMatrix) < len(MOs[0].vector)**2: - MoMatrix.append(0.) - - # ~#~#~#~#~ # - # W r i t e # - # ~#~#~#~#~ # - - ezfio.set_mo_basis_mo_num(mo_num) - ezfio.set_mo_basis_mo_coef(MoMatrix) - ezfio.set_mo_basis_mo_occ(OccNum) - ezfio.set_mo_basis_mo_symmetry(sym) - - print("OK") - - - print("Pseudos\t\t...\t", end=' ') - try: - lmax = 0 - nucl_charge_remove = [] - klocmax = 0 - kmax = 0 - nucl_num = len(res.geometry) - for ecp in res.pseudo: - lmax_local = ecp['lmax'] - lmax = max(lmax_local, lmax) - nucl_charge_remove.append(ecp['zcore']) - klocmax = max(klocmax, len(ecp[str(lmax_local)])) - for l in range(lmax_local): - kmax = max(kmax, len(ecp[str(l)])) - lmax = lmax-1 - ezfio.set_pseudo_pseudo_lmax(lmax) - ezfio.set_pseudo_nucl_charge_remove(nucl_charge_remove) - ezfio.set_pseudo_pseudo_klocmax(klocmax) - ezfio.set_pseudo_pseudo_kmax(kmax) - pseudo_n_k = [[0 for _ in range(nucl_num)] for _ in range(klocmax)] - pseudo_v_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)] - pseudo_dz_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)] - pseudo_n_kl = [[[0 for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)] - pseudo_v_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)] - pseudo_dz_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)] - for ecp in res.pseudo: - lmax_local = ecp['lmax'] - klocmax = len(ecp[str(lmax_local)]) - atom = ecp['atom']-1 - for kloc in range(klocmax): - try: - v, n, dz = ecp[str(lmax_local)][kloc] - pseudo_n_k[kloc][atom] = n-2 - pseudo_v_k[kloc][atom] = v - pseudo_dz_k[kloc][atom] = dz - except: - pass - for l in range(lmax_local): - for k in range(kmax): - try: - v, n, dz = ecp[str(l)][k] - pseudo_n_kl[l][k][atom] = n-2 - pseudo_v_kl[l][k][atom] = v - pseudo_dz_kl[l][k][atom] = dz - except: - pass - ezfio.set_pseudo_pseudo_n_k(pseudo_n_k) - ezfio.set_pseudo_pseudo_v_k(pseudo_v_k) - ezfio.set_pseudo_pseudo_dz_k(pseudo_dz_k) - 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) - - n_alpha = res.num_alpha - n_beta = res.num_beta - for i in range(nucl_num): - charge[i] -= nucl_charge_remove[i] - ezfio.set_nuclei_nucl_charge(charge) - ezfio.set_electrons_elec_alpha_num(n_alpha) - ezfio.set_electrons_elec_beta_num(n_beta) - - except: - ezfio.set_pseudo_do_pseudo(False) - else: - ezfio.set_pseudo_do_pseudo(True) - - print("OK") - - - - def get_full_path(file_path): file_path = os.path.expanduser(file_path) file_path = os.path.expandvars(file_path) -# file_path = os.path.abspath(file_path) return file_path +class NoneDict(dict): + def __missing__(self, key): + return None + +def main(FILE,EZFIO_FILE): + with tempfile.NamedTemporaryFile(mode='w+b', delete=True, errors=None) as f: + trexio_file = f.name + + args = NoneDict() + args["convert-from"] = True + args["--input"] = FILE + args["--back_end"] = "hdf5" + args["--type"] = "gaussian" + args["TREXIO_FILE"] = trexio_file + + trexio_run.main(filename=trexio_file, args=args) + write_ezfio(trexio_file, EZFIO_FILE) + trexio_run.remove_trexio_file(trexio_file, overwrite=True) + if __name__ == '__main__': ARGUMENTS = docopt(__doc__) @@ -436,20 +74,5 @@ if __name__ == '__main__': else: EZFIO_FILE = "{0}.ezfio".format(FILE) - try: - RES_FILE = getFile(FILE) - except: - raise - else: - print(FILE, 'recognized as', str(RES_FILE).split('.')[-1].split()[0]) - - write_ezfio(RES_FILE, EZFIO_FILE) - sys.stdout.flush() - if os.system("qp_run save_ortho_mos "+EZFIO_FILE) != 0: - print("""Warning: You need to run - - qp run save_ortho_mos - -to be sure your MOs will be orthogonal, which is not the case when -the MOs are read from output files (not enough precision in output).""") + main (FILE,EZFIO_FILE) diff --git a/plugins/local/basis_correction/print_routine.irp.f b/plugins/local/basis_correction/print_routine.irp.f index cc9744f2..e4f42510 100644 --- a/plugins/local/basis_correction/print_routine.irp.f +++ b/plugins/local/basis_correction/print_routine.irp.f @@ -43,7 +43,8 @@ subroutine print_basis_correction else if(mu_of_r_potential.EQ."cas_full".or. & mu_of_r_potential.EQ."cas_truncated".or. & - mu_of_r_potential.EQ."pure_act") then + mu_of_r_potential.EQ."pure_act".or. & + mu_of_r_potential.EQ."proj_cas") then print*, '' print*,'Using a CAS-like two-body density to define mu(r)' print*,'This assumes that the CAS is a qualitative representation of the wave function ' diff --git a/scripts/qp_cipsi_rsh b/scripts/qp_cipsi_rsh deleted file mode 120000 index c3d4376b..00000000 --- a/scripts/qp_cipsi_rsh +++ /dev/null @@ -1 +0,0 @@ -/home/scemama/qp2/plugins/qp_plugins_lct/stable/rsdft_cipsi/qp_cipsi_rsh \ No newline at end of file diff --git a/scripts/qp_exc_energy.py b/scripts/qp_exc_energy.py index e08866e3..e3060349 100755 --- a/scripts/qp_exc_energy.py +++ b/scripts/qp_exc_energy.py @@ -1,6 +1,9 @@ #!/usr/bin/env python3 # Computes the error on the excitation energy of a CIPSI run. +# see "QUESTDB: a database of highly-accurate excitation energies for the electronic structure community" +# doi:10.1002/wcms.1517 or arXiv:2011.14675 + def student(p,df): import scipy diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index 156c4344..72f24950 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -80,6 +80,8 @@ def generate_xyz(l): def write_ezfio(trexio_filename, filename): warnings = [] + while trexio_filename[-1] == '/': + trexio_filename = trexio_filename[:-1] trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_AUTO) ezfio.set_file(filename) @@ -258,7 +260,6 @@ def write_ezfio(trexio_filename, filename): else: raise TypeError - print(basis_type) except: basis_type = "None" print("None") @@ -275,7 +276,7 @@ def write_ezfio(trexio_filename, filename): if basis_type.lower() == "gaussian" and not cartesian: try: import trexio_tools - tmp = "cartesian_"+trexio_filename + tmp = trexio_filename+"_cartesian" retcode = subprocess.call(["trexio", "convert-to", "-t", "cartesian", "-o", tmp, trexio_filename]) trexio_file_cart = trexio.File(tmp,mode='r',back_end=trexio.TREXIO_AUTO) cartesian = trexio.read_ao_cartesian(trexio_file_cart) @@ -357,7 +358,7 @@ def write_ezfio(trexio_filename, filename): else: if basis_type.lower() == "gaussian" and not cartesian: - warnings.append(f"Spherical AOs not handled by QP. Convert the TREXIO file using trexio_tools:\n trexio convert-to -t cartesian -o cartesian_{trexio_filename} {trexio_filename}") + warnings.append(f"Spherical AOs not handled by QP. Convert the TREXIO file using trexio_tools:\n trexio convert-to -t cartesian -o {trexio_filename}_cartesian {trexio_filename}") warnings.append("Integrals should be imported using:\n qp run import_trexio_integrals") print("None") diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 02eedf53..ef871538 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -28,11 +28,15 @@ BEGIN_PROVIDER [ integer, ao_sphe_num ] ! Number of spherical AOs END_DOC integer :: n, i - ao_sphe_num=0 - do i=1,shell_num - n = shell_ang_mom(i) - ao_sphe_num += 2*n+1 - enddo + if (ao_cartesian) then + ao_sphe_num = ao_num + else + ao_sphe_num=0 + do i=1,shell_num + n = shell_ang_mom(i) + ao_sphe_num += 2*n+1 + enddo + endif END_PROVIDER BEGIN_PROVIDER [ integer, ao_sphe_shell, (ao_sphe_num) ] diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f index 523e49f7..e015c89e 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f @@ -14,45 +14,55 @@ prev = 0 ao_cart_to_sphe_coef(:,:) = 0.d0 ao_cart_to_sphe_normalization(:) = 1.d0 - ! Assume order provided by ao_power_index - i = 1 - ao_sphe_count = 0 - do while (i <= ao_num) - select case ( ao_l(i) ) - case (0) - ao_sphe_count += 1 - ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0 - ao_cart_to_sphe_normalization(i) = 1.d0 - i += 1 - BEGIN_TEMPLATE - case ($SHELL) - if (ao_power(i,1) == $SHELL) then - do k=1,size(cart_to_sphe_$SHELL,2) - do j=1,size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k) + + if (ao_cartesian) then + ! Identity matrix + do i=1,ao_sphe_num + ao_cart_to_sphe_coef(i,i) = 1.d0 + enddo + + else + ! Assume order provided by ao_power_index + i = 1 + ao_sphe_count = 0 + do while (i <= ao_num) + select case ( ao_l(i) ) + case (0) + ao_sphe_count += 1 + ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0 + ao_cart_to_sphe_normalization(i) = 1.d0 + i += 1 + BEGIN_TEMPLATE + case ($SHELL) + if (ao_power(i,1) == $SHELL) then + do k=1,size(cart_to_sphe_$SHELL,2) + do j=1,size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k) + enddo enddo - enddo - do j=1,size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j) - enddo - i += size(cart_to_sphe_$SHELL,1) - ao_sphe_count += size(cart_to_sphe_$SHELL,2) - endif - SUBST [ SHELL ] - 1;; - 2;; - 3;; - 4;; - 5;; - 6;; - 7;; - 8;; - 9;; - END_TEMPLATE - case default - stop 'Error in ao_cart_to_sphe : angular momentum too high' - end select - enddo + do j=1,size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j) + enddo + i += size(cart_to_sphe_$SHELL,1) + ao_sphe_count += size(cart_to_sphe_$SHELL,2) + endif + SUBST [ SHELL ] + 1;; + 2;; + 3;; + 4;; + 5;; + 6;; + 7;; + 8;; + 9;; + END_TEMPLATE + case default + stop 'Error in ao_cart_to_sphe : angular momentum too high' + end select + enddo + + endif if (ao_sphe_count /= ao_sphe_num) then call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num") diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f index 851f26d8..f752614c 100644 --- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -116,13 +116,14 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] alpha = ao_expo_ordered_transp(l,j) do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - c = 0.d0 - if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& < thresh) then cycle endif + + beta = ao_expo_ordered_transp(m,i) + c = 0.d0 + do k = 1, nucl_num Z = nucl_charge(k) @@ -154,6 +155,12 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] !$OMP END DO !$OMP END PARALLEL + do i=1,ao_num + do j=1,i + ao_pseudo_integrals_local(j,i) = 0.5d0*(ao_pseudo_integrals_local(i,j) + ao_pseudo_integrals_local(i,j)) + ao_pseudo_integrals_local(i,j) = ao_pseudo_integrals_local(i,j) + enddo + enddo END_PROVIDER @@ -216,14 +223,14 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] alpha = ao_expo_ordered_transp(l,j) do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - c = 0.d0 - if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& < thresh) then cycle endif + beta = ao_expo_ordered_transp(m,i) + c = 0.d0 + do k = 1, nucl_num Z = nucl_charge(k) @@ -256,6 +263,12 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] !$OMP END PARALLEL + do i=1,ao_num + do j=1,i + ao_pseudo_integrals_non_local(j,i) = 0.5d0*(ao_pseudo_integrals_non_local(i,j) + ao_pseudo_integrals_non_local(i,j)) + ao_pseudo_integrals_non_local(i,j) = ao_pseudo_integrals_non_local(i,j) + enddo + enddo END_PROVIDER BEGIN_PROVIDER [ double precision, pseudo_v_k_transp, (pseudo_klocmax,nucl_num) ] diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index a04eae69..1cb7617e 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -649,20 +649,11 @@ double precision function general_primitive_integral(dim, & ! 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 - if (ic > ib) then - do ib=0,n_pt_tmp - d1(ib:) = d1(ib:) + Iz_pol(:) * d_poly(ib) - enddo - else + do ib=0,n_pt_tmp do ic = 0,n_Iz - d1(ic:) = d1(ic:) + Iz_pol(ic) * d_poly(:) + d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib) enddo - endif -! 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 + enddo do n_pt_out = n_pt_tmp+n_Iz, 0, -1 if (d1(n_pt_out) /= 0.d0) exit diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index 2f497bd7..b087ec24 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -16,11 +16,10 @@ integer*8 function spin_det_search_key(det,Nint) integer(bit_kind), intent(in) :: det(Nint) integer(bit_kind), parameter :: unsigned_shift = 1_bit_kind-huge(1_bit_kind) ! 100...00 integer :: i - spin_det_search_key = det(1) + spin_det_search_key = det(1)+unsigned_shift do i=2,Nint spin_det_search_key = ieor(spin_det_search_key,det(i)) enddo - spin_det_search_key = spin_det_search_key+unsigned_shift end diff --git a/src/mo_basis/mo_class.irp.f b/src/mo_basis/mo_class.irp.f index 7705e414..8017e119 100644 --- a/src/mo_basis/mo_class.irp.f +++ b/src/mo_basis/mo_class.irp.f @@ -31,3 +31,42 @@ BEGIN_PROVIDER [ character*(32), mo_class , (mo_num) ] IRP_ENDIF END_PROVIDER + + + +BEGIN_PROVIDER [ integer, mo_symmetry , (mo_num) ] + implicit none + BEGIN_DOC +! MOs with the same integer belong to the same irrep. + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + if (size(mo_symmetry) == 0) return + + call ezfio_has_mo_basis_mo_symmetry(has) + if (has) then +! write(6,'(A)') '.. >>>>> [ IO READ: mo_symmetry ] <<<<< ..' + call ezfio_get_mo_basis_mo_symmetry(mo_symmetry) + else + mo_symmetry(:) = 1 + call ezfio_set_mo_basis_mo_symmetry(mo_symmetry) + 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( mo_symmetry, (mo_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_symmetry with MPI' + endif + IRP_ENDIF + +! call write_time(6) + +END_PROVIDER diff --git a/src/mo_one_e_ints/ao_to_mo.irp.f b/src/mo_one_e_ints/ao_to_mo.irp.f index 7ebc4638..72c1328d 100644 --- a/src/mo_one_e_ints/ao_to_mo.irp.f +++ b/src/mo_one_e_ints/ao_to_mo.irp.f @@ -73,3 +73,87 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_from_mo, (ao_num, ao_num)] END_DOC call mo_to_ao(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_num) END_PROVIDER + + + +! --- + + +subroutine mo_to_ao_sphe(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the AO basis + ! + ! $(S.C).A_{mo}.(S.C)^\dagger$ + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + double precision, intent(in) :: A_mo(LDA_mo,mo_num) + double precision, intent(out) :: A_ao(LDA_ao,ao_sphe_num) + double precision, allocatable :: T(:,:) + + allocate ( T(mo_num,ao_sphe_num) ) + + call dgemm('N','T', mo_num, ao_sphe_num, mo_num, & + 1.d0, A_mo,size(A_mo,1), & + S_mo_sphe_coef, size(S_mo_sphe_coef,1), & + 0.d0, T, size(T,1)) + + call dgemm('N','N', ao_sphe_num, ao_sphe_num, mo_num, & + 1.d0, S_mo_sphe_coef, size(S_mo_sphe_coef,1), & + T, size(T,1), & + 0.d0, A_ao, size(A_ao,1)) + + deallocate(T) +end + +subroutine mo_to_ao_sphe_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! $C.A_{mo}.C^\dagger$ + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + double precision, intent(in) :: A_mo(LDA_mo,mo_num) + double precision, intent(out) :: A_ao(LDA_ao,ao_sphe_num) + double precision, allocatable :: T(:,:) + + allocate ( T(mo_num,ao_sphe_num) ) + + call dgemm('N','T', mo_num, ao_sphe_num, mo_num, & + 1.d0, A_mo,size(A_mo,1), & + mo_sphe_coef, size(mo_sphe_coef,1), & + 0.d0, T, size(T,1)) + + call dgemm('N','N', ao_sphe_num, ao_sphe_num, mo_num, & + 1.d0, mo_sphe_coef, size(mo_sphe_coef,1), & + T, size(T,1), & + 0.d0, A_ao, size(A_ao,1)) + + deallocate(T) +end + +BEGIN_PROVIDER [ double precision, S_mo_sphe_coef, (ao_sphe_num, mo_num) ] + implicit none + BEGIN_DOC + ! Product S.C where S is the overlap matrix in the AO basis and C the mo_sphe_coef matrix. + END_DOC + + call dgemm('N','N', ao_sphe_num, mo_num, ao_sphe_num, & + 1.d0, ao_overlap,size(ao_overlap,1), & + mo_sphe_coef, size(mo_coef,1), & + 0.d0, S_mo_coef, size(S_mo_coef,1)) + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals_from_mo, (ao_sphe_num, ao_sphe_num)] + implicit none + BEGIN_DOC +! Integrals of the one e hamiltonian obtained from the integrals on the MO basis +! +! WARNING : this is equal to ao_one_e_integrals only if the AO and MO basis have the same number of functions + END_DOC + call mo_to_ao_sphe(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_sphe_num) +END_PROVIDER + + + diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index d13605f6..f0b81eb8 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -38,7 +38,7 @@ do ipoint = 1, n_points_final_grid mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) enddo - else if(mu_of_r_potential.EQ."proj")then + else if((mu_of_r_potential.EQ."proj").or.(mu_of_r_potential.EQ."proj_cas"))then do ipoint = 1, n_points_final_grid mu_of_r_prov(ipoint,istate) = mu_of_r_projector_mo(ipoint) enddo