diff --git a/bin/qp_set_frozen_large_core b/bin/qp_set_frozen_large_core new file mode 100755 index 00000000..cce5ca70 --- /dev/null +++ b/bin/qp_set_frozen_large_core @@ -0,0 +1,79 @@ +#!/usr/bin/env python2 + + +""" +Automatically finds n, the number of core electrons. Calls qp_set_mo_class +setting all MOs as Active, except the n/2 first ones which are set as Core. +If pseudo-potentials are used, all the MOs are set as Active. + + + +Usage: + qp_set_frozen_core [-q|--query] EZFIO_DIR + +Options: + -q --query Prints in the standard output the number of frozen MOs + +""" + +import os +import sys +import os.path + +try: + import qp_path +except ImportError: + print "source .quantum_package.rc" + raise + +from docopt import docopt +from ezfio import ezfio + + +def main(arguments): + """Main function""" + + filename = arguments["EZFIO_DIR"] + ezfio.set_filename(filename) + + n_frozen = 0 + try: + do_pseudo = ezfio.pseudo_do_pseudo + except: + do_pseudo = False + + if not do_pseudo: + for charge in ezfio.nuclei_nucl_charge: + if charge <= 2: + pass + elif charge <= 10: + n_frozen += 1 + elif charge <= 18: + n_frozen += 5 + elif charge <= 36: + n_frozen += 9 + elif charge <= 54: + n_frozen += 18 + elif charge <= 86: + n_frozen += 27 + elif charge <= 118: + n_frozen += 43 + + mo_num = ezfio.mo_basis_mo_num + + if arguments["--query"]: + print n_frozen + sys.exit(0) + + if n_frozen == 0: + os.system("""qp_set_mo_class -a "[1-%d]" %s""" % + (mo_num, sys.argv[1])) + else: + os.system("""qp_set_mo_class -c "[1-%d]" -a "[%d-%d]" %s""" % + (n_frozen, n_frozen+1, mo_num, sys.argv[1])) + + + +if __name__ == '__main__': + ARGUMENTS = docopt(__doc__) + main(ARGUMENTS) diff --git a/configure b/configure index 02705017..ffda5016 100755 --- a/configure +++ b/configure @@ -387,7 +387,7 @@ if [[ ${ZEROMQ} = $(not_found) ]] ; then fail fi -F77ZMQ=$(find_lib -lzmq -lf77zmq) +F77ZMQ=$(find_lib -lzmq -lf77zmq -lpthread) if [[ ${F77ZMQ} = $(not_found) ]] ; then error "Fortran binding of ZeroMQ (f77zmq) is not installed." fail diff --git a/etc/.gitignore b/etc/.gitignore index 6fbf4389..1f5036f0 100644 --- a/etc/.gitignore +++ b/etc/.gitignore @@ -1 +1,2 @@ 00.qp_root +local.rc diff --git a/etc/local.rc b/etc/local.rc index 83d4b34b..954b315b 100644 --- a/etc/local.rc +++ b/etc/local.rc @@ -16,6 +16,7 @@ # export OMP_NUM_THREADS=16 # Name of the network interface to be chosen +# export QP_NIC=lo # export QP_NIC=ib0 diff --git a/ocaml/MO_label.ml b/ocaml/MO_label.ml index e23b2542..9396668e 100644 --- a/ocaml/MO_label.ml +++ b/ocaml/MO_label.ml @@ -6,6 +6,7 @@ type t = | Natural | Localized | Orthonormalized +| MCSCF | None [@@deriving sexp] @@ -16,6 +17,7 @@ let to_string = function | Orthonormalized -> "Orthonormalized" | Natural -> "Natural" | Localized -> "Localized" + | MCSCF -> "MCSCF" | None -> "None" ;; @@ -26,7 +28,8 @@ let of_string s = | "natural" -> Natural | "localized" -> Localized | "orthonormalized" -> Orthonormalized + | "mcscf" -> MCSCF | "none" -> None | _ -> (print_endline s ; failwith "MO_label should be one of: -Guess | Orthonormalized | Canonical | Natural | Localized | None.") +Guess | Orthonormalized | Canonical | Natural | Localized | MCSCF | None.") ;; diff --git a/ocaml/MO_label.mli b/ocaml/MO_label.mli index 732bf1f2..dd4bb7bd 100644 --- a/ocaml/MO_label.mli +++ b/ocaml/MO_label.mli @@ -4,6 +4,7 @@ type t = | Natural | Localized | Orthonormalized + | MCSCF | None [@@deriving sexp] diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index b3e88fdb..0f63308b 100644 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -5,6 +5,8 @@ import os keywords = """ check_double_excitation copy_buffer +filter_only_connected_to_hf_single +filter_only_connected_to_hf_double declarations decls_main deinit_thread @@ -205,84 +207,84 @@ class H_apply(object): def filter_only_2h(self): self["only_2h_single"] = """ ! ! DIR$ FORCEINLINE - if (is_a_2h(hole).eqv. .False.) cycle + if (.not.is_a_2h(hole)) cycle """ self["only_2h_double"] = """ ! ! DIR$ FORCEINLINE - if ( is_a_2h(key).eqv. .False. )cycle + if (.not.is_a_2h(key))cycle """ def filter_only_1h(self): self["only_1h_single"] = """ ! ! DIR$ FORCEINLINE - if (is_a_1h(hole) .eqv. .False.) cycle + if (.not.is_a_1h(hole)) cycle """ self["only_1h_double"] = """ ! ! DIR$ FORCEINLINE - if (is_a_1h(key) .eqv. .False.) cycle + if (.not.is_a_1h(key) ) cycle """ def filter_only_1p(self): self["only_1p_single"] = """ ! ! DIR$ FORCEINLINE - if ( is_a_1p(hole) .eqv. .False.) cycle + if (.not. is_a_1p(hole) ) cycle """ self["only_1p_double"] = """ ! ! DIR$ FORCEINLINE - if ( is_a_1p(key) .eqv. .False.) cycle + if (.not. is_a_1p(key) ) cycle """ def filter_only_2h1p(self): self["only_2h1p_single"] = """ ! ! DIR$ FORCEINLINE - if ( is_a_2h1p(hole) .eqv. .False.) cycle + if (.not. is_a_2h1p(hole) ) cycle """ self["only_2h1p_double"] = """ ! ! DIR$ FORCEINLINE - if (is_a_2h1p(key) .eqv. .False.) cycle + if (.not.is_a_2h1p(key) ) cycle """ def filter_only_2p(self): self["only_2p_single"] = """ ! ! DIR$ FORCEINLINE - if (is_a_2p(hole).eqv. .False.) cycle + if (.not.is_a_2p(hole)) cycle """ self["only_2p_double"] = """ ! ! DIR$ FORCEINLINE - if (is_a_2p(key).eqv. .False.) cycle + if (.not.is_a_2p(key)) cycle """ def filter_only_1h1p(self): self["filter_only_1h1p_single"] = """ -! ! DIR$ FORCEINLINE - if (is_a_1h1p(hole).eqv..False.) cycle + if (.not.is_a_1h1p(hole)) cycle """ self["filter_only_1h1p_double"] = """ -! ! DIR$ FORCEINLINE - if (is_a_1h1p(key).eqv..False.) cycle + if (.not.is_a_1h1p(key)) cycle """ + + def filter_only_2h2p(self): self["filter_only_2h2p_single"] = """ ! ! DIR$ FORCEINLINE - if (is_a_two_holes_two_particles(hole).eqv..False.) cycle + if (.not.is_a_two_holes_two_particles(hole)) cycle """ self["filter_only_2h2p_double"] = """ ! ! DIR$ FORCEINLINE - if (is_a_two_holes_two_particles(key).eqv..False.) cycle + if (.not.is_a_two_holes_two_particles(key)) cycle """ def filter_only_1h2p(self): self["filter_only_1h2p_single"] = """ ! ! DIR$ FORCEINLINE - if (is_a_1h2p(hole).eqv..False.) cycle + if (.not.is_a_1h2p(hole)) cycle """ self["filter_only_1h2p_double"] = """ ! ! DIR$ FORCEINLINE - if (is_a_1h2p(key).eqv..False.) cycle + if (.not.is_a_1h2p(key)) cycle """ @@ -294,6 +296,16 @@ class H_apply(object): if (is_a_two_holes_two_particles(hole)) cycle """ + def filter_only_connected_to_hf(self): + self["filter_only_connected_to_hf_single"] = """ + call connected_to_hf(hole,yes_no) + if (.not.yes_no) cycle + """ + self["filter_only_connected_to_hf_double"] = """ + call connected_to_hf(key,yes_no) + if (.not.yes_no) cycle + """ + def set_perturbation(self,pert): if self.perturbation is not None: diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 2f508a08..304fec49 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -213,6 +213,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) BEGIN_DOC ! Gets multiple AO bi-electronic integral from the AO map . ! All i are retrieved for j,k,l fixed. + ! physicist convention : END_DOC implicit none integer, intent(in) :: j,k,l, sze diff --git a/src/aux_quantities/EZFIO.cfg b/src/aux_quantities/EZFIO.cfg index 2e1c5b12..6b4bd0f2 100644 --- a/src/aux_quantities/EZFIO.cfg +++ b/src/aux_quantities/EZFIO.cfg @@ -24,3 +24,17 @@ type: double precision size: (mo_basis.mo_num,mo_basis.mo_num,determinants.n_states) +[data_one_e_dm_alpha_ao] +interface: ezfio, provider +doc: Alpha one body density matrix on the |AO| basis computed with the wave function +type: double precision +size: (ao_basis.ao_num,ao_basis.ao_num,determinants.n_states) + + +[data_one_e_dm_beta_ao] +interface: ezfio, provider +doc: Beta one body density matrix on the |AO| basis computed with the wave function +type: double precision +size: (ao_basis.ao_num,ao_basis.ao_num,determinants.n_states) + + diff --git a/src/bitmask/bitmasks.irp.f b/src/bitmask/bitmasks.irp.f index 86b478d6..d425dda6 100644 --- a/src/bitmask/bitmasks.irp.f +++ b/src/bitmask/bitmasks.irp.f @@ -466,35 +466,17 @@ END_PROVIDER BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)] -&BEGIN_PROVIDER [ integer, n_core_inact_act_orb ] implicit none BEGIN_DOC ! Reunion of the core, inactive and active bitmasks END_DOC integer :: i,j - n_core_inact_act_orb = 0 do i = 1, N_int - reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1)) - reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,2,1)) - n_core_inact_act_orb +=popcnt(reunion_of_core_inact_act_bitmask(i,1)) + reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),act_bitmask(i,1)) + reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),act_bitmask(i,2)) enddo END_PROVIDER - BEGIN_PROVIDER [ integer, list_core_inact_act, (n_core_inact_act_orb)] -&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (mo_num)] - implicit none - integer :: occ_inact(N_int*bit_kind_size) - integer :: itest,i - occ_inact = 0 - call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), occ_inact(1), itest, N_int) - list_inact_reverse = 0 - do i = 1, n_core_inact_act_orb - list_core_inact_act(i) = occ_inact(i) - list_core_inact_act_reverse(occ_inact(i)) = i - enddo -END_PROVIDER - - BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] @@ -563,8 +545,8 @@ END_PROVIDER END_DOC integer :: i,j do i = 1, N_int - reunion_of_cas_inact_bitmask(i,1) = ior(cas_bitmask(i,1,1),inact_bitmask(i,1)) - reunion_of_cas_inact_bitmask(i,2) = ior(cas_bitmask(i,2,1),inact_bitmask(i,2)) + reunion_of_cas_inact_bitmask(i,1) = ior(act_bitmask(i,1),inact_bitmask(i,1)) + reunion_of_cas_inact_bitmask(i,2) = ior(act_bitmask(i,2),inact_bitmask(i,2)) enddo END_PROVIDER diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index e384de64..f830da4e 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -194,3 +194,53 @@ END_PROVIDER END_PROVIDER +BEGIN_PROVIDER [integer, n_inact_act_orb ] + implicit none + n_inact_act_orb = (n_inact_orb+n_act_orb) + +END_PROVIDER + +BEGIN_PROVIDER [integer, list_inact_act, (n_inact_act_orb)] + integer :: i,itmp + itmp = 0 + do i = 1, n_inact_orb + itmp += 1 + list_inact_act(itmp) = list_inact(i) + enddo + do i = 1, n_act_orb + itmp += 1 + list_inact_act(itmp) = list_act(i) + enddo +END_PROVIDER + +BEGIN_PROVIDER [integer, n_core_inact_act_orb ] + implicit none + n_core_inact_act_orb = (n_core_orb + n_inact_orb + n_act_orb) + +END_PROVIDER + + BEGIN_PROVIDER [integer, list_core_inact_act, (n_core_inact_act_orb)] +&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (n_core_inact_act_orb)] + integer :: i,itmp + itmp = 0 + do i = 1, n_core_orb + itmp += 1 + list_core_inact_act(itmp) = list_core(i) + enddo + do i = 1, n_inact_orb + itmp += 1 + list_core_inact_act(itmp) = list_inact(i) + enddo + do i = 1, n_act_orb + itmp += 1 + list_core_inact_act(itmp) = list_act(i) + enddo + + integer :: occ_inact(N_int*bit_kind_size) + occ_inact = 0 + call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), occ_inact(1), itest, N_int) + list_inact_reverse = 0 + do i = 1, n_core_inact_act_orb + list_core_inact_act_reverse(occ_inact(i)) = i + enddo +END_PROVIDER diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 9cb94a43..eedcd67f 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -627,6 +627,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) ! endif end do + if(pseudo_sym)then + if(dabs(mat(1, p1, p2)).lt.thresh_sym)then + sum_e_pert = 10.d0 + endif + endif if(sum_e_pert <= buf%mini) then call add_to_selection_buffer(buf, det, sum_e_pert) diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f index 48f0b9e5..816253c5 100644 --- a/src/cis/cis.irp.f +++ b/src/cis/cis.irp.f @@ -57,7 +57,11 @@ subroutine run implicit none integer :: i - call H_apply_cis + if(pseudo_sym)then + call H_apply_cis_sym + else + call H_apply_cis + endif print *, 'N_det = ', N_det print*,'******************************' print *, 'Energies of the states:' diff --git a/src/cis/h_apply.irp.f b/src/cis/h_apply.irp.f index 2505123d..95aa52e4 100644 --- a/src/cis/h_apply.irp.f +++ b/src/cis/h_apply.irp.f @@ -5,5 +5,10 @@ BEGIN_SHELL [ /usr/bin/env python2 ] from generate_h_apply import H_apply H = H_apply("cis",do_double_exc=False) print H + +H = H_apply("cis_sym",do_double_exc=False) +H.filter_only_connected_to_hf() +print H + END_SHELL diff --git a/src/cisd/cisd.irp.f b/src/cisd/cisd.irp.f index 9631f971..65f943d3 100644 --- a/src/cisd/cisd.irp.f +++ b/src/cisd/cisd.irp.f @@ -53,7 +53,11 @@ subroutine run implicit none integer :: i - call H_apply_cisd + if(pseudo_sym)then + call H_apply_cisd_sym + else + call H_apply_cisd + endif print *, 'N_det = ', N_det print*,'******************************' print *, 'Energies of the states:' diff --git a/src/cisd/h_apply.irp.f b/src/cisd/h_apply.irp.f index 1c864607..b18f6373 100644 --- a/src/cisd/h_apply.irp.f +++ b/src/cisd/h_apply.irp.f @@ -5,5 +5,9 @@ BEGIN_SHELL [ /usr/bin/env python2 ] from generate_h_apply import H_apply H = H_apply("cisd",do_double_exc=True) print H + +H = H_apply("cisd_sym",do_double_exc=True) +H.filter_only_connected_to_hf() +print H END_SHELL diff --git a/src/density_for_dft/density_for_dft.irp.f b/src/density_for_dft/density_for_dft.irp.f index 2c6aa46e..4514f111 100644 --- a/src/density_for_dft/density_for_dft.irp.f +++ b/src/density_for_dft/density_for_dft.irp.f @@ -9,6 +9,8 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s one_e_dm_mo_alpha_for_dft = data_one_e_dm_alpha_mo + damping_for_rs_dft * delta_alpha else if (density_for_dft .EQ. "input_density")then one_e_dm_mo_alpha_for_dft = data_one_e_dm_alpha_mo + else if (density_for_dft .EQ. "input_density_ao")then + call ao_to_mo(data_one_e_dm_alpha_mo,size(data_one_e_dm_alpha_mo,1),one_e_dm_mo_alpha_for_dft,size(one_e_dm_mo_alpha_for_dft,1)) else if (density_for_dft .EQ. "WFT")then provide mo_coef one_e_dm_mo_alpha_for_dft = one_e_dm_mo_alpha @@ -58,6 +60,8 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo + damping_for_rs_dft * delta_beta else if (density_for_dft .EQ. "input_density")then one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo + else if (density_for_dft .EQ. "input_density_ao")then + call ao_to_mo(data_one_e_dm_beta_mo,size(data_one_e_dm_beta_mo,1),one_e_dm_mo_beta_for_dft,size(one_e_dm_mo_beta_for_dft,1)) else if (density_for_dft .EQ. "WFT")then provide mo_coef one_e_dm_mo_beta_for_dft = one_e_dm_mo_beta @@ -119,16 +123,22 @@ END_PROVIDER one_e_dm_alpha_ao_for_dft = 0.d0 one_e_dm_beta_ao_for_dft = 0.d0 - do istate = 1, N_states - call mo_to_ao_no_overlap( one_e_dm_mo_alpha_for_dft(1,1,istate), & - size(one_e_dm_mo_alpha_for_dft,1), & - one_e_dm_alpha_ao_for_dft(1,1,istate), & - size(one_e_dm_alpha_ao_for_dft,1) ) - call mo_to_ao_no_overlap( one_e_dm_mo_beta_for_dft(1,1,istate), & - size(one_e_dm_mo_beta_for_dft,1), & - one_e_dm_beta_ao_for_dft(1,1,istate), & - size(one_e_dm_beta_ao_for_dft,1) ) - enddo + + if (density_for_dft .EQ. "input_density_ao")then + one_e_dm_alpha_ao_for_dft = data_one_e_dm_alpha_ao + one_e_dm_beta_ao_for_dft = data_one_e_dm_beta_ao + else + do istate = 1, N_states + call mo_to_ao_no_overlap( one_e_dm_mo_alpha_for_dft(1,1,istate), & + size(one_e_dm_mo_alpha_for_dft,1), & + one_e_dm_alpha_ao_for_dft(1,1,istate), & + size(one_e_dm_alpha_ao_for_dft,1) ) + call mo_to_ao_no_overlap( one_e_dm_mo_beta_for_dft(1,1,istate), & + size(one_e_dm_mo_beta_for_dft,1), & + one_e_dm_beta_ao_for_dft(1,1,istate), & + size(one_e_dm_beta_ao_for_dft,1) ) + enddo + endif END_PROVIDER diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index 12a2c2d1..0b289e99 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -89,3 +89,16 @@ doc: Weight of the states in state-average calculations. interface: ezfio size: (determinants.n_states) + +[thresh_sym] +type: Threshold +doc: Thresholds to check if a determinant is connected with HF +interface: ezfio,provider,ocaml +default: 1.e-15 + +[pseudo_sym] +type: logical +doc: If |true|, discard any Slater determinants with an interaction smaller than thresh_sym with HF. +interface: ezfio,provider,ocaml +default: False + diff --git a/src/determinants/h_apply.template.f b/src/determinants/h_apply.template.f index 169bd47e..f16d2e7f 100644 --- a/src/determinants/h_apply.template.f +++ b/src/determinants/h_apply.template.f @@ -150,6 +150,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl logical :: is_a_2h1p logical :: is_a_2h logical :: b_cycle + logical :: yes_no check_double_excitation = .True. iproc = iproc_in @@ -284,6 +285,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl $only_1h_double $only_1p_double $only_2h1p_double + $filter_only_connected_to_hf_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -339,6 +341,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl $only_1h_double $only_1p_double $only_2h1p_double + $filter_only_connected_to_hf_double key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -412,6 +415,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato logical :: is_a_1h logical :: is_a_1p logical :: is_a_2p + logical :: yes_no do k=1,N_int key_mask(k,1) = 0_bit_kind @@ -493,6 +497,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato $filter_only_1h1p_single $filter_only_1h2p_single $filter_only_2h2p_single + $filter_only_connected_to_hf_single key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = hole(k,1) diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 3f0f9335..6b164816 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -2257,3 +2257,38 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij) end +subroutine connected_to_hf(key_i,yes_no) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: key_i(N_int,2) + logical , intent(out) :: yes_no + double precision :: hij,phase + integer :: exc(0:2,2,2) + integer :: degree + integer :: m,p + yes_no = .True. + call get_excitation_degree(ref_bitmask,key_i,degree,N_int) + if(degree == 2)then + call i_H_j(ref_bitmask,key_i,N_int,hij) + if(dabs(hij) .lt. thresh_sym)then + yes_no = .False. + endif + else if(degree == 1)then + call get_single_excitation(ref_bitmask,key_i,exc,phase,N_int) + ! Single alpha + if (exc(0,1,1) == 1) then + m = exc(1,1,1) + p = exc(1,2,1) + ! Single beta + else + m = exc(1,1,2) + p = exc(1,2,2) + endif + hij = mo_one_e_integrals(m,p) + if(dabs(hij) .lt. thresh_sym)then + yes_no = .False. + endif + else if(degree == 0)then + yes_no = .True. + endif +end diff --git a/src/tools/rotate_mos.irp.f b/src/tools/rotate_mos.irp.f new file mode 100644 index 00000000..3e79c4f7 --- /dev/null +++ b/src/tools/rotate_mos.irp.f @@ -0,0 +1,18 @@ +program rotate_mos + implicit none + integer :: iorb,jorb + read(5,*)iorb,jorb + double precision, allocatable :: mo_coef_tmp(:,:) + allocate(mo_coef_tmp(ao_num,mo_num)) + mo_coef_tmp = mo_coef + integer :: i,j + double precision :: dsqrt2_inv + dsqrt2_inv = 1.d0/dsqrt(2.d0) + do i = 1, ao_num + mo_coef(i,iorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) + mo_coef_tmp(i,jorb) ) + mo_coef(i,jorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) - mo_coef_tmp(i,jorb) ) + enddo + touch mo_coef + call save_mos + +end diff --git a/src/tools/save_one_e_dm.irp.f b/src/tools/save_one_e_dm.irp.f index e850131e..c888f55c 100644 --- a/src/tools/save_one_e_dm.irp.f +++ b/src/tools/save_one_e_dm.irp.f @@ -1,15 +1,15 @@ program save_one_e_dm implicit none BEGIN_DOC -! Program that computes the one body density on the |MO| basis +! Program that computes the one body density on the |MO| and |AO| basis ! for $\alpha$ and $\beta$ electrons from the wave function ! stored in the |EZFIO| directory, and then saves it into the ! :ref:`module_aux_quantities`. ! ! Then, the global variable :option:`aux_quantities data_one_e_dm_alpha_mo` -! and :option:`aux_quantities data_one_e_dm_beta_mo` will automatically -! read this density in the next calculation. This can be used to perform -! damping on the density in |RSDFT| calculations (see +! and :option:`aux_quantities data_one_e_dm_beta_mo` (and the corresponding for |AO|) +! will automatically ! read this density in the next calculation. +! This can be used to perform damping on the density in |RSDFT| calculations (see ! :ref:`module_density_for_dft`). END_DOC read_wf = .True. @@ -25,4 +25,6 @@ subroutine routine_save_one_e_dm END_DOC call ezfio_set_aux_quantities_data_one_e_dm_alpha_mo(one_e_dm_mo_alpha) call ezfio_set_aux_quantities_data_one_e_dm_beta_mo(one_e_dm_mo_beta) + call ezfio_set_aux_quantities_data_one_e_dm_alpha_ao(one_e_dm_ao_alpha) + call ezfio_set_aux_quantities_data_one_e_dm_beta_ao(one_e_dm_ao_beta) end